structure {{group.name}} = struct structure D = Dynamics structure SEventPriority = struct type priority = real fun compare (x,y) = Real.compare (x,y) type item = real * int fun priority (x : item) = (#1(x)) end structure SEQ = TEventQueue (structure P = SEventPriority type value = int fun value (x : P.item) = (#2(x))) {% if (trace) %} val trace = true {% else %} val trace = false {% endif %} fun putStrLn out str = (TextIO.output (out, str); TextIO.output (out, "\n")) fun putStr out str = (TextIO.output (out, str)) fun showBoolean b = (if b then "1" else "0") fun showReal n = let open StringCvt open Real in (if n < 0.0 then "-" else "") ^ (fmt (FIX (SOME 12)) (abs n)) end fun showRealArray v = (String.concatWith ", " (Array.foldr (fn (x, ax) => (showReal x)::ax) [] v)) fun showBoolArray v = (String.concatWith ", " (Array.foldr (fn (x, ax) => (showBoolean x)::ax) [] v)) fun foldl1 f lst = let val v = List.hd lst val lst' = List.tl lst in List.foldl f v lst' end fun dimVals [m,n] = (m,n) | dimVals _ = raise BitSparseMatrix.Shape fun fromDiag (m, n, offset, a, dflt) = if Index.validShape [m,n] then (let val na = BitTensor.Array.length a val na' = na-1 val te = BitTensor.new ([m,n], dflt) fun diag (i, j, ia) = let val v = BitTensor.Array.sub (a, ia) val ia' = (BitTensor.update (te, [i,j], v); if ia = na' then 0 else ia+1) in if (i (i := (!i)+1; if (!i) <= n then update (sample, !i, x) else (let val k = randomInt (0, !i) in if k < n then update(sample, randomInt (0, n), x) else () end))) end fun index_foldri (a, b, f, init) = if a <= b then index_foldri (a, b-1, f, f (b, init)) else init fun assoc ((k,v)::rst) x = if k=x then SOME v else assoc rst x | assoc [] x = NONE fun member (k,lst) = case List.find (fn(x)=>k=x) lst of SOME _ => true | NONE => false {% if CheckBounds %} val getindex = Array.sub {% else %} val getindex = Unsafe.Array.sub {% endif %} val RandomInitFromInt = RandomMTZig.fromInt val RandomInit = RandomMTZig.fromEntropy val ZigInit = RandomMTZig.ztnew fun optApply fopt args = case fopt of SOME f => SOME (f args) | NONE => NONE fun optApplyF fopt fargs = case fopt of SOME f => SOME (fn(x) => f (fargs x)) | NONE => NONE fun make_initial_vector (make_real_state,make_bool_state,n,nev,ndsc,nregime,h0,rs) (i,condfun,initfun,SOME evinitial,dinitial,SOME rinitial,extinitial,extevinitial) = let val initial = initfun rs val dinitial': real array -> real array = case dinitial of SOME v => v rs | NONE => (fn v => v) val err0 = case (!(D.abstol)) of SOME tol => tol | NONE => 0.0 in D.condApply condfun (D.RegimeState (0.0, 0.0, initial(make_real_state n), evinitial (make_real_state nev), dinitial' (make_real_state ndsc), rinitial (make_bool_state nregime), extinitial(), extevinitial(), make_real_state n, make_real_state n, make_real_state nev, D.Right {h=h0, cst=h0, r=err0, tevent=Real.posInf}, D.RootBefore)) end | make_initial_vector (make_real_state,make_bool_state,n,nev,ndsc,nregime,h0,rs) (i,condfun,initfun,SOME evinitial,NONE,NONE,extinitial,extevinitial) = let val initial = initfun rs val err0 = case (!(D.abstol)) of SOME tol => tol | NONE => 0.0 in D.condApply condfun (D.EventState (0.0, 0.0, (initfun rs)(make_real_state n), evinitial (make_real_state nev), extinitial (), extevinitial (), initial(make_real_state n), initial(make_real_state n), make_real_state nev, D.Right {h=h0, cst=h0, r=err0, tevent=Real.posInf}, D.RootBefore)) end | make_initial_vector (make_real_state,make_bool_state,n,nev,ndsc,nregime,h0,rs) (i,condfun,initfun,NONE,NONE,NONE,extinitial,extevinitial) = let val initial = initfun rs val err0 = case (!(D.abstol)) of SOME tol => tol | NONE => 0.0 in D.ContState (0.0, 0.0, (initfun rs)(make_real_state n), extinitial(), extevinitial(), make_real_state n, D.Right {h=h0, cst=h0, r=err0, tevent=Real.posInf}) end | make_initial_vector (make_real_state,make_bool_state,n,nev,ndsc,nregime,h0,rs) (i,condfun,initial,_,_,_,extinitial,extevinitial) = (raise Fail "make_initial_vector: invalid domain") fun state_sub (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root), i) = getindex(y, i) | state_sub (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root), i) = getindex(y, i) | state_sub (D.ContState (x, cx, y, ext, extev, ynext, _), i) = getindex(y, i) fun state_y (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = y | state_y (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = y | state_y (D.ContState (x, cx, y, ext, extev, ynext, _)) = y fun state_indep (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = x | state_indep (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = x | state_indep (D.ContState (x, cx, y, ext, extev, ynext, _)) = x fun state_ext (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = ext | state_ext (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = ext | state_ext (D.ContState (x, cx, y, ext, extev, ynext, _)) = ext fun state_extev (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = extev | state_extev (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = extev | state_extev (D.ContState (x, cx, y, ext, extev, ynext, _)) = extev fun state_regime (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = r | state_regime (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = Array.fromList [] | state_regime (D.ContState (x, cx, y, ext, extev, ynext, _)) = Array.fromList[] fun state_event (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = ev | state_event (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = ev | state_event (D.ContState (x, cx, y, ext, extev, ynext, _)) = Array.fromList [] fun state_discrete (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = d | state_discrete (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = Array.fromList [] | state_discrete (D.ContState (x, cx, y, ext, extev, ynext, _)) = Array.fromList [] fun ext_sub (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root), i) = getindex(ext, i) | ext_sub (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root), i) = getindex(ext, i) | ext_sub (D.ContState (x, cx, y, ext, extev, ynext, _), i) = getindex(ext, i) fun event_sub (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root), i) = getindex(ev, i) | event_sub (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root), i) = getindex(ev, i) | event_sub (D.ContState (x, cx, y, ext, extev, ynext, _), i) = raise Fail "event_sub: ContState does not have event information" fun state_root (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, _, root)) = root | state_root (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, _, root)) = root | state_root (D.ContState (x, cx, y, ext, extev, ynext, _)) = D.RootBefore fun state_h (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, cst, root)) = D.controller_h cst | state_h (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, cst, root)) = D.controller_h cst | state_h (D.ContState (x, cx, y, ext, extev, ynext, cst)) = D.controller_h cst fun update_state_h (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, cst, root), h') = D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, D.controller_update_h(cst,h'), root) | update_state_h (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, cst, root), h') = D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, D.controller_update_h(cst,h'), root) | update_state_h (D.ContState (x, cx, y, ext, extev, ynext, cst), h') = D.ContState (x, cx, y, ext, extev, ynext, D.controller_update_h(cst,h')) fun copy_array a = let val len = Array.length a val b = if len>0 then Array.array (len, Array.sub(a,0)) else a in if len>0 then Array.copy {di=0, dst=b, src=a} else (); b end fun copy_state (D.RegimeState (x, cx, y, ev, d, r, ext, extev, ynext, yrsp, enext, cst, root)) = D.RegimeState(x, cx, copy_array y, copy_array ev, copy_array d, copy_array r, copy_array ext, copy_array extev, ynext, yrsp, enext, cst, root) | copy_state (D.EventState (x, cx, y, ev, ext, extev, ynext, yrsp, enext, cst, root)) = D.EventState (x, cx, copy_array y, copy_array ev, copy_array ext, copy_array extev, ynext, yrsp, enext, cst, root) | copy_state (D.ContState (x, cx, y, ext, extev, ynext, cst)) = D.ContState (x, cx, copy_array y, copy_array ext, copy_array extev, ynext, cst) fun random_normal (rs,zt) = RandomMTZig.randNormal(rs,zt) fun random_uniform (rs) = RandomMTZig.randUniform(rs) fun random_uniform_vector (rs) = RandomMTZig.randUniform(rs) fun sample_without_replacement (rs, (* random state *) populationSize, (* size of set sampling from *) sampleSize, (* size of sample *) f) = Loop.foldsi (0, sampleSize, fn (m, t) => (* t: total inputs processed *) let val u = random_uniform rs in if Real.>= (Real.*(Real.fromInt (populationSize - t), u), Real.fromInt (sampleSize - m)) then Loop.SeqRecur (t+1) else (f(m, t); Loop.SeqAdvance (t+1)) end, 0) val global_h = ref 0.1 val label = "{{group.name}}" val N = {{group.order}} (* total population size *) {% if not (group.properties == []) %} {% for p in dict (group.properties) %} val {{p.name}} = {{p.value.exprML}} {% endfor %} {% endif %} (* network propagation delays for each projection *) val D: real list = [{% for pr in dict (group.projections) %}{{pr.value.delay.exprML}}{% if not loop.last %},{% endif %}{% endfor %}] val minDelay = List.foldl Real.min (hd D) (tl D) {% for pop in dict (group.populations) %} val N_{{pop.name}} = {{ pop.value.size }} val {{pop.name}}_n0 = {{pop.value.start}} val {{pop.name}}_parameters = {{pop.value.prototype.name}}.paramfun() val {{pop.name}}_initial_vector = let val parameters = {{pop.name}}_parameters val field_vector = Vector.tabulate (N_{{pop.name}}, fn (i) => {{pop.value.prototype.name}}.fieldfun()) val make_initial_vector' = make_initial_vector ({{pop.value.prototype.name}}.make_real_state,{{pop.value.prototype.name}}.make_bool_state, {{pop.value.prototype.name}}.n,{{pop.value.prototype.name}}.nev,{{pop.value.prototype.name}}.ndsc, {{pop.value.prototype.name}}.nregime, Real.min((!global_h), minDelay), RandomInitFromInt ({{pop.name}}_n0) ) in List.tabulate (N_{{pop.name}}, fn(i) => let val fields = Vector.sub(field_vector,i) val cond = optApply {{pop.value.prototype.name}}.condfun (parameters,fields) val init = fn(rs) => {{pop.value.prototype.name}}.initfun(parameters,fields,rs,ZigInit()) val evinitial = {{pop.value.prototype.name}}.initcondfun val dinitial = optApplyF {{pop.value.prototype.name}}.dinitfun (fn(rs) => (parameters,fields,rs,ZigInit())) val rinitial = {{pop.value.prototype.name}}.initregfun val extinitial = {{pop.value.prototype.name}}.initextfun (parameters,fields) val extevinitial = {{pop.value.prototype.name}}.initextevfun (parameters,fields) in (i, fields, make_initial_vector' (i,cond,init,evinitial,dinitial,rinitial,extinitial,extevinitial)) end) end val {{pop.name}}_ext_ports = [{% for portspec in dict (group.projectionPorts) %}{% if portspec.name == pop.name %}{% for port in portspec.value %}{{port}}{% if not loop.last %},{% endif %}{% endfor %}{% endif %}{% endfor %}] fun {{pop.name}}_f (i,p,fld,rs_vector,rszt_vector) = let val rs = Vector.sub(rs_vector, i) val rszt = Vector.sub(rszt_vector, i) in D.integral({{pop.value.prototype.name}}.odefun (p,fld,rs,rszt), {{pop.value.prototype.name}}.interpfun, optApply {{pop.value.prototype.name}}.condfun (p,fld), optApply {{pop.value.prototype.name}}.posfun (p,fld,rs,rszt), optApply {{pop.value.prototype.name}}.negfun (p,fld,rs,rszt), optApply {{pop.value.prototype.name}}.dposfun (p,fld,rs,rszt), {{pop.value.prototype.name}}.regfun) end (* {% if pop.value.prototype.fieldExprML %} val fieldV = getindex ({{pop.name}}_field_vector,i) {% endif %} *) fun {{pop.name}}_run f (rs,rszt) (statesample,extsample,evsample) n0 = let {% if CheckBounds %} val update = Array.update {% else %} val update = Unsafe.Array.update {% endif %} fun f'(i,fld) = f (i,{{pop.name}}_parameters,fld,rs,rszt) val ext_buf = Array.array({{pop.value.prototype.name}}.next, 0.0) fun merge_extevent xs = case xs of [] => xs | [(tx, wx)] => xs | ((tx, wx) :: (ty, wy) :: rst) => let val diff = Real.-(tx, ty) in if Real.<= (Real.abs(diff), Real.*(2.0,D.float_eps)) then merge_extevent ((tx, wx+wy) :: rst) else (if Real.<(diff, 0.0) then (tx, wx) :: (merge_extevent ((ty, wy) :: rst)) else (ty, wy) :: (merge_extevent ((tx, wx) :: rst))) end fun show_extev_map extev_map = List.app (fn(index,lst) => putStrLn TextIO.stdOut ("external event: index=" ^ (Int.toString index) ^ ": " ^ (String.concat (map (fn(t,v) => "[" ^ (showReal t) ^ ": " ^ (showReal v) ^ "]") lst)))) extev_map fun extevent_clear_at (extev_i,ext_i,tstart) = (List.foldl (fn(port,index) => if getindex(extev_i, index) <= tstart then (update (ext_i, index, 0.0); update (extev_i, index, Real.posInf); index+1) else index+1) 0 {{pop.name}}_ext_ports) before (ignore ({{pop.value.prototype.name}}.linkextevfun(extev_i,extev_i))) fun extevent loc extev_map (i, extev_i, ext_i, tstart, tend, iter) = (List.foldl (fn ((index,lst), ax) => (if trace then (putStrLn TextIO.stdOut ("extevent " ^ loc ^ ": cell " ^ (Int.toString (i+n0)) ^ ": " ^ (Int.toString (List.length lst)) ^ " tstart = " ^ (showReal tstart) ^ " tend = " ^ (showReal tend)); List.app (fn(t,v) => putStr TextIO.stdOut (" " ^ (showReal t) ^ " : " ^ (showReal v))) lst; putStrLn TextIO.stdOut "") else (); let fun recur lst = case lst of (t',v)::rst => (if trace then putStrLn TextIO.stdOut ("extevent " ^ loc ^ ": cell " ^ (Int.toString (i+n0)) ^ ": " ^ " t' = " ^ (Real.toString t') ^ " tprev = " ^ (Real.toString (getindex(extev_i, index)))) else (); if tstart-t' > D.float_eps then raise Fail ("extevent " ^ loc ^ ": cell " ^ (Int.toString (i+n0)) ^ ": event time out of range: tstart=" ^ (showReal tstart) ^ " tprev=" ^ (showReal (getindex(extev_i, index))) ^ " t'=" ^ (showReal t') ^ " tstart-t' = " ^ (showReal (tstart - t'))) else (if t' <= tend andalso (let val tprev = getindex(extev_i, index) in ((not (Real.isFinite(tprev))) orelse (tprev < tstart)) end) then (if trace then putStrLn TextIO.stdOut ("extevent " ^ loc ^ ": cell " ^ (Int.toString (i+n0)) ^ ": " ^ " updating event with t' = " ^ (Real.toString t')) else (); update (ext_buf, index, v); update (extev_i, index, t'); if List.null rst then ax else (case rst of (t'',v')::rst' => (if t'' <= t' then raise Fail ("extevent " ^ loc ^ ": cell " ^ (Int.toString (i+n0)) ^ ": successive events with incorrect order: t'=" ^ (showReal t') ^ " t''=" ^ (showReal t'')) else ()) | _ => (); (index,rst)::ax)) else (if List.null lst then ax else (index,lst)::ax))) | [] => if getindex(extev_i, index) < tstart then (update (ext_i, index, 0.0); update (extev_i, index, Real.posInf); ax) else ax in recur lst end)) [] extev_map) before (ignore ({{pop.value.prototype.name}}.linkextevfun(extev_i,extev_i))) fun recur (i,fld,nstate,extev_map,spikes,statelog: (int * D.model_state) list,iter,t_next) = let val t = state_indep nstate val h = state_h nstate val ext_i = state_ext nstate val extev_i = state_extev nstate val spikes' = case state_root nstate of D.RootFound _ => if (Real.>= (event_sub(nstate,0), 0.0)) then (i+n0,state_indep nstate,1.0)::spikes else spikes | _ => spikes val statelog' = case state_root nstate of D.RootFound _ => (if (i (if (i statelog val _ = if trace then (putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": recur: " ^ " t = " ^ (showReal (state_indep nstate)) ^ " t_next = " ^ (showReal t_next) ^ " root state is " ^ (D.showRoot (state_root nstate)) ^ " discrete state is " ^ (showRealArray (state_discrete nstate)) ^ " regime state is " ^ (showBoolArray (state_regime nstate)) ^ " h is " ^ (showReal (state_h nstate)) ^ " extev is " ^ (showRealArray (state_extev nstate)) ); show_extev_map extev_map) else () val _ = if (state_h nstate) < D.float_eps then raise Fail ("cell " ^ (Int.toString (i+n0)) ^ ": time step too small (" ^ (Real.toString (state_h nstate)) ^ ")") else () in case state_root nstate of D.RootFound _ => (if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": root state is " ^ (D.showRoot (state_root nstate))) else (); let val _ = Array.modifyi (fn(ei,_) => if abs(getindex(extev_i, ei)-t) <= D.float_eps then getindex(ext_buf, ei) else 0.0) ext_i val nstate' = f' (i,fld) nstate handle Fail s => raise Fail ("exception in function evaluation for cell " ^ (Int.toString (i+n0)) ^ ": " ^ s) in recur (i, fld, nstate', extev_map, spikes', statelog', iter+1, t_next) end) | D.RootStep [] => (if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": root state is " ^ (D.showRoot (state_root nstate))) else (); if List.null extev_map then (nstate, spikes', statelog') else raise Fail ("RootStep: cell " ^ (Int.toString (i+n0)) ^ " has a non-empty external event map at t = " ^ (showReal (state_indep nstate)) ^ " t_next = " ^ (showReal t_next))) | D.RootStep (h::hs) => (if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": root state is " ^ (D.showRoot (state_root nstate))) else (); let val t = state_indep nstate val extev_map' = extevent "RootStep" extev_map (i, extev_i, ext_i, t, t_next, iter) val nstate' = f' (i,fld) nstate handle Fail s => raise Fail ("exception in function evaluation for cell " ^ (Int.toString (i+n0)) ^ ": " ^ s) in recur (i, fld, nstate', extev_map', spikes', statelog', iter+1, t_next) end) | D.RootBefore => (if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": root state is " ^ (D.showRoot (state_root nstate)) ^ " t = " ^ (showReal (state_indep nstate)) ^ " t_next = " ^ (showReal t_next) ^ " t_next - t = " ^ (showReal (t_next - (state_indep nstate))) ^ " t_next - t > D.float_eps = " ^ (showBoolean ((t_next - (state_indep nstate)) > D.float_eps))) else (); if (t_next - (state_indep nstate)) > D.float_eps then (let val t_next_event = if List.null extev_map then t_next else (foldl (fn((i,lst),ax) => Real.min(#1(hd lst), ax)) t_next extev_map) val t = state_indep nstate val h = state_h nstate val nstate' = if ((h+t)-t_next_event) > D.float_eps then update_state_h (nstate,Real.max(t_next-t,D.float_eps)) else nstate val _ = extevent_clear_at (extev_i,ext_i,t) val nstate'' = f' (i,fld) nstate' handle Fail s => raise Fail ("exception in function evaluation for cell " ^ (Int.toString (i+n0)) ^ ": " ^ s) in recur (i, fld, nstate'', extev_map, spikes', statelog', iter+1, t_next) end) else (let val t = state_indep nstate val _ = extevent_clear_at (extev_i,ext_i,t) in if List.null extev_map then (nstate, spikes', statelog') else (show_extev_map extev_map; raise Fail ("RootBefore: cell " ^ (Int.toString (i+n0)) ^ " has a non-empty external event map at t = " ^ (showReal (state_indep nstate)) ^ " t_next = " ^ (showReal t_next))) end)) | D.RootAfter (ei,_) => (if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": root state is " ^ (D.showRoot (state_root nstate)) ^ " t = " ^ (showReal (state_indep nstate))) else (); let val t = state_indep nstate val _ = extevent_clear_at (extev_i,ext_i,t) val extev_map' = extevent "RootAfter" extev_map (i, extev_i, ext_i, t, t_next, iter) val _ = if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": root state is " ^ (D.showRoot (state_root nstate)) ^ " t = " ^ (showReal (state_indep nstate)) ^ " y = " ^ (showReal (getindex (state_y nstate, 0))) ^ " e = " ^ (showReal (getindex (state_event nstate, 0))) ^ " d = " ^ (showRealArray (state_discrete nstate)) ^ " r = " ^ (showBoolArray (state_regime nstate)) ^ " extev = " ^ (showRealArray (state_extev nstate))) else () val nstate' = f' (i,fld) nstate handle Fail s => raise Fail ("exception in function evaluation for cell " ^ (Int.toString (i+n0)) ^ ": " ^ s) val _ = if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": (RootAfter) root state is " ^ (D.showRoot (state_root nstate')) ^ " t = " ^ (showReal (state_indep nstate')) ^ " y = " ^ (showReal (getindex (state_y nstate', 0))) ^ " e = " ^ (showReal (getindex (state_event nstate', 0))) ^ " e > float_eps = " ^ (showBoolean (abs(getindex (state_event nstate', 0)) > D.float_eps)) ^ " d = " ^ (showRealArray (state_discrete nstate)) ^ " r = " ^ (showBoolArray (state_regime nstate)) ^ " extev = " ^ (showRealArray (state_extev nstate'))) else () in recur (i, fld, nstate', extev_map', spikes', statelog', iter+1, t_next) end) end in fn (Wmap,t_next) => fn ((i,fld,input),(nstates,spikes,statelog)) => let val t = state_indep input val h = state_h input val ext_i = state_ext input val extev_i = state_extev input val extev_map = case IntMap.find (Wmap, i+n0) of SOME spm => (if trace then (putStrLn TextIO.stdOut ("extev_map: cell " ^ (Int.toString (i+n0)) ^ " external events: "); IntMap.appi (fn (port, lst) => (putStr TextIO.stdOut (" port " ^ (Int.toString port) ^ ": "); List.app (fn(t,v) => putStr TextIO.stdOut ("[" ^ (showReal t) ^ " : " ^ (showReal v) ^ "]")) lst; putStrLn TextIO.stdOut "") ) spm) else (); #2(List.foldl (fn(port,(index,ax)) => case IntMap.find (spm, port) of SOME lst => (index+1, (index,merge_extevent lst)::ax) | NONE => (index+1, ax)) (0,[]) {{pop.name}}_ext_ports)) | NONE => [] val statelog' = if (i Real.min(#1(hd lst), ax)) t_next extev_map) val input' = if ((h+t)-t_next) > D.float_eps then update_state_h (input,Real.max(t_next_event-t,D.float_eps)) else input val _ = if trace then putStrLn TextIO.stdOut ("cell " ^ (Int.toString (i+n0)) ^ ": run: " ^ "t = " ^ (showReal t) ^ " t_next = " ^ (showReal t_next) ^ " h = " ^ (showReal h) ^ " root state is " ^ (D.showRoot (state_root input))) else () val (nstate,spikes',statelog') = recur (i, fld, input', extev_map, spikes, statelog', 0, t_next) in ((i,fld,nstate)::nstates, spikes', statelog') end end {% endfor %} val popmap = [ {% for s in dict (group.sets) %} ("{{s.name}}"{% for p in s.value.populations %}::"{{p.name}}"{% endfor %}::[]){% if not loop.last %},{% endif %} {% endfor %} ] {% if group.conntypes %}{% for conn in dict (group.conntypes) %} {% if conn.value.sysFn %} val {{conn.name}}_initial = {{conn.value.initialExprML}} val {{conn.name}}_f = Model_{{conn.name}}.{{conn.value.sysFn}} {% endif %} {% endfor %}{% endif %} {% if group.projections %}{% for prj in dict (group.projections) %} val {{prj.name}}_rs = ref (RandomInitFromInt({{loop.index}})) {% endfor %}{% endif %} val initial = ( {% for pop in dict (group.populations) %} {{pop.name}}_initial_vector{% if not loop.last %},{% endif %} {% endfor %} ) fun update_projection_rs (seeds) = ( {% for prj in dict (group.projections) %} (let val seed = List.nth (seeds, {{loop.index0}}) in {{prj.name}}_rs := RandomInitFromInt(seed) end){% if not loop.last %};{% endif %} {% endfor %} ) val rszt_vector = ( {% for pop in dict (group.populations) %} Vector.tabulate (N_{{pop.name}}, fn (i) => ZigInit()){% if not loop.last %},{% endif %} {% endfor %} ) fun make_rs_vector NONE = ( {% for pop in dict (group.populations) %} Vector.tabulate (N_{{pop.name}}, fn (i) => RandomInitFromInt(i+{{pop.name}}_n0)){% if not loop.last %},{% endif %} {% endfor %} ) | make_rs_vector (SOME seeds) = ( {% for pop in dict (group.populations) %} (let val seed = List.nth (seeds, {{loop.index0}}) in Vector.tabulate (N_{{pop.name}}, fn (i) => RandomInitFromInt(i+seed)) end){% if not loop.last %},{% endif %} {% endfor %} ) val P = {{length(dict(group.populations))}} val C = {{length(dict(group.conntypes))}} val Pn = [ {% for pop in dict (group.populations) %} {{pop.name}}_n0{% if not loop.last %},{% endif %} {% endfor %} ] val Plabels = [ {% for pop in dict (group.populations) %} "{{pop.name}}"{% if not loop.last %},{% endif %} {% endfor %} ] fun frun (statesample,extsample,evsample) = fn (I,t_next, ( {% for pop in dict (group.populations) %} {{pop.name}}_state_vector{% if not loop.last %},{% endif %} {% endfor %}), ( {% for pop in dict (group.populations) %} {{pop.name}}_rs_vector{% if not loop.last %},{% endif %} {% endfor %}), ( {% for pop in dict (group.populations) %} {{pop.name}}_rszt_vector{% if not loop.last %},{% endif %} {% endfor %}) ) => let {% for pop in dict (group.populations) %} val {{pop.name}}_run' = {{pop.name}}_run {{pop.name}}_f ({{pop.name}}_rs_vector,{{pop.name}}_rszt_vector) (statesample,extsample,evsample) {{pop.name}}_n0 {% endfor %} {% for pop in dict (group.populations) %} val ({{pop.name}}_states,{{pop.name}}_spikes,{{pop.name}}_statelog) = List.foldl ({{pop.name}}_run' (I,t_next)) ([],[],[]) {{pop.name}}_state_vector {% endfor %} in (( {% for pop in dict (group.populations) %} {{pop.name}}_states{% if not loop.last %},{% endif %} {% endfor %} ), ( {% for pop in dict (group.populations) %} {{pop.name}}_spikes{% if not loop.last %},{% endif %} {% endfor %} ), [ {% for pop in dict (group.populations) %} {{pop.name}}_statelog{% if not loop.last %},{% endif %} {% endfor %} ]) end fun ftime ( {% for pop in dict (group.populations) %} {{pop.name}}_state_vector: (int * real array * D.model_state) list{% if not loop.last %},{% endif %} {% endfor %} ) = {% with pop = first (dict (group.populations)) %} state_indep (#3 (hd {{pop.name}}_state_vector)) {% endwith %} fun fspikeidxs spikesetName = let val spikesetPops = case List.find (fn(x) => spikesetName = (hd x)) popmap of SOME lst => tl lst | NONE => raise Fail ("invalid population set: " ^ spikesetName) val pops = [{% for pop in dict (group.populations) %}"{{pop.name}}"{% if not loop.last %},{% endif %}{% endfor %}] val popIdxs = ListPair.zip (pops, List.tabulate (List.length pops, fn (i) => i)) in List.map (valOf o assoc popIdxs) spikesetPops end fun fspikes spikelogIdxs ( {% for pop in dict (group.populations) %} {{pop.name}}_spike_i{% if not loop.last %},{% endif %} {% endfor %} ) = let fun enqueueSpikes xs = List.foldl (fn(xs,ax) => List.foldl (fn((i,t,nv),ax) => SEQ.addEvent ((t,i),ax)) ax xs) SEQ.empty xs val all_spike_i = [{% for pop in dict (group.populations) %}{{pop.name}}_spike_i{% if not loop.last %},{% endif %}{% endfor %}] val labeled_spike_i = ListPair.zip (List.tabulate (List.length all_spike_i, fn (i) => i), all_spike_i) val (log_spike_i,ext_spike_i) = List.partition (fn(i,x) => member (i, spikelogIdxs)) labeled_spike_i in (enqueueSpikes all_spike_i, enqueueSpikes (map #2 log_spike_i)) end {% macro cartesian_product(sp, tp) %} {% for s,t in allCombs(sp,tp) %} {{ caller(s,t) }}{% if not loop.last %},{% endif %} {% endfor %} {% endmacro %} {% macro for_each(name, sp, tp, component, cstate) %} val Pr_{{name}} = BitSparseMatrix.fromGeneratorList [N,N] [ {% call cartesian_product (sp,tp) %}{offset=[{{t.start}},{{s.start}}], fshape=[{{t.size}},{{s.size}}], f=(fn (i) => #{{cstate}}({{component}}_f {{component}}_initial) )}{% endcall %} ] {% endmacro %} {% macro all_to_all(name, sp, tp) %} val Pr_{{name}} = BitSparseMatrix.fromTensorList [N,N] [ {% call cartesian_product (sp,tp) %}{offset=[{{t.start}},{{s.start}}], tensor=(BitTensor.new ([{{t.size}},{{s.size}}],true)), sparse=false}{% endcall %} ] {% endmacro %} {% macro one_to_one(name, sp, tp) %} val Pr_{{name}} = BitSparseMatrix.fromTensorList [N,N] [ {% call cartesian_product (sp,tp) %}{offset=[{{t.start}},{{s.start}}], tensor=(fromDiag ({{t.size}},{{s.size}},{% if haskey(t,"relativeStart") %}{{t.relativeStart}}{% else %}0{% endif %},BitArray.fromList [true],false)), sparse=true}{% endcall %} ] {% endmacro %} {% macro from_file(name, sp, tp, filename) %} val Pr_{{name}} = let val infile = TextIO.openIn "{{filename}}" val S = TensorFile.realTensorRead (infile) val _ = TextIO.closeIn infile in SparseMatrix.fromTensorSliceList [N,N] [ {% with %} {% set soffset = 0 %} {% for s in sp %} {% set toffset = 0 %} {% for t in tp %} {offset=[{{t.start}},{{s.start}}], slice=(BitTensorSlice.slice ([([{{toffset}},{{soffset}}],[{{toffset}}+{{t.size}}-1,{{soffset}}+{{s.size}}-1])],S)), sparse=false}{% if not loop.last %},{% endif %} {% set toffset = toffset + t.size %} {% endfor %}{% if not loop.last %},{% endif %} {% set soffset = soffset + s.size %} {% endfor %} {% endwith %} ] end {% endmacro %} {% macro preorder(sp, tp) %} {% for s in sp %} {% for t in tp %} {{ caller(s,t) }}{% if not loop.last %},{% endif %} {% endfor %}{% if not loop.last %},{% endif %} {% endfor %} {% endmacro %} {% macro postorder(sp, tp) %} {% for t in tp %} {% for s in sp %} {{ caller(s,t) }}{% if not loop.last %},{% endif %} {% endfor %}{% if not loop.last %},{% endif %} {% endfor %} {% endmacro %} {% macro random_fan_in(name, sp, tp, prob) %} val Pr_{{name}} = let val f = (fn (i) => if (Real.> ({{prob}}, random_uniform (!{{name}}_rs))) then true else false) in BitSparseMatrix.fromGeneratorList [N,N] [ {% call postorder (sp,tp) %} {offset=[{{t.start}},{{s.start}}], fshape=[{{t.size}},{{s.size}}], f=f}{% endcall %} ] end {% endmacro %} {% macro random_fan_out(name, sp, tp, prob) %} val Pr_{{name}} = let val f = (fn (i) => if (Real.> ({{prob}}, random_uniform (!{{name}}_rs))) then true else false) in BitSparseMatrix.fromGeneratorList [N,N] [ {% call preorder (sp,tp) %} {offset=[{{t.start}},{{s.start}}], fshape=[{{t.size}},{{s.size}}], f=f}{% endcall %} ] end {% endmacro %} {% macro random_fan_in_num(name, sp, tp, num) %} val Pr_{{name}} = BitSparseMatrix.fromVectors [N,N] [ {% call postorder (sp,tp) %} (let val sample = Array.array (Int.min({{num | int}},{{s.size}}) * {{t.size}}, (0,0,false)) {% if CheckBounds %} val update = Array.update {% else %} val update = Unsafe.Array.update {% endif %} val _ = Loop.foldi (0,{{t.size}}, fn (i,k) => (sample_without_replacement (!{{name}}_rs,{{s.size}},{{num | int}}, fn (j,t) => (update (sample,k+j,(i,t,true)))); k+{{num | int}}), 0) val v = Array.vector sample in {offset=[{{t.start}},{{s.start}}], v=v, shape_v=[{{t.size}},{{s.size}}]} end){% endcall %} ] {% endmacro %} {% macro random_fan_out_num(name, sp, tp, num) %} val Pr_{{name}} = BitSparseMatrix.fromVectors [N,N] [ {% call preorder (sp,tp) %} (let val sample = Array.array (Int.min({{num | int}},{{s.size}}) * {{t.size}}, (0,0,false)) {% if CheckBounds %} val update = Array.update {% else %} val update = Unsafe.Array.update {% endif %} val _ = Loop.foldi (0,{{s.size}}, fn (i,k) => (sample_without_replacement (!{{name}}_rs,{{t.size}},{{num | int}}, fn (j,t) => (update (sample,k+j,(i,t,true)))); k+{{num | int}}), 0) val v = Array.vector sample in {offset=[{{t.start}},{{s.start}}], v=v, shape_v=[{{t.size}},{{s.size}}]} end){% endcall %} ] {% endmacro %} {% macro range_map(name, sp, tp) %} val srangemap_{{name}} = [ {% with %} {% set soffset = 0 %} {% for s in sp %} {size={{s.size}} localStart={{soffset}}, globalStart={{s.start}} }{% if not loop.last %},{% endif %} {% set soffset = soffset + s.size %} {% endfor %} {% endwith %} ] val trangemap_{{name}} = [ {% with %} {% set toffset = 0 %} {% for t in tp %} {size={{t.size}} localStart={{toffset}}, globalStart={{t.start}} }{% if not loop.last %},{% endif %} {% set toffset = toffset + t.size %} {% endfor %} {% endwith %} ] {% endmacro %} fun fprojection () = (let {% for pr in dict (group.projections) %} val _ = putStrLn TextIO.stdOut "constructing {{pr.name}}" {% if pr.value.connectivity.type.sysFn %} {% call for_each(pr.name, pr.value.source.populations, pr.value.destination.populations, pr.value.connectivity.name, pr.value.connectivity.port) %} {% endcall %} {% elseif ((pr.value.connectivity.type.stdlib == "OneToOne") || (pr.value.connectivity.type.stdlib == "http://nineml.net/9ML/1.0/connectionrules/OneToOne")) %} {% call one_to_one(pr.name, pr.value.source.populations, pr.value.destination.populations) %} {% endcall %} {% elseif ((pr.value.connectivity.type.stdlib == "AllToAll") || (pr.value.connectivity.type.stdlib == "http://nineml.net/9ML/1.0/connectionrules/AllToAll")) %} {% call all_to_all(pr.name, pr.value.source.populations, pr.value.destination.populations) %} {% endcall %} {% elseif ((pr.value.connectivity.type.stdlib == "FromFile") || (pr.value.connectivity.type.stdlib == "http://nineml.net/9ML/1.0/connectionrules/OneToOne")) %} {% call from_file(pr.name, pr.value.source.populations, pr.value.destination.populations, pr.value.rule.properties.filename.exprML) %} {% endcall %} {% elseif ((pr.value.connectivity.type.stdlib == "RandomFanIn") || (pr.value.connectivity.type.stdlib == "http://nineml.net/9ML/1.0/connectionrules/RandomFanIn")) && pr.value.connectivity.type.probability %} {% call random_fan_in(pr.name, pr.value.source.populations, pr.value.destination.populations, pr.value.connectivity.type.probability) %} {% endcall %} {% elseif ((pr.value.connectivity.type.stdlib == "RandomFanIn") || (pr.value.connectivity.type.stdlib == "http://nineml.net/9ML/1.0/connectionrules/RandomFanIn")) && pr.value.connectivity.type.number %} {% call random_fan_in_num(pr.name, pr.value.source.populations, pr.value.destination.populations, pr.value.connectivity.type.number) %} {% endcall %} {% elseif ((pr.value.connectivity.type.stdlib == "RandomFanOut") || (pr.value.connectivity.type.stdlib == "http://nineml.net/9ML/1.0/connectionrules/RandomFanOut")) && pr.value.connectivity.type.probability %} {% call random_fan_out(pr.name, pr.value.source.populations, pr.value.destination.populations, pr.value.connectivity.type.probability) %} {% endcall %} {% elseif ((pr.value.connectivity.type.stdlib == "RandomFanOut") || (pr.value.connectivity.type.stdlib == "http://nineml.net/9ML/1.0/connectionrules/RandomFanOut")) && pr.value.connectivity.type.number %} {% call random_fan_out_num(pr.name, pr.value.source.populations, pr.value.destination.populations, pr.value.connectivity.type.number) %} {% endcall %} {% endif %} {% endfor %} val S = ([ {% for pr in dict (group.projections) %} {% if pr.value.type == "event" %} Pr_{{pr.name}}{% if not loop.last %},{% endif %} {% endif %} {% endfor %} ]) {% for pr in dict (group.projections) %} {% if pr.value.type == "continuous" %} {% call range_map(pr.name, pr.value.source.populations, pr.value.destination.populations) %} {% endcall %} {% endif %} {% endfor %} val Elst = [ {% for pr in dict (group.projections) %} {% if pr.value.type == "continuous" %} ElecGraph.junctionMatrix ([N,N],ElecGraph.elecGraph ({{pr.name}}(srangemap_{{pr.name}},trangemap_{{pr.name}})) Pr_{{pr.name}}){% if not loop.last %},{% endif %} {% endif %} {% endfor %} ] val E = if List.null Elst then NONE else SOME Elst in S end) end