;; ;; An extension for translating NEMO models to NEST code. ;; ;; Copyright 2011-2013 Ivan Raikov and the Okinawa Institute of Science and Technology ;; ;; This program is free software: you can redistribute it and/or ;; modify it under the terms of the GNU General Public License as ;; published by the Free Software Foundation, either version 3 of the ;; License, or (at your option) any later version. ;; ;; This program is distributed in the hope that it will be useful, but ;; WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ;; General Public License for more details. ;; ;; A full copy of the GPL license can be found at ;; . ;; (module nemo-nest (nemo:nest-translator) (import scheme chicken utils data-structures ports extras srfi-1 srfi-13 srfi-69) (require-extension lolevel posix matchable strictly-pretty varsubst datatype nemo-core nemo-utils nemo-units nemo-geometry nemo-defaults nemo-constraints nemo-gate-complex nemo-synapse nemo-currents ) (define C++-ops `(+ - * / > < <= >= =)) (define (nest-name s) (let ((cs (string->list (->string s)))) (let loop ((lst (list)) (cs cs)) (if (null? cs) (string->symbol (list->string (reverse lst))) (let* ((c (car cs)) (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c) (else #\_)))) (loop (cons c1 lst) (cdr cs))))))) (define builtin-consts (append `(params) (map (lambda (x) (nest-name (s+ "M_" (first x)))) nemo:math-constants))) (define builtin-fns `(+ - * / pow neg abs atan asin acos sin cos exp ln sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube > < <= >= = and or round ceiling floor max min )) (define (rhsexpr/C++ expr) (match expr (('if . es) `(if . ,(map (lambda (x) (rhsexpr/C++ x)) es))) (('pow x y) (if (and (integer? y) (positive? y)) (if (> y 1) (let ((tmp (gensym "x"))) `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp))))) x) (if (and (number? y) (zero? y)) 1.0 expr))) ((s . es) (if (symbol? s) (cons (if (member s builtin-fns) s (nest-name s)) (map (lambda (x) (rhsexpr/C++ x)) es)) expr)) (id (if (symbol? id) (if (assoc id nemo:math-constants) (nest-name (s+ "M_" id)) (nest-name id)) id)) )) (define (nest-state-name n s) (nest-name (s+ n s))) (define-syntax pp (syntax-rules () ((pp indent val ...) (ppf indent (quasiquote val) ...)))) (define group/C++ (doc:block 2 (doc:text "(") (doc:text ")"))) (define block/C++ (doc:block 2 (doc:text "{") (doc:text "}"))) (define (stmt/C++ x) (match x (($ doc 'DocCons _ ($ doc 'DocText ";")) x) (else (doc:cons x (doc:text ";"))))) (define (ifthen/C++ c e1 e2) (doc:nest 2 (doc:connect (doc:group (doc:connect (doc:text "if") c)) (doc:connect (doc:nest 2 e1) (doc:nest 2 (doc:connect (doc:text "else") e2)))) )) (define (letblk/C++ e1 e2) (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2))) ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1))) (else (doc:connect (doc:group (doc:nest 2 (stmt/C++ e1))) (doc:group (doc:nest 2 e2)))))) (define (format-op/C++ indent op args) (let ((op1 (doc:text (->string op)))) (if (null? args) op1 (match args ((x) (doc:concat (list op1 x))) ((x y) (doc:concat (intersperse (list x op1 y) (doc:space)))) ((x y z) (doc:concat (intersperse (list x op1 y op1 z) (doc:space)))) (lst (let* ((n (length lst)) (n/2 (inexact->exact (round (/ n 2))))) (doc:concat (intersperse (list (format-op/C++ indent op (take lst n/2 )) op1 (format-op/C++ indent op (drop lst n/2 ))) (doc:space))))))))) (define (format-fncall/C++ indent op args) (let ((op1 (doc:text (->string op)))) (doc:cons op1 (group/C++ ((doc:list indent identity (lambda () (doc:text ", "))) args))))) (define (name-normalize expr) (match expr (('if c t e) `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e))) (('let bs e) `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e))) ((f . es) (cons (if (member f builtin-fns) f (nest-name f)) (map name-normalize es))) ((? symbol? ) (nest-name expr)) ((? atom? ) expr))) (define (canonicalize-expr/C++ expr) (let ((subst-convert (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term))) (let* ((expr1 (if-convert expr)) (expr2 (subst-convert expr1 subst-empty)) (expr3 (let-lift expr2)) (expr4 (name-normalize expr3))) expr4))) (define (format-expr/C++ indent expr . rest) (let-optionals rest ((rv #f)) (let ((indent+ (+ 2 indent))) (match expr (('let bindings body) (letblk/C++ (fold-right (lambda (x ax) (letblk/C++ (match (second x) (('if c t e) (ifthen/C++ (group/C++ (format-expr/C++ indent c)) (block/C++ (format-expr/C++ indent t (first x))) (block/C++ (format-expr/C++ indent e (first x))))) (else (stmt/C++ (format-op/C++ indent+ " = " (list (format-expr/C++ indent (first x) ) (format-expr/C++ indent (second x))))))) ax)) (doc:empty) bindings) (match body (('let _ _) (format-expr/C++ indent body rv)) (else (let ((body1 (doc:nest indent (format-expr/C++ indent body)))) (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) body1))) body1)))))) (('if . rest) (error 'format-expr/C++ "invalid if statement " expr)) ((op . rest) (let ((op (case op ((ln) 'log) ((abs) 'fabs) (else op)))) (let ((fe (if (member op C++-ops) (let ((mdiv? (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest)) (mul? (any (lambda (x) (match x (('* . _) #t) (else #f))) rest)) (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest))) (case op ((/) (format-op/C++ indent op (map (lambda (x) (let ((fx (format-expr/C++ indent+ x))) (if (or (symbol? x) (number? x)) fx (if (or mul? plmin?) (group/C++ fx) fx)))) rest))) ((*) (format-op/C++ indent op (map (lambda (x) (let ((fx (format-expr/C++ indent+ x))) (if (or (symbol? x) (number? x)) fx (if plmin? (group/C++ fx) fx)))) rest))) (else (format-op/C++ indent op (map (lambda (x) (let ((fx (format-expr/C++ indent+ x))) fx)) rest))))) (let ((op (case op ((neg) '-) (else op)))) (format-fncall/C++ indent op (map (lambda (x) (format-expr/C++ indent+ x)) rest)))))) (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe))) fe)))) (else (let ((fe (doc:text (->string expr)))) (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe))) fe))))))) (define (expr->string/C++ x . rest) (let-optionals rest ((rv #f) (width 72)) (sdoc->string (doc:format width (format-expr/C++ 2 x rv))))) (define (make-define-fn sysname ) (lambda (indent n proc) (let ((lst (procedure-data proc)) (indent+ (+ 2 indent))) (let ((rt (or (lookup-def 'rt lst) 'double)) (formals (lookup-def 'formals lst)) (vars (lookup-def 'vars lst)) (consts (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst))) (body (lookup-def 'body lst)) (rv (gensym 'rv))) (let ((argument-list (append (if (null? vars) '() (map (lambda (x) (sprintf "double ~A" (nest-name x))) vars)) '("const void* params")))) (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) "{" )) (let* ((body0 (rhsexpr/C++ body)) (body1 (canonicalize-expr/C++ (add-params-to-fncall body0 builtin-fns))) (lbs (enum-bnds body1 (list)))) (pp indent+ (double ,rv ";")) (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs) ";"))) (if (not (null? consts)) (begin (pp indent+ (double ,(slp ", " (map (compose nest-name car) consts)) ";") (,(sprintf "const ~A::Parameters_" sysname) "& p = " ,(sprintf " *(reinterpret_cast< const ~A::Parameters_ *>(params));" sysname))) (for-each (lambda (x) (let ((n (car x))) (pp indent+ (,(nest-name n) = ,(sprintf "p.~A;" (nest-name n)))))) consts) )) (pp indent+ ,(expr->string/C++ body1 (nest-name rv))) (pp indent+ ,(s+ "return " rv ";")) )) (pp indent "}")) ))) (define (make-define-fn-header sysname ) (lambda (indent n proc) (let ((lst (procedure-data proc)) (indent+ (+ 2 indent))) (let ((rt (or (lookup-def 'rt lst) 'double)) (formals (lookup-def 'formals lst)) (vars (lookup-def 'vars lst)) (consts (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst))) (body (lookup-def 'body lst)) (rv (gensym 'rv))) (let ((argument-list (append (if (null? vars) '() (map (lambda (x) (sprintf "double ~A" (nest-name x))) vars)) '("const void* params")))) (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) ";" )) )) ))) (define (ith v i) (sprintf "Ith(~A,~A)" v i)) (define (output-prelude sysname steady-state-index-map indent) (pp indent ,#<#EOF ##include "#{sysname}.h" ##include "exceptions.h" ##include "network.h" ##include "dict.h" ##include "integerdatum.h" ##include "doubledatum.h" ##include "dictutils.h" ##include "numerics.h" ##include ##include "universal_data_logger_impl.h" ##include ##include ##include ##include ##include EOF )) (define cvode-prelude-header #<= 0 * opt == 2 means function allocates memory so check if returned * NULL pointer */ static int check_flag(void *flagvalue, const char *funcname, int opt) { int *errflag; /* Check if CVode function returned NULL pointer - no memory allocated */ if (opt == 0 && flagvalue == NULL) { fprintf(stderr, "\nCVode error: %s() failed - returned NULL pointer\n\n", funcname); return(1); } /* Check if flag < 0 */ else if (opt == 1) { errflag = (int *) flagvalue; if (*errflag < 0) { fprintf(stderr, "\nCVode error: %s() failed with flag = %d\n\n", funcname, *errflag); return(1); } } /* Check if function returned NULL pointer - no memory allocated */ else if (opt == 2 && flagvalue == NULL) { fprintf(stderr, "\nMemory error: %s() failed - returned NULL pointer\n\n", funcname); return(1); } return(0); } void adjust_zero_crossings (N_Vector v, double abstol) { int i; for (i = 0; i < NV_LENGTH_S(v); i++) { if (fabs(NV_Ith_S(v,i)) < abstol) NV_Ith_S(v,i) = 0; } return; } EOF ) (define (output-event-handle sysname state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info indent indent+) (let ((isyns (lookup-def 'i-synapses synapse-info)) (pscs (lookup-def 'post-synaptic-conductances synapse-info))) (pp indent ,nl (,#<#EOF void #{sysname}::handle(SpikeEvent & e) { int flag; assert(e.get_delay() > 0); flag = 0; EOF )) (for-each (lambda (isyn psc) (let ((wscale (fourth isyn)) (wthreshold (fifth isyn))) (pp indent+ (,#<#EOF if ((flag==0) && (e.get_rport() == #{(nest-name (first psc))}_SPIKE_RECEPTOR ) && (e.get_weight() > #(or (and wthreshold (sprintf "P_.~A" (nest-name wthreshold))) 0.0))) { B_.spike_#(nest-name (second psc)).add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), fabs(e.get_weight()) * e.get_multiplicity()); flag = 1; } EOF )) )) isyns pscs) (pp indent ,nl #\} ,nl) (pp indent (,#<#EOF void #{sysname}::handle(CurrentEvent& e) { assert(e.get_delay() > 0); const double_t c = e.get_current(); const double_t w = e.get_weight(); B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), w * c); } void #{sysname}::handle(DataLoggingRequest& e) { B_.logger_.handle(e); } EOF )) )) (define (output-header sysname method state-index-map steady-state-index-map defaults external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info indent indent+) (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info))) (pp indent ("#include \"nest.h\"") ("#include \"event.h\"") ("#include \"archiving_node.h\"") ("#include \"ring_buffer.h\"") ("#include \"connection.h\"") ("#include \"universal_data_logger.h\"") ("#include \"recordables_map.h\"") ) (case method ((leapfrog) (pp indent ("#define Ith(v,i) (v[i]) /* Ith component in a vector */") )) ((cvode) (pp indent ("#include /* prototypes for CVODE fcts., consts. */") ("#include /* serial N_Vector types, fcts., macros */") ("#include /* prototype for CVDiag */") ("#include /* definition of type realtype */") ("#define Ith(v,i) NV_Ith_S(v,i) /* Ith component in a vector */") )) ((ida) (pp indent ("#include /* prototypes for IDA fcts., consts. */") ("#include ") ("#include /* serial N_Vector types, fcts., macros */") ("#define Ith(v,i) NV_Ith_S(v,i) /* Ith component in a vector */") )) ((gsl) (pp indent ("#include ") ("#include ") ("#include ") ("#include ") ("#define Ith(v,i) (v[i])") )) ) (pp indent ,nl) (pp indent "namespace nest {" ,nl) (case method ((cvode) (pp indent ,cvode-prelude-header)) ((ida) (pp indent ,ida-prelude-header))) (case method ((leapfrog) (pp indent ("typedef int (*rhsfn_t)(double t, const double *y, double *ydot, void *user_data);") ,nl (extern "\"C\"" int ,(sprintf "~A_dynamics" sysname) "(double, const double *, double *, void*)" #\;) ,nl)) ((cvode) (pp indent (extern "\"C\"" int ,(sprintf "~A_dynamics" sysname) "(double, const N_Vector, N_Vector, void*)" #\;) ,nl (extern "\"C\"" int ,(sprintf "~A_event" sysname) "(double, N_Vector, double *, void*)" #\;) ,nl )) ((ida) (pp indent (extern "\"C\"" int ,(sprintf "~A_residual" sysname) "(double, N_Vector, N_Vector, N_Vector, void*)" #\;) ,nl (extern "\"C\"" int ,(sprintf "~A_event" sysname) "(double, N_Vector, N_Vector, double *, void*)" #\;) ,nl )) (else (pp indent (extern "\"C\"" int ,(sprintf "~A_dynamics" sysname) "(double, const double*, double*, void*)" #\;) ,nl)) ) (pp indent (,(sprintf "class ~A:" sysname) "public Archiving_Node { ")) (pp indent+ ( ,#<#EOF public: ~#{sysname} (); #{sysname} (const #{sysname} &); #{sysname} (); EOF )) (pp indent ( #<; friend class UniversalDataLogger<#{sysname}>; EOF )) (pp indent+ ,nl "struct Parameters_ { ") (for-each (lambda (x) (pp indent+ (,(sprintf "double ~A;" x)))) (append (map (compose nest-name first) const-defs) (map (compose nest-name first) defaults))) (pp indent+ ("Parameters_();") ("void get(DictionaryDatum&) const;") ("void set(const DictionaryDatum&);")) (pp indent+ "};") (pp indent+ ,nl "struct State_ { ") (pp indent+ ,nl "enum StateVecElems {" (,(slp ", " (map (lambda (x) (let ((n (string-upcase (->string (first x)))) (ni (second x))) (sprintf "~A = ~A" n ni))) state-index-map))) "};") (pp indent+ ,nl (,(sprintf "double y_[~A];" (length state-index-map))) "State_(const Parameters_& p);" "State_(const State_& s);" "State_& operator=(const State_& s);" "void get(DictionaryDatum&) const;" "void set(const DictionaryDatum&, const Parameters_&);") (case method ((gsl leapfrog) (pp indent+ ,nl "int_t r_; /* refractory counter */"))) (pp indent+ "};") (case method ((cvode ida) (pp indent+ ,nl "struct Variables_ {};")) (else (pp indent+ ,nl "struct Variables_ { int_t RefractoryCounts_; double U_old_; /* for spike-detection */ };")) ) (pp indent+ ,nl "struct Buffers_ { ") (pp indent+ ,nl (,(sprintf "Buffers_(~A&);" sysname)) (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname)) (,(sprintf "UniversalDataLogger<~A> logger_;" sysname))) (pp indent+ ,nl (,(case method ((leapfrog) #<" "double_t get_y_elem_() const { return S_.y_[elem]; }" "Parameters_ P_;" "State_ S_;" "Variables_ V_;" "Buffers_ B_;" (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname)) "}; ") (pp indent+ (,#<#EOF inline port #{sysname}::check_connection(Connection& c, port receptor_type) { SpikeEvent e; e.set_sender(*this); c.check_event(e); return c.get_target()->connect_sender(e, receptor_type); } inline port #{sysname}::connect_sender(SpikeEvent&, port receptor_type) { if ( receptor_type < MIN_SPIKE_RECEPTOR || receptor_type >= SUP_SPIKE_RECEPTOR ) { if ( receptor_type < 0 || receptor_type >= SUP_SPIKE_RECEPTOR ) throw UnknownReceptorType(receptor_type, get_name()); else throw IncompatibleReceptorType(receptor_type, get_name(), "SpikeEvent"); } return receptor_type; } inline port #{sysname}::connect_sender(CurrentEvent&, port receptor_type) { if (receptor_type != 0) throw UnknownReceptorType(receptor_type, get_name()); return 0; } inline port #{sysname}::connect_sender(DataLoggingRequest& dlr, port receptor_type) { if (receptor_type != 0) throw UnknownReceptorType(receptor_type, get_name()); return B_.logger_.connect_logging_device(dlr, recordablesMap_); } inline void #{sysname}::get_status(DictionaryDatum &d) const { P_.get(d); S_.get(d); Archiving_Node::get_status(d); (*d)[names::recordables] = recordablesMap_.get_list(); def(d, names::t_spike, get_spiketime_ms()); DictionaryDatum receptor_dict_ = new Dictionary(); #{(string-concatenate (map (lambda (psc) (sprintf " (*receptor_dict_)[Name(\"~A\")] = ~A_SPIKE_RECEPTOR;~%" (first psc) (first psc))) pscs))} (*d)[names::receptor_types] = receptor_dict_; } inline void #{sysname}::set_status(const DictionaryDatum &d) { Parameters_ ptmp = P_; // temporary copy in case of errors ptmp.set(d); // throws if BadProperty State_ stmp = S_; // temporary copy in case of errors stmp.set(d, ptmp); // throws if BadProperty // We now know that (ptmp, stmp) are consistent. We do not // write them back to (P_, S_) before we are also sure that // the properties to be set in the parent class are internally // consistent. Archiving_Node::set_status(d); // if we get here, temporaries contain consistent set of properties P_ = ptmp; S_ = stmp; calibrate(); } EOF )) (pp indent "}" ,nl) ) ) (define (output-jac sysname method state-index-map pool-ions i-eqs v-eq indent indent+) ;; Diagonal Jacobian approximation: (f(s+.01) - f(s))/.001 ;; CVODE: ;; typedef int (*CVDlsDenseJacFn)(long int N, realtype t, ;; N_Vector y, N_Vector fy, ;; DlsMat Jac, void *user_data, ;; N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); ;; GSL: ;; int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params); (case method ((gsl) (pp indent ,nl (,(sprintf "extern \"C\" int ~A_jacobian" sysname) (,(slp ", " `("double t" "const double y[]" "double *dfdy" "double dfdt[]" "void* pnode" )))) #\{ ("// cast the node ptr to an object of the proper type") ("assert(pnode);") ("const " ,sysname "& node = *(reinterpret_cast<" ,sysname "*>(pnode));") (,sysname "& vnode = *(reinterpret_cast<" ,sysname "*>(pnode));") ("// state is a reference to the model state ") (struct ,(sprintf "~A::Buffers_" sysname) "*b;") ("b = &(vnode.B_);") (,#<#EOF for (int i = 0; i < b->N; i++) { b->u[i] = y[i] + 0.01; } #{sysname}_dynamics(t, b->u, b->jac, pnode); for (int i = 0; i < b->N; i++) { dfdt[i*b->N + i] = (b->jac[i] - dfdy[i]) / .001; } return 0; EOF ) #\} )) )) (define (output-cvode-event sysname method imports const-defs state-index-map external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs defaults i-eqs v-eq indent indent+) (let ((vi (lookup-def 'v state-index-map))) (pp indent ,nl (,(sprintf "extern \"C\" int ~A_event" sysname) (,(slp ", " `("double t" "N_Vector y" "double *g" "void* pnode" ))) #\{ )) (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;") ,nl ("// S is shorthand for the type that describes the model state ") (typedef ,(sprintf "~A::State_" sysname) "S;") ,nl ("// cast the node ptr to an object of the proper type") ("assert(pnode);") ("const " ,sysname "& node = *(reinterpret_cast<" ,sysname "*>(pnode));") ,nl ("// params is a reference to the model parameters ") (const struct ,(sprintf "~A::Parameters_" sysname) "*params;") ("params = &(node.P_);") ,nl ("// state is a reference to the model state ") (const struct ,(sprintf "~A::State_" sysname) "*state;") ("state = &(node.S_);") ) (let ((vt (lookup-def 'V_t defaults))) (if (and vi vt) (pp indent+ (,(expr->string/C++ (ith 'y vi) 'v)) (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt)) (,(expr->string/C++ '(- v vt) "g[0]")) ,nl ("return 0; ")) (pp indent+ (,(expr->string/C++ -1.0 "g[0]"))) )) (pp indent (#\})) )) (define (output-ida-event sysname method imports const-defs state-index-map external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs defaults i-eqs v-eq indent indent+) (let ((vi (lookup-def 'v state-index-map))) (pp indent ,nl (,(sprintf "extern \"C\" int ~A_event" sysname) (,(slp ", " `("double t" "N_Vector y" "N_Vector yp" "double *g" "void* pnode" ))) #\{ )) (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;") ,nl ("// S is shorthand for the type that describes the model state ") (typedef ,(sprintf "~A::State_" sysname) "S;") ,nl ("// cast the node ptr to an object of the proper type") ("assert(pnode);") ("const " ,sysname "& node = *(reinterpret_cast<" ,sysname "*>(pnode));") ,nl ("// params is a reference to the model parameters ") (const struct ,(sprintf "~A::Parameters_" sysname) "*params;") ("params = &(node.P_);") ,nl ("// state is a reference to the model state ") (const struct ,(sprintf "~A::State_" sysname) "*state;") ("state = &(node.S_);") ) (let ((vt (lookup-def 'V_t defaults))) (if (and vi vt) (pp indent+ (,(expr->string/C++ (ith 'y vi) 'v)) (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt)) (,(expr->string/C++ '(- v vt) "g[0]")) ,nl ("return 0; ")) (pp indent+ (,(expr->string/C++ -1.0 "g[0]"))) )) (pp indent (#\})) )) (define (cvode-solve spike-handle) #<#EOF int N = NV_LENGTH_S (B_.y); tout = (current_steps+lag)*h; // adaptive step integration while (tt < tout) { const int status = CVode (B_.sys_, tout, B_.y, &tt, CV_NORMAL); switch (status) { case CV_SUCCESS: continue; case CV_ROOT_RETURN: #{spike-handle}; default: throw CVodeSolverFailure (get_name(), 0); } } for (int i = 0; i < N; i++) { S_.y_[i] = Ith(B_.y,i); } EOF ) (define (ida-solve spike-handle) #<#EOF int N = NV_LENGTH_S (B_.y); tout = (current_steps+lag)*h; // adaptive step integration while (tt < tout) { const int status = IDASolve (B_.sys_, tout, &tt, B_.y, B_.yp, IDA_NORMAL); switch (status) { case IDA_SUCCESS: continue; case IDA_ROOT_RETURN: #{spike-handle}; case IDA_TSTOP_RETURN: break; default: throw IDASolverFailure (get_name(), 0); } } for (int i = 0; i < N; i++) { S_.y_[i] = Ith(B_.y,i); } EOF ) (define (gsl-solve vi spike-handle) #<#EOF V_.U_old_ = S_.y_[#{vi}]; while (tt < h) { const int status = gsl_odeiv2_evolve_apply (B_.e_, B_.c_, B_.s_, &B_.sys_, // system of ODE &tt, // from t... h, // ...to t=t+h &B_.IntegrationStep_ , // integration window (written on!) S_.y_); // neuron state if ( status != GSL_SUCCESS ) throw GSLSolverFailure(get_name(), status); } #{spike-handle} EOF ) (define (leapfrog-solve vi spike-handle) #<#EOF V_.U_old_ = S_.y_[#{vi}]; tt = (current_steps+lag)*h; tout = (current_steps+lag+1)*h; while (tt < tout) { B_.sys_ (tt, S_.y_, B_.dvdt, reinterpret_cast(this)); for (int i = 0; i < B_.N; i++) { B_.u[i] = S_.y_[i] + 0.5 * h * B_.dvdt[i]; } B_.sys_ (tt+0.5*h, B_.u, B_.dvdt, reinterpret_cast(this)); for (int i = 0; i < B_.N; i++) { B_.v[i] = S_.y_[i] + h * B_.dvdt[i]; } tt += h; } for (int i = 0; i < B_.N; i++) { S_.y_[i] = B_.v[i]; } #{spike-handle} EOF ) (define (output-synapse-transients sysname method state-index-map const-defs transient-event-defs event-external-eq-defs synapse-info imports constraints indent indent+) (let ( (c-units (map (lambda (x) (let ((n (first x)) (v (second x))) (list (nest-name n) v))) (lookup-def 'c-units constraints))) (isyns (lookup-def 'i-synapses synapse-info)) (pscs (lookup-def 'post-synaptic-conductances synapse-info)) (psc-transients (map (lambda (lst) (map nest-name lst)) (lookup-def 'psc-transients synapse-info))) (getstate (case method ((cvode ida) (lambda (i) (sprintf "Ith(B_.y,~A)" i))) (else (lambda (i) (sprintf "Ith(S_.y_,~A)" i))))) ) (if (not (= (length event-external-eq-defs) (length pscs))) (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" event-external-eq-defs pscs)) (for-each (lambda (psc transients) (let* ( (ltransient-event-defs (filter (lambda (x) (member (first x) transients)) transient-event-defs)) (levent-external-eq-def (car (fold (lambda (def ax) (let* ((n (nest-name (first def)) ) (b (second def)) (events (let ((fvs (enum-freevars b '() '()))) (filter (lambda (x) (member (first x) fvs)) event-external-eq-defs))) ) (append events ax))) '() ltransient-event-defs))) (lconsts (delete-duplicates (fold (lambda (def ax) (let* ((n (nest-name (first def)) ) (b (second def)) (consts (let ((fvs (enum-freevars b '() '()))) (filter (lambda (x) (member (first x) fvs)) const-defs))) ) (append consts ax))) '() ltransient-event-defs) (lambda (x y) (equal? (first x) (first y)))) )) (pp indent (,(sprintf "inline int ~A::~A_transients (long_t lag) {" sysname (first psc)))) (if (not (null? lconsts)) (pp indent+ (double ,(slp ", " (map (compose nest-name car) lconsts)) ";"))) (let ((n (second levent-external-eq-def))) (pp indent+ (,(sprintf "double_t ~A;" n)) (,(sprintf "~A = B_.spike_~A.get_value(lag);" n (nest-name (second psc)))) ) ) (pp indent+ (,(sprintf "if (~A > 0.0)" (second levent-external-eq-def))) (#\{)) (let* ((n (nest-name (first levent-external-eq-def))) (nu (lookup-def n c-units)) (nscale (and nu (nemo:unit-scale nu))) (b (second levent-external-eq-def)) ) (pp (+ 2 indent+) (,(sprintf "double_t ~A;" n)) (,(expr->string/C++ (if nscale `(* ,nscale ,b) b) n))) ) (for-each (lambda (def) (let* ((n (nest-name (first def)) ) (ni (lookup-def n state-index-map)) (b (second def)) (consts (let ((fvs (enum-freevars b '() '()))) (filter (lambda (x) (member (first x) fvs)) const-defs))) ) (if (not (null? consts)) (for-each (lambda (x) (let ((n (first x))) (pp (+ 2 indent+) (,(nest-name n) = ,(sprintf "P_.~A;" (nest-name n)))))) consts) ) (pp (+ 2 indent+ ) (,(sprintf "double_t ~A;" n)) ,(if ni (expr->string/C++ (getstate ni) n) '()) ,(expr->string/C++ b n) ,(if ni (expr->string/C++ n (getstate ni)) '()) ) )) ltransient-event-defs) (pp (+ 2 indent+) "return 1;"); (pp indent+ (#\})) (pp indent+ "return 0;"); (pp indent (#\})) )) pscs psc-transients) )) (define (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports indent) (let ( (isyns (lookup-def 'i-synapses synapse-info)) (pscs (lookup-def 'post-synaptic-conductances synapse-info)) ) (pp indent (,(sprintf "int events = 0;"))) (if (not (= (length event-external-eq-defs) (length pscs))) (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" event-external-eq-defs pscs)) (for-each (lambda (psc) (pp indent (,(sprintf "events = events + ~A_transients (lag);" (first psc)))) ) pscs) (case method ((cvode) (pp indent ( "/* Reinitializes CVode state if any synaptic events have occurred */") ("if (events > 0)") (#\{) (" int status = CVodeReInit (B_.sys_, tt, B_.y);") (" if (check_flag(&status, \"CVodeReInit\", 1)) throw CVodeSolverFailure (get_name(), status);") (#\}) )) ((ida) (pp indent ( "/* Reinitializes IDA state if any synaptic events have occurred */") ("if (events > 0)") (#\{) (" int status = IDAReInit (B_.sys_, tt, B_.y, B_.yp);") (" if (check_flag(&status, \"IDAReInit\", 1)) throw IDASolverFailure (get_name(), status);") (#\}) )) ((gsl) (pp indent ( "/* Resets the GSL stepping function if any synaptic events have occurred */") ("if (events > 0)") (#\{) (" gsl_odeiv2_step_reset (B_.s_);") (#\}) )) ) )) (define (spike-handle method defaults vi) (case method ((cvode ida) #<send(*this, se, lag); adjust_zero_crossings(B_.y, 1e-7); continue; } EOF ) (else (if (lookup-def 'V_t defaults) (sprintf #< 0 ) S_.r_--; else if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A]) { S_.r_ = V_.RefractoryCounts_; set_spiketime(Time::step(origin.get_steps()+lag+1)); SpikeEvent se; network()->send(*this, se, lag); } EOF vi vi) )) )) (define (output-update sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs transient-event-defs i-eqs pool-ions perm-ions synapse-info event-external-eq-defs defaults imports indent indent+) (let ((vi (lookup-def 'v state-index-map))) (pp indent (,(sprintf "void ~A::update(Time const & origin, const long_t from, const long_t to)" sysname) #\{)) (pp indent+ #<= 0 && (delay) from < Scheduler::get_min_delay()); assert(from < to); double tout; long_t current_steps = origin.get_steps(); for ( long_t lag = from ; lag < to ; ++lag ) { double h = B_.step_; double tt = 0.0 ; EOF ) (case method ((leapfrog) (begin (pp indent+ (,(leapfrog-solve vi (spike-handle method defaults vi)))) )) ((cvode) (begin (pp indent+ (,(cvode-solve (spike-handle method defaults vi) )) ,nl) (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) )) ((ida) (begin (pp indent+ (,(ida-solve (spike-handle method defaults vi) )) ,nl) (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) )) ((gsl) (begin (pp indent+ (,(gsl-solve vi (spike-handle method defaults vi))) ,nl) (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) )) ) (pp (+ 6 indent+) ("B_.I_stim_ = B_.currents_.get_value(lag);") ("B_.logger_.record_data(current_steps + lag);")) (pp indent+ #\}) (pp indent #\}) )) (define (output-nodes-leapfrog sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info defaults indent indent+) (let ((N (length state-index-map)) (pscs (lookup-def 'post-synaptic-conductances synapse-info)) ) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname}() : Archiving_Node(), P_(), S_(P_), B_(*this) { recordablesMap_.create(); } EOF )) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname} (const #{sysname}& n) : Archiving_Node(n), P_(n.P_), S_(n.S_), B_(n.B_, *this) { } EOF )) (pp indent #<#EOF #{sysname}::~#{sysname} () { if ( B_.u != NULL ) free (B_.u); if ( B_.v != NULL ) free (B_.v); if ( B_.dvdt != NULL ) free (B_.dvdt); } EOF ) (pp indent ,nl (,#<#EOF void #{sysname}::init_node_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); P_ = pr.P_; S_ = State_(P_); } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_state_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); S_ = State_(pr.P_); } EOF )) )) (define (output-nodes-cvode sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info defaults indent indent+) (let ( (N (length state-index-map)) (pscs (lookup-def 'post-synaptic-conductances synapse-info)) ) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname}() : Archiving_Node(), P_(), S_(P_), B_(*this) { recordablesMap_.create(); } EOF )) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname}(const #{sysname}& n) : Archiving_Node(n), P_(n.P_), S_(n.S_), B_(n.B_, *this) { } EOF )) (pp indent (,#<#EOF #{sysname}::~#{sysname}() { if ( B_.sys_ ) { /* Free y vector */ N_VDestroy_Serial(B_.y); /* Free integrator memory */ // CVodeFree(&B_.sys_); } } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_node_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); P_ = pr.P_; S_ = State_(P_); } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_state_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); S_ = State_(pr.P_); } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_buffers_() { EOF )) (for-each (lambda (psc) (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc))))) pscs) (pp indent+ (#<(this)); if (check_flag(&status, "CVodeSetUserData", 1)) throw CVodeSolverFailure (get_name(), status); /* Initializes diagonal linear solver. */ status = CVDiag (B_.sys_); if (check_flag(&status, "CVDiag", 1)) throw CVodeSolverFailure (get_name(), status); } EOF )) (pp indent (,#<#EOF void #{sysname}::calibrate() { B_.logger_.init(); } EOF )) )) (define (output-nodes-ida sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info defaults indent indent+) (let ( (N (length state-index-map)) (pscs (lookup-def 'post-synaptic-conductances synapse-info)) ) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname}() : Archiving_Node(), P_(), S_(P_), B_(*this) { recordablesMap_.create(); } EOF )) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname}(const #{sysname}& n) : Archiving_Node(n), P_(n.P_), S_(n.S_), B_(n.B_, *this) { } EOF )) (pp indent (,#<#EOF #{sysname}::~#{sysname}() { if ( B_.sys_ ) { /* Free y vector */ N_VDestroy_Serial(B_.y); N_VDestroy_Serial(B_.yp); /* Free integrator memory */ // IDAFree(&B_.sys_); } } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_node_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); P_ = pr.P_; S_ = State_(P_); } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_state_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); S_ = State_(pr.P_); } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_buffers_() { EOF )) (for-each (lambda (psc) (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc))))) pscs) (pp indent+ (#<(this)); if (check_flag(&status, "IDASetUserData", 1)) throw IDASolverFailure (get_name(), status); /* Initializes dense linear solver. */ status = IDADense (B_.sys_); if (check_flag(&status, "IDADense", 1)) throw IDASolverFailure (get_name(), status); } EOF )) (pp indent (,#<#EOF void #{sysname}::calibrate() { B_.logger_.init(); } EOF )) )) (define (output-nodes-gsl sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info defaults indent indent+) (let ((N (length state-index-map)) (iv (lookup-def 'v state-index-map)) (pscs (lookup-def 'post-synaptic-conductances synapse-info))) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname}() : Archiving_Node(), P_(), S_(P_), B_(*this) { recordablesMap_.create(); } EOF )) (pp indent ,nl (,#<#EOF #{sysname}::#{sysname}(const #{sysname}& n) : Archiving_Node(n), P_(n.P_), S_(n.S_), B_(n.B_, *this) { } EOF )) (pp indent ,#<#EOF #{sysname}::~#{sysname} () { // GSL structs only allocated by init_nodes_(), so we need to protect destruction if ( B_.s_ != NULL) gsl_odeiv2_step_free (B_.s_); if ( B_.c_ != NULL) gsl_odeiv2_control_free (B_.c_); if ( B_.e_ != NULL) gsl_odeiv2_evolve_free (B_.e_); if ( B_.u != NULL) free (B_.u); if ( B_.jac != NULL) free (B_.jac); } EOF ) (pp indent ,nl (,#<#EOF void #{sysname}::init_node_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); P_ = pr.P_; S_ = State_(P_); } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_state_(const Node& proto) { const #{sysname}& pr = downcast<#{sysname}>(proto); S_ = State_(pr.P_); } EOF )) (pp indent ,nl (,#<#EOF void #{sysname}::init_buffers_() { EOF )) (for-each (lambda (psc) (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc))))) pscs) (pp indent+ (#<(this); B_.u = (double *)malloc(sizeof(double) * B_.N); assert (B_.u); B_.jac = (double *)malloc(sizeof(double) * B_.N); assert (B_.jac); EOF )) (pp indent ,nl #\}) (pp indent (,#<#EOF void #{sysname}::calibrate() { B_.logger_.init(); #{(if iv (sprintf "V_.U_old_ = S_.y_[~A];" iv) "")} #{(if (lookup-def 't_ref defaults) "V_.RefractoryCounts_ = Time(Time::ms(P_.t_ref)).get_steps();" "V_.RefractoryCounts_ = 0;")} } EOF )) )) (define (output-buffers sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info defaults indent indent+) (let ((N (length state-index-map))) (case method ((leapfrog) (pp indent ,nl (,#<#EOF #{sysname}::Buffers_::Buffers_(#{sysname}& n) : logger_(n), sys_(0), N(0), u(0), v(0), dvdt(0) { // Initialization of the remaining members is deferred to // init_buffers_(). } EOF ))) ((cvode ida) (pp indent ,nl (,#<#EOF #{sysname}::Buffers_::Buffers_(#{sysname}& n) : logger_(n), sys_(0) { // Initialization of the remaining members is deferred to // init_buffers_(). } EOF ) (,#<#EOF #{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n) : logger_(n), sys_(0) { // Initialization of the remaining members is deferred to // init_buffers_(). } EOF ))) ((gsl) (pp indent ,nl ( ,#<#EOF #{sysname}::Buffers_::Buffers_(#{sysname}& n) : logger_(n), s_(0), c_(0), e_(0), N(0), u(0), jac(0) { // Initialization of the remaining members is deferred to // init_buffers_(). } EOF ) ( ,#<#EOF #{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n) : logger_(n), s_(0), c_(0), e_(0), N(0), u(0), jac(0) { // Initialization of the remaining members is deferred to // init_buffers_(). } EOF ))) )) ) (define (output-accessors+modifiers sysname imports defaults state-index-map const-defs asgn-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions constraints indent indent+) (let ((c-eqs (lookup-def 'c-eqs constraints)) (c-units (map (lambda (x) (let ((n (first x)) (v (second x))) (list (nest-name n) v))) (lookup-def 'c-units constraints))) ) (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) )) (pp indent #\{) (for-each (lambda (def) (let ((n (first def))) (pp indent+ (,(sprintf "def(d, ~S, ~A);" (->string n) n))))) (append const-defs defaults)) (pp indent #\}) (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) )) (pp indent #\{) (for-each (lambda (def) (let* ((n (first def)) (nu (lookup-def n c-units)) (scale (and nu (nemo:unit-scale nu))) ) (pp indent+ (,(sprintf "updateValue(d, ~S, ~A);" (->string n) n))) (if scale (pp indent+ (,(sprintf "~A = ~A * ~A;" n scale n )))) )) (append const-defs defaults)) (pp indent #\}) (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) )) (pp indent #\{) (for-each (lambda (def) (let* ( (n (first def)) (i (second def)) ) (pp indent+ (,(sprintf "def(d, ~S, y_[~A]);" (->string n) i))) )) state-index-map) (let ((vi (lookup-def 'v state-index-map))) (if vi (pp indent+ (,(sprintf "def(d, names::V_m, y_[~A]);" vi) )) )) (pp indent #\}) (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) )) (pp indent #\{) (for-each (lambda (def) (let* ((n (first def)) (i (second def)) (nu (lookup-def n c-units)) (scale (and nu (nemo:unit-scale nu))) ) (pp indent+ (,(sprintf "updateValue(d, ~S, y_[~A]);" (->string n) i))) (if scale (pp indent+ (,(sprintf "y_[~A] = ~A * y_[~A];" i scale i )))) )) state-index-map) (let ((vi (lookup-def 'v state-index-map))) (if vi (pp indent+ (,(sprintf "updateValue(d, names::V_m, y_[~A]);" vi) )) )) (pp indent #\}) )) (define (output-init sysname method state-index-map steady-state-index-map imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+) (let* ((N (length state-index-map)) (c-eqs (lookup-def 'c-eqs constraints)) (c-units (map (lambda (x) (let ((n (first x)) (v (second x))) (list (nest-name n) v))) (lookup-def 'c-units constraints))) (i-eqs (map (lambda (def) (list (first def) (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns))) i-eqs)) (init-eqs (append (map (lambda (def) (let ((n (first def)) (b (second def))) (list (nest-name n) (nest-name b)))) external-eq-defs) asgn-eq-defs init-eq-defs (map (lambda (pool-ion) (let ((n (pool-ion-in pool-ion)) (b (pool-ion-inq pool-ion))) (list n b))) pool-ions) (map (lambda (pool-ion) (let ((n (pool-ion-out pool-ion)) (b (pool-ion-outq pool-ion))) (list n b))) pool-ions) )) (init-dag (map (lambda (def) (cons (first def) (enum-freevars (second def) '() '()))) init-eqs)) (init-order (reverse (topological-sort init-dag (lambda (x y) (string=? (->string x) (->string y)))))) (init-locals (find-locals (map second (append init-eqs i-eqs reaction-eq-defs)))) ) (pp indent ,nl (,(sprintf "~A::State_::State_" sysname) (const Parameters_& p)) ( #\{)) (pp indent+ (double ,(slp ", " (delete-duplicates (map ->string (filter (lambda (x) (not (member x builtin-consts))) (append init-locals init-order (map first external-eq-defs) (map pool-ion-in pool-ions) (map pool-ion-out pool-ions) (map first i-eqs) (map first state-index-map) (map first const-defs) (map first reaction-eq-defs) ))) string=?)) ";") ("const Parameters_ *params = &p;") ) (case method ((gsl leapfrog) (pp indent+ "r_ = 0;" ,nl))) (pp indent+ ,(sprintf "memset(y_,0,~A*sizeof(double));" N)) (for-each (lambda (x) (let* ((n (first x))) (pp indent+ ,(expr->string/C++ (sprintf "p.~A" n) n)))) const-defs) (for-each (lambda (eq) (match-let (((op left right) eq)) (pp indent+ ,(sprintf "if (!((~A) ~A (~A))) { throw BadProperty (~S); }; " (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ left))) op (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ right))) (sprintf "Constraint ~A is not satisfied." eq) )) )) c-eqs) (let ((vi (lookup-def 'v state-index-map)) (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0))) (if (and vi vrest) (pp indent+ ,(expr->string/C++ vrest 'v )))) (for-each (lambda (n) (let ((b (lookup-def n init-eqs))) (if b (pp indent+ ,(expr->string/C++ b (nest-name n)))))) init-order) (for-each (lambda (def) (let* ((n (first def)) (ni (lookup-def n state-index-map))) (if ni (pp indent+ ,(expr->string/C++ n (sprintf "y_[~A]" ni)))))) init-eq-defs) #| (if (not (null? steady-state-index-map)) (begin (gsl-fminimizer sysname N indent+) (for-each (lambda (def) (let* ((n (first def)) (si (lookup-def n steady-state-index-map)) (ni (lookup-def n state-index-map))) (if si (begin (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n)) (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")") (s+ "y_[" ni "]"))))))) rate-eq-defs) (pp indent+ "gsl_vector_free (ys);") )) |# (for-each (lambda (def) (let ((n (first def)) (b (second def))) (if (not (lookup-def n init-eq-defs)) (pp indent+ ,(expr->string/C++ b n))))) reaction-eq-defs) (for-each (lambda (def) (pp indent+ ,(expr->string/C++ (second def) (first def)))) i-eqs) (let ((vi (lookup-def 'v state-index-map))) (if vi (pp indent+ ,(expr->string/C++ 'v (sprintf "y_[~A]" vi))))) (pp indent #\}) (pp indent ,nl (,(sprintf "~A::State_::State_" sysname) (const State_& s) ) ( #\{) ,(case method ((gsl leapfrog) `("r_ = s.r_;")) (else "")) (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N)) ( #\}) ,nl ) (pp indent ,nl (,(sprintf "~A::State_& ~A::State_::operator=(const State_& s)" sysname sysname)) ( #\{) ) (case method ((gsl leapfrog) (pp indent+ "r_ = s.r_;" ,nl))) (pp indent+ ( ,#<#EOF assert(this != &s); for ( size_t i = 0 ; i < #{N} ; ++i ) y_[i] = s.y_[i]; EOF )) (pp indent+ "return *this;") (pp indent #\}) )) (define (output-parameters sysname imports const-defs external-eq-defs defaults indent indent+) (let* ((parameter-defs (map (lambda (def) (list (first def) (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns))) const-defs)) (parameter-locals (find-locals (map second parameter-defs))) (const-exprs (map (lambda (def) (let* ((n (first def)) (b (second def))) (s+ (nest-name n) " (" (expr->string/C++ b) ")"))) const-defs) ) (default-exprs (map (lambda (def) (let ((n (first def)) (b (second def))) (expr->string/C++ (nest-name b) (nest-name n)))) defaults)) ) (pp indent ,nl (,(sprintf "~A::Parameters_::Parameters_" sysname) () ":")) (if (not (null? parameter-locals)) (pp indent+ (double ,(slp ", " parameter-locals) ";"))) (pp indent+ ,(slp ",\n" const-exprs)) (pp indent #\{ ,(slp "\n" default-exprs) #\}) )) (define (output-recordables sysname state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions indent indent+) (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname))) (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname))) (pp indent #\{) (for-each (lambda (def) (let ((n (first def)) (i (second def))) (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);" (->string n) sysname sysname (string-upcase (->string n))))) )) state-index-map) (if (lookup-def 'v state-index-map) (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);" sysname sysname (string-upcase (->string 'v)))))) (pp indent #\}) ) (define (output-dy sysname method imports const-defs state-index-map external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs pool-ions i-eqs v-eq constraints indent indent+) (let* ( (c-units (map (lambda (x) (let ((n (first x)) (v (second x))) (list (nest-name n) v))) (lookup-def 'c-units constraints))) (i-eqs (map (lambda (def) (list (first def) (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns))) i-eqs)) (v-eq (and v-eq (list (first v-eq) (add-params-to-fncall (canonicalize-expr/C++ (second v-eq)) builtin-fns)))) (eqs (append external-eq-defs asgn-eq-defs reaction-eq-defs (map (lambda (pool-ion) (let ((n (pool-ion-in pool-ion)) (b (pool-ion-inq pool-ion))) (list n b))) pool-ions) (map (lambda (pool-ion) (let ((n (pool-ion-out pool-ion)) (b (pool-ion-outq pool-ion))) (list n b))) pool-ions) i-eqs )) (eq-dag (map (lambda (def) (cons (first def) (enum-freevars (second def) '() '()))) eqs)) (eq-order (reverse (topological-sort eq-dag (lambda (x y) (string=? (->string x) (->string y)))))) (eq-locals (find-locals (map second (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs)) (append i-eqs rate-eq-defs eqs))))) ) (case method ((leapfrog) (pp indent ,nl ( ,(sprintf "extern \"C\" int ~A_dynamics" sysname) (,(slp ", " `("double t" "const double y[]" "double f[]" "void* pnode" ))) #\{ ))) ((cvode) (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname) (,(slp ", " `("double t" "N_Vector y" "N_Vector f" "void* pnode" ))) #\{ ))) ((ida) (pp indent ,nl (,(sprintf "extern \"C\" int ~A_residual" sysname) (,(slp ", " `("double t" "N_Vector y" "N_Vector yp" "N_Vector f" "void* pnode" ))) #\{ ))) ((gsl) (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname) (,(slp ", " `("double t" "const double y[]" "double f[]" "void* pnode" ))) #\{ ))) ) (pp indent+ (double ,(slp ", " (delete-duplicates (map (compose ->string nest-name) (filter (lambda (x) (not (member x builtin-consts))) (append eq-locals eq-order (map first i-eqs) (map first external-eq-defs) (map first state-index-map) (map first const-defs) ))) string=?)) ";")) (pp indent+ ,nl ("// S is shorthand for the type that describes the model state ") (typedef ,(sprintf "~A::State_" sysname) "S;") ,nl ("// cast the node ptr to an object of the proper type") ("assert(pnode);") ("const " ,sysname "& node = *(reinterpret_cast<" ,sysname "*>(pnode));") (,sysname "& vnode = *(reinterpret_cast<" ,sysname "*>(pnode));") ,nl ("// params is a reference to the model parameters ") (const struct ,(sprintf "~A::Parameters_" sysname) "*params;") ("params = &(node.P_);") ,nl ("// state is a reference to the model state ") (struct ,(sprintf "~A::State_" sysname) "*state;") ("state = &(vnode.S_);") ,nl ("// y[] must be the state vector supplied by the integrator, ") ("// not the state vector in the node, node.S_.y[]. ") ,nl ) (for-each (lambda (def) (let ((n (first def)) ) (pp indent+ ,(expr->string/C++ (sprintf "params->~A" n) n)))) const-defs) (let ((vi (lookup-def 'v state-index-map))) (if vi (pp indent+ ,(expr->string/C++ (ith 'y vi) 'v)))) (for-each (lambda (def) (let* ((n (first def)) (ni (lookup-def n state-index-map))) (pp indent+ ,(expr->string/C++ (ith 'y ni) (nest-name n))))) rate-eq-defs) (for-each (lambda (n) (let ((b (lookup-def n eqs))) (if b (pp indent+ ,(expr->string/C++ b (nest-name n)))))) eq-order) (for-each (lambda (def) (let* ((n (first def)) (b (second def)) (fv (ith 'f (lookup-def n state-index-map))) ) (pp indent+ ,(s+ "// rate equation " n) ,(expr->string/C++ b fv)))) rate-eq-defs) (if v-eq (let ((vi (lookup-def 'v state-index-map))) (if vi (pp indent+ ,(expr->string/C++ (second v-eq) (ith 'f vi))) ))) (case method ((leapfrog) (pp indent+ ,nl ("return 0;"))) ((cvode) (pp indent+ ,nl ("return 0;"))) ((gsl) (pp indent+ ,nl ("return GSL_SUCCESS;"))) ) (pp indent #\}) )) (define (nemo:nest-translator sys . rest) (define (cid x) (second x)) (define (cn x) (first x)) (let-optionals rest ((dirname ".") (method #f)) (let ((method (or method 'gsl))) (if (not (member method '(gsl cvode leapfrog))) (nemo:error 'nemo:nest-translator ": unknown method " method)) (match-let ((($ nemo:quantity 'DISPATCH dis) (hash-table-ref sys (nemo-intern 'dispatch)))) (let ((imports ((dis 'imports) sys)) (exports ((dis 'exports) sys))) (let* ((indent 0) (indent+ (+ 2 indent )) (sysname (nest-name ((dis 'sysname) sys))) (prefix (->string sysname)) (deps* ((dis 'depgraph*) sys)) (consts ((dis 'consts) sys)) (asgns ((dis 'asgns) sys)) (states ((dis 'states) sys)) (reactions ((dis 'reactions) sys)) (defuns ((dis 'defuns) sys)) (components ((dis 'components) sys)) (g (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g)) (poset (vector->list ((dis 'depgraph->bfs-dist-poset) g))) (const-defs (filter-map (lambda (nv) (and (not (member (first nv) builtin-consts)) (let ((v1 (canonicalize-expr/C++ (second nv)))) (list (nest-name (first nv)) v1)))) consts)) (defaults (nemo:defaults-query sys)) (geometry (nemo:geometry-query sys)) (gate-complex-info (nemo:gate-complex-query sys)) (perm-ions (map (match-lambda ((comp i e erev val) `(,comp ,(nest-name i) ,(nest-name e) ,erev))) (lookup-def 'perm-ions gate-complex-info))) (acc-ions (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out))))) (lookup-def 'acc-ions gate-complex-info))) (epools (lookup-def 'pool-ions gate-complex-info)) (pool-ions (pool-ion-name-map nest-name epools)) (comprc (any (match-lambda ((name 'membrane-tau id) (list name id)) (else #f)) components)) (compcap (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components)) (mrc (or (and comprc (car ((dis 'component-exports) sys (cid comprc)))) (lookup-def 'membrane-tau defaults) (lookup-def 'tau_m defaults) (and compcap (car ((dis 'component-exports) sys (cid compcap)))) (lookup-def 'membrane-capacitance defaults) (lookup-def 'C_m defaults) )) (soma-geometry (lookup-def 'soma geometry)) (marea (and soma-geometry (third soma-geometry))) (gate-complexes (lookup-def 'gate-complexes gate-complex-info)) (synapse-info (nemo:post-synaptic-conductance-query sys)) (pscs (lookup-def 'post-synaptic-conductances synapse-info)) (i-syns (lookup-def 'i-synapses synapse-info)) (i-gates (lookup-def 'i-gates gate-complex-info)) (i-defs (nemo:ionic-current-definitions gate-complexes i-gates i-syns pscs marea (lambda (x) (state-power sys x)) (lambda (x) ((dis 'component-exports) sys x)) (lambda (x) ((dis 'component-subcomps) sys x)) nest-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)) (i-eqs (lookup-def 'i-eqs i-defs)) (i-names (lookup-def 'i-names i-defs)) (constraints (nemo:constraint-definitions gate-complexes i-gates i-syns pscs marea imports (lambda (x) (state-power sys x)) (lambda (x) (quantity-unit sys x)) (lambda (x) ((dis 'component-exports) sys x)) (lambda (x) ((dis 'component-subcomps) sys x)) nest-name)) (external-eq-defs (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++ namespace-filter: (lambda (x) (not (equal? x 'event))))) (event-external-eq-defs (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++ namespace-filter: (lambda (x) (equal? x 'event)))) (asgn-eq-defs (poset->asgn-eq-defs* poset sys nest-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)) (rate-eq-defs (reverse (poset->rate-eq-defs* poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))) (reaction-eq-defs (poset->reaction-eq-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++)) (transient-event-defs (poset->transient-event-defs poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)) (init-eq-defs (poset->init-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)) (conserve-eq-defs (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) (poset->state-conserve-eq-defs poset sys nest-name nest-state-name))) (imports-sans-v (filter (lambda (x) (not (equal? 'v (first x)))) imports)) (v-eq (and (not (null? i-names)) (let ((istim "(node.B_.I_stim_)" )) (cond ((and mrc marea) (list 'v (rhsexpr/C++ `(/ (+ (* ,istim (/ 100. ,marea)) (* -1e3 ,(sum i-names))) ,mrc)))) (marea (list 'v (rhsexpr/C++ `(+ (* ,istim (/ 100. ,marea)) (* -1e-3 ,(sum i-names)))))) (mrc (list 'v (rhsexpr/C++ `(/ (+ ,istim (* -1e-3 ,(sum i-names))) ,mrc)))) (else (list 'v (rhsexpr/C++ `(+ ,istim (* -1e-3 ,(sum i-names)))))) )) )) (state-index-map (let ((acc (fold (lambda (def ax) (let ((st-name (first def))) (list (+ 1 (first ax)) (cons `(,st-name ,(first ax)) (second ax))))) (list 0 (list)) (if (lookup-def 'v imports) (cons (list 'v) rate-eq-defs) rate-eq-defs) ))) (second acc))) (steady-state-index-map (let ((acc (fold (lambda (def ax) (let ((st-name (first def))) (if (not (alist-ref st-name init-eq-defs)) (list (+ 1 (first ax)) (cons `(,st-name ,(first ax)) (second ax))) ax))) (list 0 (list)) rate-eq-defs))) (second acc))) (dfenv (map (lambda (x) (let ((n (first x))) (list n (nest-name (sprintf "d_~A" n ))))) defuns)) ) (for-each (lambda (a) (let ((acc-ion (car a))) (if (assoc acc-ion perm-ions) (nemo:error 'nemo:nest-translator ": ion species " acc-ion " cannot be declared as both accumulating and permeating")))) acc-ions) (let ((cpp-output (open-output-file (make-output-fname dirname prefix ".cpp"))) (hpp-output (open-output-file (make-output-fname dirname prefix ".h")))) ;; include files and other prelude definitions (with-output-to-port cpp-output (lambda () (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl)) (output-prelude sysname steady-state-index-map indent) )) (with-output-to-port cpp-output (lambda () (pp indent ,nl "namespace nest {" ,nl))) (case method ((cvode) (with-output-to-port cpp-output (lambda () (pp indent ,sundials-prelude))))) ;; user-defined functions (let ((define-fn (make-define-fn sysname)) (define-fn-header (make-define-fn-header sysname))) (for-each (lambda (fndef) (if (not (member (car fndef) builtin-fns)) (with-output-to-port cpp-output (lambda () (apply define-fn-header (cons indent fndef)) (pp indent ,nl))) )) defuns) (for-each (lambda (fndef) (if (not (member (car fndef) builtin-fns)) (with-output-to-port cpp-output (lambda () (apply define-fn (cons indent fndef)) (pp indent ,nl))) )) defuns) ) ;; derivative function (with-output-to-port cpp-output (lambda () (output-dy sysname method imports-sans-v const-defs state-index-map external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs pool-ions i-eqs v-eq constraints indent indent+) )) ;; Jacobian function (with-output-to-port cpp-output (lambda () (output-jac sysname method state-index-map pool-ions i-eqs v-eq indent indent+) )) ;; system recordables (with-output-to-port cpp-output (lambda () (output-recordables sysname state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions indent indent+) (pp indent ,nl) )) ;; system parameters (with-output-to-port cpp-output (lambda () (output-parameters sysname imports-sans-v const-defs external-eq-defs defaults indent indent+) )) ;; initial values function (with-output-to-port cpp-output (lambda () (output-init sysname method state-index-map steady-state-index-map imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+) (pp indent ,nl) )) ;; accessors and modifiers for states and parameters (with-output-to-port cpp-output (lambda () (output-accessors+modifiers sysname imports-sans-v defaults state-index-map const-defs asgn-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions constraints indent indent+) (pp indent ,nl) )) ;; buffer and node initialization (with-output-to-port cpp-output (lambda () (output-buffers sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info defaults indent indent+) (pp indent ,nl) )) (with-output-to-port cpp-output (lambda () ((case method ((leapfrog) output-nodes-leapfrog) ((gsl) output-nodes-gsl) ((cvode) output-nodes-cvode) ((ida) output-nodes-ida) ) sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info defaults indent indent+) (pp indent ,nl) )) ;; spike threshold detect (cvode method only) (case method ((cvode) (with-output-to-port cpp-output (lambda () (output-cvode-event sysname method imports const-defs state-index-map external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs defaults i-eqs v-eq indent indent+) (pp indent ,nl) ))) ) ;; synaptic event transients (with-output-to-port cpp-output (lambda () (output-synapse-transients sysname method state-index-map const-defs transient-event-defs event-external-eq-defs synapse-info imports constraints indent indent+) (pp indent ,nl) )) ;; state update (with-output-to-port cpp-output (lambda () (output-update sysname method state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs transient-event-defs i-eqs pool-ions perm-ions synapse-info event-external-eq-defs defaults imports indent indent+) (pp indent ,nl) )) ;; spike handler (with-output-to-port cpp-output (lambda () (output-event-handle sysname state-index-map steady-state-index-map const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info indent indent+) (pp indent ,nl) )) (with-output-to-port cpp-output (lambda () (pp indent "}" ,nl))) (with-output-to-port hpp-output (lambda () (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl)) (output-header sysname method state-index-map steady-state-index-map defaults external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs reaction-eq-defs i-eqs pool-ions perm-ions synapse-info indent indent+) (pp indent ,nl) )) (close-output-port cpp-output) (close-output-port hpp-output) )) )) )) ))