Changeset 27726 in project


Ignore:
Timestamp:
10/29/12 10:00:36 (8 years ago)
Author:
Ivan Raikov
Message:

nemo: current status of nest backend

File:
1 edited

Legend:

Unmodified
Added
Removed
  • release/4/nemo/trunk/nemo-nest.scm

    r27687 r27726  
    482482
    483483
    484 (define (output-dy sysname const-defs state-index-map
    485                    external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
    486                    pool-ions mcap i-eqs v-eq indent indent+)
    487 
    488 
    489   (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
    490                       (,(slp ", " `("double t"
    491                                     "const double y[]"
    492                                     "double f[]"
    493                                     "void* pnode"
    494                                     )))
    495                       #\{
    496                       ))
    497 
    498 
    499   (let* (
    500          (asgn-eq-defs
    501           (map
    502            (lambda (def) (list (first def)
    503                                (add-params-to-fncall
    504                                 (canonicalize-expr/C++ (second def)))))
    505            asgn-eq-defs))
    506 
    507          (rate-eq-defs
    508           (map
    509            (lambda (def) (list (first def)
    510                                (add-params-to-fncall
    511                                 (canonicalize-expr/C++ (second def)))))
    512            rate-eq-defs))
    513 
    514          (reaction-eq-defs
    515           (map
    516            (lambda (def) (list (first def)
    517                                (add-params-to-fncall
    518                                 (canonicalize-expr/C++ (second def)))))
    519            reaction-eq-defs))
    520 
    521          (i-eqs
    522           (map
    523            (lambda (def) (list (first def)
    524                                (add-params-to-fncall (canonicalize-expr/C++ (second def)))))
    525            i-eqs))
    526 
    527          (v-eq
    528           (and v-eq
    529                (list (first v-eq)
    530                      (add-params-to-fncall (canonicalize-expr/C++ (second v-eq))))))
    531 
    532          (eqs
    533           (append
    534            
    535            asgn-eq-defs
    536            reaction-eq-defs
    537                  
    538            (map (lambda (pool-ion)
    539                   (let ((n (third pool-ion))
    540                         (b (first pool-ion)))
    541                     (list n b)))
    542                 pool-ions)))
    543          
    544          (eq-dag
    545           (map (lambda (def)
    546                  (cons (first def) (enum-freevars (second def) '() '())))
    547                eqs))
    548 
    549          (eq-order
    550           (reverse
    551            (topological-sort eq-dag
    552                              (lambda (x y) (string=? (->string x) (->string y))))))
    553 
    554          (eq-locals  (find-locals
    555                       (map second
    556                            (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
    557                                (append i-eqs rate-eq-defs eqs)))))
    558          )
    559 
    560 
    561 
    562     (pp indent+ (double ,(slp ", " (delete-duplicates
    563                                     (map ->string
    564                                          (filter (lambda (x) (not (member x nest-builtin-consts)))
    565                                                  (append
    566                                                   eq-locals
    567                                                   eq-order
    568                                                   (map first external-eq-defs)
    569                                                   (map first state-index-map)
    570                                                   (map first i-eqs)
    571                                                   (map first const-defs))))
    572                                     string=?))
    573                         ";"))
    574 
    575     (pp indent+ ,nl
    576         ("// S is shorthand for the type that describes the model state ")
    577         (typedef ,(s+ sysname "::State_") "S;")
    578        
    579         ,nl
    580        
    581         ("// cast the node ptr to an object of the proper type")
    582         ("assert(pnode);")
    583         ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
    584        
    585         ,nl
    586 
    587         ("// params is a reference to the model parameters ")
    588         (const struct ,(s+ sysname "::Parameters_") "*params;")
    589         ("params = &(node.P_);")
    590        
    591         ,nl
    592 
    593         ("// y[] must be the state vector supplied by the integrator, ")
    594         ("// not the state vector in the node, node.S_.y[]. ")
    595         ,nl
    596         )
    597 
    598     (for-each (lambda (def)
    599                 (let ((n (first def))
    600                       (b (second def)))
    601                   (pp indent+ ,(expr->string/C++ b (nest-name n)))))
    602               external-eq-defs)
    603    
    604     (for-each (lambda (def)
    605                 (let ((n (first def)) )
    606                   (pp indent+
    607                       ,(expr->string/C++
    608                         (s+ "params->" n) n))))
    609               const-defs)
    610 
    611     (let ((vi (lookup-def 'v state-index-map)))
    612       (if vi (pp indent+ ,(expr->string/C++ (s+ "y[" vi "]") 'v))))
    613 
    614     (for-each (lambda (def)
    615                 (let ((n (first def)) )
    616                   (pp indent+
    617                       ,(expr->string/C++
    618                         (s+ "y[" (lookup-def n state-index-map) "]") n))))
    619               rate-eq-defs)
    620 
    621     (for-each (lambda (n)
    622                 (let ((b  (lookup-def n eqs)))
    623                   (if b (pp indent+
    624                             ,(expr->string/C++ b (nest-name n))))))
    625               eq-order)
    626 
    627     (for-each (lambda (def)
    628                 (let ((n (s+ "f[" (lookup-def (first def) state-index-map) "]"))
    629                       (b (second def)))
    630                   (pp indent+ ,(s+ "// rate equation " (first def))
    631                       ,(expr->string/C++ b n))))
    632               rate-eq-defs)
    633    
    634     (for-each (lambda (def)
    635                 (pp indent+ ,(expr->string/C++ (second def) (first def))))
    636             i-eqs)
    637    
    638     (if v-eq
    639         (let ((vi (lookup-def 'v state-index-map)))
    640           (if vi
    641               (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "f[" vi "]")))
    642               )))
    643    
    644     (pp indent+ ,nl ("return GSL_SUCCESS;"))
    645     (pp indent  #\})
    646 
    647     ))
    648 
    649 
    650 
    651 (define (output-parameters sysname const-defs indent indent+)
    652 
    653   (let* ((parameter-defs
    654           (map
    655            (lambda (def) (list (first def)
    656                                (add-params-to-fncall (canonicalize-expr/C++ (second def)))))
    657            const-defs))
    658          (parameter-locals  (find-locals (map second parameter-defs))))
    659 
    660     (pp indent ,nl (,(s+ sysname "::Parameters_::Parameters_") () ":"))
    661    
    662     (let ((const-exprs
    663            (slp ",\n"
    664                 (map (lambda (def)
    665                        (let* ((n  (first def)) (b (second def)))
    666                          (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
    667                      const-defs) )))
    668      
    669       (if (not (null? parameter-locals))
    670           (pp indent+ (double ,(slp ", " parameter-locals) ";")))
    671 
    672       (pp indent+ ,const-exprs)
    673      
    674       (pp indent ("{}"))
    675      
    676       )))
    677 
    678 
    679 (define (output-init sysname state-index-map steady-state-index-map
    680                      external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    681                      reaction-eq-defs i-eqs v-eq pool-ions perm-ions indent indent+)
    682 
    683  
    684   (let* ((N (length state-index-map))
    685 
    686          (asgn-eq-defs
    687           (map
    688            (lambda (def) (list (first def)
    689                                (add-params-to-fncall
    690                                 (canonicalize-expr/C++ (second def)))))
    691            asgn-eq-defs))
    692          
    693          (init-eq-defs
    694           (map
    695            (lambda (def) (list (first def)
    696                                (add-params-to-fncall
    697                                 (canonicalize-expr/C++ (second def)))))
    698            init-eq-defs))
    699 
    700          (reaction-eq-defs
    701           (map
    702            (lambda (def) (list (first def)
    703                                (add-params-to-fncall
    704                                 (canonicalize-expr/C++ (second def)))))
    705            reaction-eq-defs))
    706 
    707          (i-eqs
    708           (map
    709            (lambda (def) (list (first def)
    710                                (add-params-to-fncall
    711                                 (canonicalize-expr/C++ (second def)))))
    712            i-eqs))
    713 
    714          (v-eq
    715           (and v-eq
    716                (list (first v-eq)
    717                      (add-params-to-fncall
    718                       (canonicalize-expr/C++ (second v-eq))))))
    719          
    720          (init-eqs
    721           (append
    722            
    723            asgn-eq-defs
    724            init-eq-defs
    725 
    726            (map (lambda (pool-ion)
    727                   (let ((n (third pool-ion))
    728                         (b (first pool-ion)))
    729                     (list n b)))
    730                 pool-ions)))
    731 
    732          (init-dag
    733           (map (lambda (def)
    734                  (cons (first def) (enum-freevars (second def) '() '())))
    735                init-eqs))
    736 
    737          (init-order
    738           (reverse
    739            (topological-sort init-dag
    740                              (lambda (x y) (string=? (->string x) (->string y))))))
    741 
    742          (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
    743          )
    744    
    745     (pp indent ,nl (,(s+ sysname "::State_::State_") (const Parameters_& p) ": r_(0)"))
    746    
    747     (pp indent  #\{)
    748 
    749     (pp indent+
    750         (double ,(slp ", " (delete-duplicates
    751                             (map ->string
    752                                  (filter (lambda (x) (not (member x nest-builtin-consts)))
    753                                          (append
    754                                           init-locals
    755                                           init-order
    756                                           (map first external-eq-defs)
    757                                           (map first state-index-map)
    758                                           (map first i-eqs)
    759                                           (map first const-defs)
    760                                           (map first reaction-eq-defs)
    761                                           )))
    762                                  string=?))
    763                       ";")
    764         ("const Parameters_ *params = &p;")
    765         )
    766 
    767     (for-each (lambda (def)
    768                 (let ((n (first def)) )
    769                   (pp indent+
    770                       ,(expr->string/C++
    771                         (s+ "p." n) n))))
    772               const-defs)
    773 
    774     (let ((vi (lookup-def 'v state-index-map))
    775           (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
    776       (if (and vi vrest)
    777           (pp indent+ ,(expr->string/C++  vrest 'v ))))
    778 
    779     (for-each (lambda (def)
    780                 (let ((n (first def))
    781                       (b (second def)))
    782                   (pp indent+ ,(expr->string/C++ b (nest-name n)))))
    783               external-eq-defs)
    784 
    785     (for-each (lambda (n)
    786                 (let ((b  (lookup-def n init-eqs)))
    787                   (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
    788               init-order)
    789    
    790     (for-each (lambda (def)
    791                 (let* ((n  (first def))
    792                        (ni (lookup-def n state-index-map)))
    793                   (if ni (pp indent+ ,(expr->string/C++ n  (s+ "y_[" ni "]"))))))
    794               init-eq-defs)
    795 #|
    796     (if (not (null? steady-state-index-map))
    797         (begin
    798           (gsl-fminimizer sysname N indent+)
    799           (for-each
    800            (lambda (def)
    801              (let* ((n (first def))
    802                     (si (lookup-def n steady-state-index-map))
    803                     (ni (lookup-def n state-index-map)))
    804                (if si (begin
    805                         (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n))
    806                         (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")")
    807                                                        (s+ "y_[" ni "]")))))))
    808            rate-eq-defs)
    809       (pp indent+ "gsl_vector_free (ys);")
    810       ))
    811 |#
    812 
    813     (for-each
    814      (lambda (def)
    815        (let ((n (first def)) (b (second def)))
    816          (if (not (lookup-def n init-eq-defs))
    817              (pp indent+ ,(expr->string/C++ b n)))))
    818      reaction-eq-defs)
    819    
    820     (for-each
    821      (lambda (def)
    822        (pp indent+ ,(expr->string/C++ (second def) (first def))))
    823      i-eqs)
    824 
    825     (if v-eq
    826         (let ((vi (lookup-def 'v state-index-map)))
    827           (if vi
    828               (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "y_[" vi "]"))))))
    829      
    830     (pp indent  #\})
    831    
    832     (pp indent ,nl (,(s+ sysname "::State_::State_") (const State_& s) ": r_(s.r_)"))
    833     (pp indent  #\{)
    834     (pp indent+ (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N)))
    835     (pp indent  #\})
    836 
    837     (pp indent ,nl (,(s+ sysname "::State_& " sysname "::State_::operator=(const State_& s)")))
    838 
    839     (pp indent  #\{)
    840 
    841 
    842     (pp indent+ (,(sprintf #<<EOF
    843    assert(this != &s); 
    844      for ( size_t i = 0 ; i < ~A ; ++i )
    845        y_[i] = s.y_[i];
    846      r_ = s.r_;
    847      return *this;
    848 EOF
    849 N)))
    850 
    851      (pp indent  #\})
    852 
    853 ))
    854 
    855 
    856 
    857 (define (gsl-fminimizer sysname N indent)
    858 
    859   (pp indent (,(sprintf #<<EOF
    860 
    861   int minimizer_status, minimizer_iterc;
    862 
    863   const gsl_multimin_fdfminimizer_type *T;
    864   gsl_multimin_fdfminimizer *S;
    865 
    866   gsl_vector *ys;
    867   gsl_multimin_function_fdf F;
    868 
    869   F.f   = &~A_dynamics;
    870   F.df  = &~A_dynamics;
    871   F.fdf = &~A_dynamics;
    872   F.n   = ~A;
    873   F.params = (void *)&p;
    874 
    875   ys = gsl_vector_alloc (~A);
    876 
    877   for (int i = 0; i < ~A; i++)
    878   {
    879      gsl_vector_set (ys, i, 0.0);
    880   }
    881 
    882   T = gsl_multimin_fdfminimizer_conjugate_fr;
    883   S = gsl_multimin_fdfminimizer_alloc (T, ~A);
    884 
    885   gsl_multimin_fdfminimizer_set (S, &F, ys, 0.01, 1e-4);
    886 
    887   do
    888     {
    889       minimizer_iterc++;
    890       minimizer_status = gsl_multimin_fdfminimizer_iterate (S);
    891 
    892       if (minimizer_status)
    893         break;
    894 
    895       minimizer_status = gsl_multimin_test_gradient (S->gradient, 1e-3);
    896 
    897     }
    898   while (minimizer_status == GSL_CONTINUE && minimizer_iterc < 100);
    899 
    900   gsl_multimin_fdfminimizer_free (S);
    901 
    902 EOF
    903 sysname sysname sysname N N N N)))
    904 )
    905 
    906 
    907 
    908 (define (output-accessors+modifiers
    909          sysname state-index-map const-defs asgn-eq-defs rate-eq-defs
    910          reaction-eq-defs i-eqs v-eq pool-ions perm-ions indent indent+)
    911 
    912   (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
    913   (pp indent  #\{)
    914 
    915   (for-each
    916    (lambda (def)
    917      (let ((n (first def)))
    918        (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
    919    const-defs)
    920 
    921   (pp indent  #\})
    922 
    923   (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
    924   (pp indent  #\{)
    925 
    926   (for-each
    927    (lambda (def)
    928      (let ((n (first def)))
    929        (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))))
    930    const-defs)
    931 
    932   (pp indent  #\})
    933 
    934   (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
    935   (pp indent  #\{)
    936 
    937   (for-each
    938    (lambda (def)
    939      (let ((n (first def)) (i (second def)))
    940        (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
    941    state-index-map)
    942 
    943   (pp indent  #\})
    944 
    945   (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
    946   (pp indent  #\{)
    947 
    948   (for-each
    949    (lambda (def)
    950      (let ((n (first def)) (i (second def)))
    951        (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
    952    state-index-map)
    953 
    954 
    955   (pp indent  #\})
    956 
    957 )
    958 
    959 (define (output-buffers+nodes
    960          sysname state-index-map steady-state-index-map
    961          const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    962          reaction-eq-defs i-eqs pool-ions perm-ions
    963          synapse-info
    964          indent indent+)
    965 
    966   (let ((N (length state-index-map)))
    967 
    968   (pp indent ,nl  (,(sprintf #<<EOF
    969 ~A::Buffers_::Buffers_(~A& n)
    970     : logger_(n),
    971       s_(0),
    972       c_(0),
    973       e_(0)
    974 {
    975     // Initialization of the remaining members is deferred to
    976     // init_buffers_().
    977 }
    978 EOF
    979   sysname sysname)))
    980 
    981   (pp indent  ,nl  (,(sprintf #<<EOF
    982 ~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
    983     : logger_(n),
    984       s_(0),
    985       c_(0),
    986       e_(0)
    987 {
    988     // Initialization of the remaining members is deferred to
    989     // init_buffers_().
    990 }
    991 EOF
    992   sysname sysname)))
    993 
    994 
    995   (pp indent  ,nl  (,(sprintf #<<EOF
    996 ~A::~A()
    997     : Archiving_Node(),
    998       P_(),
    999       S_(P_),
    1000       B_(*this)
    1001 {
    1002     recordablesMap_.create();
    1003 }
    1004 EOF
    1005   sysname sysname)))
    1006 
    1007 
    1008   (pp indent  ,nl  (,(sprintf #<<EOF
    1009 ~A::~A(const ~A& n)
    1010     : Archiving_Node(n),
    1011       P_(n.P_),
    1012       S_(n.S_),
    1013       B_(n.B_, *this)
    1014 {
    1015 }
    1016 EOF
    1017   sysname sysname sysname)))
    1018 
    1019   (pp indent (,(s+ sysname "::~" sysname "()")))
    1020   (pp indent #<<EOF
    1021 {
    1022     // GSL structs only allocated by init_nodes_(), so we need to protect destruction
    1023     if ( B_.s_ ) gsl_odeiv_step_free(B_.s_);
    1024     if ( B_.c_ ) gsl_odeiv_control_free(B_.c_);
    1025     if ( B_.e_ ) gsl_odeiv_evolve_free(B_.e_);
    1026 }
    1027 EOF
    1028   )
    1029 
    1030   (pp indent ,nl (,(sprintf #<<EOF
    1031   void ~A::init_node_(const Node& proto)
    1032 {
    1033     const ~A& pr = downcast<~A>(proto);
    1034     P_ = pr.P_;
    1035     S_ = pr.S_;
    1036 }
    1037 EOF
    1038   sysname sysname sysname)))
    1039 
    1040   (pp indent  ,nl (,(sprintf #<<EOF
    1041 void ~A::init_state_(const Node& proto)
    1042 {
    1043     const ~A& pr = downcast<~A>(proto);
    1044     S_ = pr.S_;
    1045 }
    1046 EOF
    1047   sysname sysname sysname)))
    1048 
    1049   (pp indent ,nl (,(sprintf #<<EOF
    1050 void ~A::init_buffers_()
    1051 {
    1052 EOF
    1053  sysname)))
    1054 
    1055      (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
    1056       (for-each
    1057        (lambda (psc)
    1058          (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
    1059        pscs))
    1060 
    1061     (pp indent+ (,(sprintf #<<EOF
    1062     B_.currents_.clear();           
    1063     Archiving_Node::clear_history();
    1064 
    1065     B_.logger_.reset();
    1066 
    1067     B_.step_ = Time::get_resolution().get_ms();
    1068     B_.IntegrationStep_ = B_.step_;
    1069 
    1070     B_.I_stim_ = 0.0;
    1071 
    1072     static const gsl_odeiv_step_type* T1 = gsl_odeiv_step_rkf45;
    1073  
    1074     if ( B_.s_ == 0 )
    1075       B_.s_ = gsl_odeiv_step_alloc (T1, ~A);
    1076     else
    1077       gsl_odeiv_step_reset(B_.s_);
    1078    
    1079     if ( B_.c_ == 0 ) 
    1080       B_.c_ = gsl_odeiv_control_y_new (1e-3, 0.0);
    1081     else
    1082       gsl_odeiv_control_init(B_.c_, 1e-3, 0.0, 1.0, 0.0);
    1083    
    1084     if ( B_.e_ == 0 ) 
    1085       B_.e_ = gsl_odeiv_evolve_alloc(~A);
    1086     else
    1087       gsl_odeiv_evolve_reset(B_.e_);
    1088  
    1089     B_.sys_.function  = ~A_dynamics;
    1090     B_.sys_.jacobian  = 0;
    1091     B_.sys_.dimension = ~A;
    1092     B_.sys_.params    = reinterpret_cast<void*>(this);
    1093 }
    1094 EOF
    1095   N N sysname N)))
    1096 
    1097 
    1098   (let ((iv (lookup-def 'v state-index-map)))
    1099     (if iv
    1100         (pp indent ,nl (,(sprintf #<<EOF
    1101 void ~A::calibrate()
    1102 {
    1103     B_.logger_.init(); 
    1104     V_.RefractoryCounts_ = 20;
    1105     V_.U_old_ = S_.y_[~A];
    1106 }
    1107 EOF
    1108   sysname iv)))
    1109   ))
    1110 
    1111 ))
    1112 
    1113 
    1114 (define (output-update
    1115          sysname state-index-map steady-state-index-map
    1116          const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1117          reaction-eq-defs i-eqs pool-ions perm-ions
    1118          synapse-info
    1119          indent indent+)
    1120  
    1121   (let* ((vi (lookup-def 'v state-index-map))
    1122          (threshold (sprintf #<<EOF
    1123             // (threshold && maximum)
    1124             if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
    1125               {
    1126                 S_.r_ = V_.RefractoryCounts_;
    1127                
    1128                 set_spiketime(Time::step(origin.get_steps()+lag+1));
    1129                
    1130                 SpikeEvent se;
    1131                 network()->send(*this, se, lag);
    1132               }
    1133 EOF
    1134         vi vi)))
    1135    
    1136     (pp indent  (,(sprintf #<<EOF
    1137 void ~A::update(Time const & origin, const long_t from, const long_t to)
    1138 EOF
    1139   sysname)))
    1140 
    1141     (pp indent  #\{)
    1142 
    1143     (pp indent+ (,(sprintf #<<EOF
    1144 assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
    1145     assert(from < to);
    1146 
    1147     for ( long_t lag = from ; lag < to ; ++lag )
    1148       {
    1149    
    1150         double tt = 0.0 ; //it's all relative!
    1151         V_.U_old_ = S_.y_[~A];
    1152 
    1153    
    1154         // adaptive step integration
    1155         while (tt < B_.step_)
    1156         {
    1157           const int status = gsl_odeiv_evolve_apply(B_.e_, B_.c_, B_.s_,
    1158                                  &B_.sys_,              // system of ODE
    1159                                  &tt,                   // from t...
    1160                                   B_.step_,             // ...to t=t+h
    1161                                  &B_.IntegrationStep_ , // integration window (written on!)
    1162                                   S_.y_);               // neuron state
    1163 
    1164           if ( status != GSL_SUCCESS )
    1165             throw GSLSolverFailure(get_name(), status);
    1166         }
    1167 EOF
    1168 vi)))
    1169 
    1170     (let ((isyns (lookup-def 'i-synapses synapse-info))
    1171           (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
    1172 
    1173       (for-each (lambda (isyn psc)
    1174                   (let ((gi (lookup-def (nest-name (fourth isyn)) state-index-map)))
    1175                     (if gi
    1176                         (pp indent+ (,(sprintf "      S_.y_[~A] += B_.spike_~A.get_value(lag);"
    1177                                                gi (second psc)))
    1178                         ))))
    1179                 isyns pscs))
    1180 
    1181     (pp indent+ (,(sprintf #<<EOF
    1182         // sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
    1183         // refractory?
    1184         if (S_.r_)
    1185           {
    1186             --S_.r_;
    1187           }
    1188         else
    1189           {
    1190            ~A
    1191           }
    1192    
    1193         // set new input current
    1194         B_.I_stim_ = B_.currents_.get_value(lag);
    1195 
    1196         // log state data
    1197         B_.logger_.record_data(origin.get_steps() + lag);
    1198 
    1199       }
    1200 EOF
    1201   (if (lookup-def 'V_t const-defs) threshold "")
    1202   )))
    1203  
    1204   (pp indent  #\})
    1205 
    1206 ))
    1207 
    1208 
    1209 (define (output-handle
    1210          sysname state-index-map steady-state-index-map
    1211          const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1212          reaction-eq-defs i-eqs pool-ions perm-ions
    1213          synapse-info
    1214          indent indent+)
    1215  
    1216   (pp indent  ,nl (,(sprintf #<<EOF
    1217 void ~A::handle(SpikeEvent & e)
    1218   {
    1219     int flag;
    1220     assert(e.get_delay() > 0);
    1221     flag = 0;
    1222 EOF
    1223    sysname)))
    1224 
    1225     (let ((isyns (lookup-def 'i-synapses synapse-info))
    1226           (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
    1227       (for-each
    1228        (lambda (isyn psc)
    1229          (pp indent+
    1230              (,(sprintf #<<EOF
    1231     if ((flag==0) && (e.get_weight() > P_.~A))
    1232       {
    1233         B_.spike_~A.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
    1234                               abs(e.get_weight()) * e.get_multiplicity());
    1235         flag = 1;
    1236       }
    1237 EOF
    1238              (nest-name (second isyn))
    1239              (second psc)))))
    1240        isyns pscs))
    1241 
    1242   (pp indent ,nl #\} ,nl)
    1243 
    1244   (pp indent  ,nl (,(sprintf #<<EOF
    1245 void ~A::handle(CurrentEvent& e)
    1246   {
    1247     assert(e.get_delay() > 0);
    1248 
    1249     const double_t c=e.get_current();
    1250     const double_t w=e.get_weight();
    1251 
    1252     B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
    1253                         w *c);
    1254   }
    1255 
    1256 void ~A::handle(DataLoggingRequest& e)
    1257   {
    1258     B_.logger_.handle(e);
    1259   }
    1260 EOF
    1261   sysname sysname)))
    1262 
    1263 )
    1264 
    1265 
    1266 (define (output-recordables
    1267          sysname state-index-map steady-state-index-map
    1268          const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1269          reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
    1270  
    1271   (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
    1272   (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
    1273   (pp indent #\{)
    1274   (for-each
    1275    (lambda (def)
    1276      (let ((n (first def)) (i (second def)))
    1277        (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);"
    1278                               (->string n) sysname sysname (string-upcase (->string n)))))
    1279        ))
    1280    state-index-map)
    1281 
    1282   (if (lookup-def 'v state-index-map)
    1283       (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);"
    1284                              sysname sysname (string-upcase (->string 'v))))))
    1285 
    1286   (pp indent #\})
    1287   )
    1288 
    1289 
    1290 (define (output-prelude sysname steady-state-index-map indent)
    1291   (pp indent (,(sprintf "#include \"~A.h\"
    1292 #include \"exceptions.h\"
    1293 #include \"network.h\"
    1294 #include \"dict.h\"
    1295 #include \"integerdatum.h\"
    1296 #include \"doubledatum.h\"
    1297 #include \"dictutils.h\"
    1298 #include \"numerics.h\"
    1299 #include <limits>
    1300 
    1301 #include \"universal_data_logger_impl.h\"
    1302 
    1303 #include <iomanip>
    1304 #include <iostream>
    1305 #include <cstdio>
    1306 #include <cmath>
    1307 
    1308 #include <gsl/gsl_errno.h>
    1309 #include <gsl/gsl_matrix.h>
    1310 #include <gsl/gsl_sf_exp.h>
    1311 "  sysname)))
    1312   (if (not (null? steady-state-index-map))
    1313       (pp indent ("#include <gsl/gsl_multimin.h>")))
    1314   )
    1315 
    1316 
    1317 (define (output-header
    1318          sysname state-index-map steady-state-index-map
    1319          const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1320          reaction-eq-defs i-eqs pool-ions perm-ions
    1321          synapse-info
    1322          indent indent+)
    1323 
    1324   (pp indent ("#include \"nest.h\"
    1325 #include \"event.h\"
    1326 #include \"archiving_node.h\"
    1327 #include \"ring_buffer.h\"
    1328 #include \"connection.h\"
    1329 #include \"universal_data_logger.h\"
    1330 #include \"recordables_map.h\"
    1331 #include <gsl/gsl_odeiv.h>
    1332 "))
    1333  
    1334   (pp indent "namespace nest {" ,nl)
    1335 
    1336   (pp indent (extern "\"C\""
    1337                      int ,(s+ sysname "_dynamics")
    1338                      "(double, const double*, double*, void*)" #\;) ,nl)
    1339 
    1340   (pp indent
    1341       (,(sprintf "class ~A:" sysname)
    1342        "public Archiving_Node { "))
    1343   (pp indent+
    1344        ("public: ")
    1345        (,(s+ "~" sysname "();"))
    1346        (,(s+ sysname "(const " sysname "&);"))
    1347        (,(s+ sysname "();")))
    1348 
    1349   (pp indent (#<<EOF
    1350     using Node::connect_sender;
    1351     using Node::handle;
    1352 
    1353     port check_connection(Connection&, port);
    1354    
    1355     void handle(SpikeEvent &);
    1356     void handle(CurrentEvent &);
    1357     void handle(DataLoggingRequest &);
    1358    
    1359     port connect_sender(SpikeEvent &, port);
    1360     port connect_sender(CurrentEvent &, port);
    1361     port connect_sender(DataLoggingRequest &, port);
    1362    
    1363     void get_status(DictionaryDatum &) const;
    1364     void set_status(const DictionaryDatum &);
    1365    
    1366     void init_node_(const Node& proto);
    1367     void init_state_(const Node& proto);
    1368     void init_buffers_();
    1369     void calibrate();
    1370    
    1371     void update(Time const &, const long_t, const long_t);
    1372 EOF
    1373 ))
    1374 
    1375   (pp indent (,(sprintf #<<EOF
    1376     // make dynamics function quasi-member
    1377     friend int ~A_dynamics(double, const double*, double*, void*);
    1378 
    1379     // The next two classes need to be friends to access the State_ class/member
    1380     friend class RecordablesMap<~A>;
    1381     friend class UniversalDataLogger<~A>;
    1382 EOF
    1383 sysname sysname sysname)))
    1384 
    1385  
    1386   (pp indent+ ,nl "struct Parameters_ { ")
    1387 
    1388   (for-each
    1389    (lambda (x) (pp indent+ (,(sprintf "double ~A;" (first x)))))
    1390    const-defs)
    1391 
    1392   (pp indent+
    1393       ("Parameters_();")
    1394       ("void get(DictionaryDatum&) const;")
    1395       ("void set(const DictionaryDatum&);"))
    1396 
    1397   (pp indent+ "};")
    1398 
    1399   (pp indent+ ,nl "struct State_ { ")
    1400 
    1401   (pp indent+ ,nl
    1402       "enum StateVecElems {"
    1403       (,(slp ", " (map (lambda (x)
    1404                          (let ((n (string-upcase (->string (first x))))
    1405                                (ni (second x)))
    1406                            (s+ n " = " ni)))
    1407                        state-index-map)))
    1408       "};"
    1409       (,(sprintf "double y_[~A];" (length state-index-map)))
    1410       "int_t     r_;"
    1411       "State_(const Parameters_& p);"
    1412       "State_(const State_& s);"
    1413      
    1414       "State_& operator=(const State_& s);"
    1415 
    1416       "void get(DictionaryDatum&) const;"
    1417       "void set(const DictionaryDatum&, const Parameters_&);")
    1418   (pp indent+ "};")
    1419 
    1420 
    1421   (pp indent+ ,nl #<<EOF
    1422     struct Variables_ {
    1423       int_t    RefractoryCounts_;
    1424       double   U_old_; // for spike-detection
    1425     };
    1426 EOF
    1427 )
    1428 
    1429   (pp indent+ ,nl "struct Buffers_ { ")
    1430 
    1431   (pp indent+ ,nl
    1432       (,(sprintf "Buffers_(~A&);" sysname))
    1433       (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
    1434       (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))
    1435 
    1436   (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
    1437     (for-each
    1438      (lambda (psc)
    1439        (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
    1440      pscs))
    1441 
    1442   (pp indent+ ,nl
    1443 #<<EOF
    1444 
    1445   RingBuffer currents_;
    1446 
    1447   gsl_odeiv_step*    s_;    //!< stepping function
    1448   gsl_odeiv_control* c_;    //!< adaptive stepsize control function
    1449   gsl_odeiv_evolve*  e_;    //!< evolution function
    1450   gsl_odeiv_system   sys_;  //!< struct describing system
    1451 
    1452   double_t step_;           //!< step size in ms
    1453   double   IntegrationStep_;//!< current integration time step, updated by GSL
    1454 
    1455   /**
    1456    * Input current injected by CurrentEvent.
    1457    * This variable is used to transport the current applied into the
    1458    * _dynamics function computing the derivative of the state vector.
    1459    * It must be a part of Buffers_, since it is initialized once before
    1460    * the first simulation, but not modified before later Simulate calls.
    1461    */
    1462   double_t I_stim_;
    1463 EOF
    1464 )
    1465   (pp indent+ "};")
    1466 
    1467   (pp indent+
    1468       "template <State_::StateVecElems elem>"
    1469       "double_t get_y_elem_() const { return S_.y_[elem]; }"
    1470       "Parameters_ P_;"
    1471       "State_      S_;"
    1472       "Variables_  V_;"
    1473       "Buffers_    B_;"
    1474       (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
    1475       "}; ")
    1476 
    1477   (pp indent+ (,(sprintf #<<EOF
    1478   inline
    1479   port ~A::check_connection(Connection& c, port receptor_type)
    1480   {
    1481     SpikeEvent e;
    1482     e.set_sender(*this);
    1483     c.check_event(e);
    1484     return c.get_target()->connect_sender(e, receptor_type);
    1485   }
    1486 
    1487   inline
    1488   port ~A::connect_sender(SpikeEvent&, port receptor_type)
    1489   {
    1490     if (receptor_type != 0)
    1491       throw UnknownReceptorType(receptor_type, get_name());
    1492     return 0;
    1493   }
    1494  
    1495   inline
    1496   port ~A::connect_sender(CurrentEvent&, port receptor_type)
    1497   {
    1498     if (receptor_type != 0)
    1499       throw UnknownReceptorType(receptor_type, get_name());
    1500     return 0;
    1501   }
    1502 
    1503   inline
    1504   port ~A::connect_sender(DataLoggingRequest& dlr,
    1505                                       port receptor_type)
    1506   {
    1507     if (receptor_type != 0)
    1508       throw UnknownReceptorType(receptor_type, get_name());
    1509     return B_.logger_.connect_logging_device(dlr, recordablesMap_);
    1510   }
    1511 
    1512   inline
    1513     void ~A::get_status(DictionaryDatum &d) const
    1514   {
    1515     P_.get(d);
    1516     S_.get(d);
    1517     Archiving_Node::get_status(d);
    1518 
    1519     (*d)[names::recordables] = recordablesMap_.get_list();
    1520 
    1521     def<double_t>(d, names::t_spike, get_spiketime_ms());
    1522   }
    1523 
    1524   inline
    1525     void ~A::set_status(const DictionaryDatum &d)
    1526   {
    1527     Parameters_ ptmp = P_;  // temporary copy in case of errors
    1528     ptmp.set(d);                       // throws if BadProperty
    1529     State_      stmp = S_;  // temporary copy in case of errors
    1530     stmp.set(d, ptmp);                 // throws if BadProperty
    1531 
    1532     // We now know that (ptmp, stmp) are consistent. We do not
    1533     // write them back to (P_, S_) before we are also sure that
    1534     // the properties to be set in the parent class are internally
    1535     // consistent.
    1536     Archiving_Node::set_status(d);
    1537 
    1538     // if we get here, temporaries contain consistent set of properties
    1539     P_ = ptmp;
    1540     S_ = stmp;
    1541 
    1542     calibrate();
    1543   }
    1544 EOF
    1545   sysname sysname sysname
    1546   sysname sysname sysname)))
    1547 
    1548   (pp indent "}" ,nl)
    1549   )
    1550484
    1551485(define (nemo:nest-translator sys . rest)
     
    1554488  (define (cn x)   (first x))
    1555489
    1556   (let-optionals rest ((dirname "."))
     490  (let-optionals rest ((dirname ".") (method 'gsl))
     491
     492    (if (not (member method '(gsl cvode)))
     493        (nemo:error 'nemo:nest-translator ": unknown method " method))
     494
     495
    1557496  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
    1558497    (let ((imports  ((dis 'imports)  sys))
     
    1861800             (with-output-to-port cpp-output
    1862801               (lambda ()
    1863                  (output-dy sysname const-defs state-index-map
     802                 (output-dy sysname method const-defs state-index-map
    1864803                            external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
    1865804                            pool-ions mcap i-eqs v-eq
     
    1910849               (lambda ()
    1911850                 (output-buffers+nodes
    1912                   sysname state-index-map steady-state-index-map
     851                  sysname method state-index-map steady-state-index-map
    1913852                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1914853                  reaction-eq-defs i-eqs pool-ions perm-ions
     
    1922861               (lambda ()
    1923862                 (output-update
    1924                   sysname state-index-map steady-state-index-map
     863                  sysname method state-index-map steady-state-index-map
    1925864                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1926865                  reaction-eq-defs i-eqs pool-ions perm-ions
     
    1948887               (lambda ()
    1949888                 (output-header
    1950                   sysname state-index-map steady-state-index-map
     889                  sysname method state-index-map steady-state-index-map
    1951890                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1952891                  reaction-eq-defs i-eqs pool-ions perm-ions
     
    1964903
    1965904
     905(define (output-dy sysname method const-defs state-index-map
     906                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
     907                   pool-ions mcap i-eqs v-eq indent indent+)
     908
     909
     910  (case method
     911    ((cvode)
     912     (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
     913                     (,(slp ", " `("double t"
     914                                   "N_Vector y"
     915                                   "N_Vector f"
     916                                   "void* pnode"
     917                                   )))
     918                     #\{
     919                     )))
     920    ((gsl)
     921     (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
     922                     (,(slp ", " `("double t"
     923                                   "const double y[]"
     924                                   "double f[]"
     925                                   "void* pnode"
     926                                   )))
     927                     #\{
     928                     )))
     929     )
     930
     931  (let* (
     932         (asgn-eq-defs
     933          (map
     934           (lambda (def) (list (first def)
     935                               (add-params-to-fncall
     936                                (canonicalize-expr/C++ (second def)))))
     937           asgn-eq-defs))
     938
     939         (rate-eq-defs
     940          (map
     941           (lambda (def) (list (first def)
     942                               (add-params-to-fncall
     943                                (canonicalize-expr/C++ (second def)))))
     944           rate-eq-defs))
     945
     946         (reaction-eq-defs
     947          (map
     948           (lambda (def) (list (first def)
     949                               (add-params-to-fncall
     950                                (canonicalize-expr/C++ (second def)))))
     951           reaction-eq-defs))
     952
     953         (i-eqs
     954          (map
     955           (lambda (def) (list (first def)
     956                               (add-params-to-fncall (canonicalize-expr/C++ (second def)))))
     957           i-eqs))
     958
     959         (v-eq
     960          (and v-eq
     961               (list (first v-eq)
     962                     (add-params-to-fncall (canonicalize-expr/C++ (second v-eq))))))
     963
     964         (eqs
     965          (append
     966           
     967           asgn-eq-defs
     968           reaction-eq-defs
     969                 
     970           (map (lambda (pool-ion)
     971                  (let ((n (third pool-ion))
     972                        (b (first pool-ion)))
     973                    (list n b)))
     974                pool-ions)))
     975         
     976         (eq-dag
     977          (map (lambda (def)
     978                 (cons (first def) (enum-freevars (second def) '() '())))
     979               eqs))
     980
     981         (eq-order
     982          (reverse
     983           (topological-sort eq-dag
     984                             (lambda (x y) (string=? (->string x) (->string y))))))
     985
     986         (eq-locals  (find-locals
     987                      (map second
     988                           (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
     989                               (append i-eqs rate-eq-defs eqs)))))
     990         )
     991
     992
     993
     994    (pp indent+ (double ,(slp ", " (delete-duplicates
     995                                    (map ->string
     996                                         (filter (lambda (x) (not (member x nest-builtin-consts)))
     997                                                 (append
     998                                                  eq-locals
     999                                                  eq-order
     1000                                                  (map first external-eq-defs)
     1001                                                  (map first state-index-map)
     1002                                                  (map first i-eqs)
     1003                                                  (map first const-defs))))
     1004                                    string=?))
     1005                        ";"))
     1006
     1007    (pp indent+ ,nl
     1008        ("// S is shorthand for the type that describes the model state ")
     1009        (typedef ,(s+ sysname "::State_") "S;")
     1010       
     1011        ,nl
     1012       
     1013        ("// cast the node ptr to an object of the proper type")
     1014        ("assert(pnode);")
     1015        ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
     1016       
     1017        ,nl
     1018
     1019        ("// params is a reference to the model parameters ")
     1020        (const struct ,(s+ sysname "::Parameters_") "*params;")
     1021        ("params = &(node.P_);")
     1022       
     1023        ,nl
     1024
     1025        ("// y[] must be the state vector supplied by the integrator, ")
     1026        ("// not the state vector in the node, node.S_.y[]. ")
     1027        ,nl
     1028        )
     1029
     1030    (for-each (lambda (def)
     1031                (let ((n (first def))
     1032                      (b (second def)))
     1033                  (pp indent+ ,(expr->string/C++ b (nest-name n)))))
     1034              external-eq-defs)
     1035   
     1036    (for-each (lambda (def)
     1037                (let ((n (first def)) )
     1038                  (pp indent+
     1039                      ,(expr->string/C++
     1040                        (s+ "params->" n) n))))
     1041              const-defs)
     1042
     1043    (let ((vi (lookup-def 'v state-index-map)))
     1044      (if vi (pp indent+ ,(expr->string/C++ (s+ "y[" vi "]") 'v))))
     1045
     1046    (for-each (lambda (def)
     1047                (let ((n (first def)) )
     1048                  (pp indent+
     1049                      ,(expr->string/C++
     1050                        (s+ "y[" (lookup-def n state-index-map) "]") n))))
     1051              rate-eq-defs)
     1052
     1053    (for-each (lambda (n)
     1054                (let ((b  (lookup-def n eqs)))
     1055                  (if b (pp indent+
     1056                            ,(expr->string/C++ b (nest-name n))))))
     1057              eq-order)
     1058
     1059    (for-each (lambda (def)
     1060                (let ((n (s+ "f[" (lookup-def (first def) state-index-map) "]"))
     1061                      (b (second def)))
     1062                  (pp indent+ ,(s+ "// rate equation " (first def))
     1063                      ,(expr->string/C++ b n))))
     1064              rate-eq-defs)
     1065   
     1066    (for-each (lambda (def)
     1067                (pp indent+ ,(expr->string/C++ (second def) (first def))))
     1068            i-eqs)
     1069   
     1070    (if v-eq
     1071        (let ((vi (lookup-def 'v state-index-map)))
     1072          (if vi
     1073              (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "f[" vi "]")))
     1074              )))
     1075
     1076    (case method
     1077      ((cvode)
     1078        (pp indent+ ,nl ("return 0;")))
     1079      ((gsl)
     1080        (pp indent+ ,nl ("return GSL_SUCCESS;")))
     1081      )
     1082    (pp indent  #\})
     1083
     1084    ))
     1085
     1086
     1087
     1088(define (output-parameters sysname const-defs indent indent+)
     1089
     1090  (let* ((parameter-defs
     1091          (map
     1092           (lambda (def) (list (first def)
     1093                               (add-params-to-fncall (canonicalize-expr/C++ (second def)))))
     1094           const-defs))
     1095         (parameter-locals  (find-locals (map second parameter-defs))))
     1096
     1097    (pp indent ,nl (,(s+ sysname "::Parameters_::Parameters_") () ":"))
     1098   
     1099    (let ((const-exprs
     1100           (slp ",\n"
     1101                (map (lambda (def)
     1102                       (let* ((n  (first def)) (b (second def)))
     1103                         (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
     1104                     const-defs) )))
     1105     
     1106      (if (not (null? parameter-locals))
     1107          (pp indent+ (double ,(slp ", " parameter-locals) ";")))
     1108
     1109      (pp indent+ ,const-exprs)
     1110     
     1111      (pp indent ("{}"))
     1112     
     1113      )))
     1114
     1115
     1116(define (output-init sysname state-index-map steady-state-index-map
     1117                     external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     1118                     reaction-eq-defs i-eqs v-eq pool-ions perm-ions indent indent+)
     1119
     1120 
     1121  (let* ((N (length state-index-map))
     1122
     1123         (asgn-eq-defs
     1124          (map
     1125           (lambda (def) (list (first def)
     1126                               (add-params-to-fncall
     1127                                (canonicalize-expr/C++ (second def)))))
     1128           asgn-eq-defs))
     1129         
     1130         (init-eq-defs
     1131          (map
     1132           (lambda (def) (list (first def)
     1133                               (add-params-to-fncall
     1134                                (canonicalize-expr/C++ (second def)))))
     1135           init-eq-defs))
     1136
     1137         (reaction-eq-defs
     1138          (map
     1139           (lambda (def) (list (first def)
     1140                               (add-params-to-fncall
     1141                                (canonicalize-expr/C++ (second def)))))
     1142           reaction-eq-defs))
     1143
     1144         (i-eqs
     1145          (map
     1146           (lambda (def) (list (first def)
     1147                               (add-params-to-fncall
     1148                                (canonicalize-expr/C++ (second def)))))
     1149           i-eqs))
     1150
     1151         (v-eq
     1152          (and v-eq
     1153               (list (first v-eq)
     1154                     (add-params-to-fncall
     1155                      (canonicalize-expr/C++ (second v-eq))))))
     1156         
     1157         (init-eqs
     1158          (append
     1159           
     1160           asgn-eq-defs
     1161           init-eq-defs
     1162
     1163           (map (lambda (pool-ion)
     1164                  (let ((n (third pool-ion))
     1165                        (b (first pool-ion)))
     1166                    (list n b)))
     1167                pool-ions)))
     1168
     1169         (init-dag
     1170          (map (lambda (def)
     1171                 (cons (first def) (enum-freevars (second def) '() '())))
     1172               init-eqs))
     1173
     1174         (init-order
     1175          (reverse
     1176           (topological-sort init-dag
     1177                             (lambda (x y) (string=? (->string x) (->string y))))))
     1178
     1179         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
     1180         )
     1181   
     1182    (pp indent ,nl (,(s+ sysname "::State_::State_") (const Parameters_& p) ": r_(0)"))
     1183   
     1184    (pp indent  #\{)
     1185
     1186    (pp indent+
     1187        (double ,(slp ", " (delete-duplicates
     1188                            (map ->string
     1189                                 (filter (lambda (x) (not (member x nest-builtin-consts)))
     1190                                         (append
     1191                                          init-locals
     1192                                          init-order
     1193                                          (map first external-eq-defs)
     1194                                          (map first state-index-map)
     1195                                          (map first i-eqs)
     1196                                          (map first const-defs)
     1197                                          (map first reaction-eq-defs)
     1198                                          )))
     1199                                 string=?))
     1200                      ";")
     1201        ("const Parameters_ *params = &p;")
     1202        )
     1203
     1204    (for-each (lambda (def)
     1205                (let ((n (first def)) )
     1206                  (pp indent+
     1207                      ,(expr->string/C++
     1208                        (s+ "p." n) n))))
     1209              const-defs)
     1210
     1211    (let ((vi (lookup-def 'v state-index-map))
     1212          (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
     1213      (if (and vi vrest)
     1214          (pp indent+ ,(expr->string/C++  vrest 'v ))))
     1215
     1216    (for-each (lambda (def)
     1217                (let ((n (first def))
     1218                      (b (second def)))
     1219                  (pp indent+ ,(expr->string/C++ b (nest-name n)))))
     1220              external-eq-defs)
     1221
     1222    (for-each (lambda (n)
     1223                (let ((b  (lookup-def n init-eqs)))
     1224                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
     1225              init-order)
     1226   
     1227    (for-each (lambda (def)
     1228                (let* ((n  (first def))
     1229                       (ni (lookup-def n state-index-map)))
     1230                  (if ni (pp indent+ ,(expr->string/C++ n  (s+ "y_[" ni "]"))))))
     1231              init-eq-defs)
     1232#|
     1233    (if (not (null? steady-state-index-map))
     1234        (begin
     1235          (gsl-fminimizer sysname N indent+)
     1236          (for-each
     1237           (lambda (def)
     1238             (let* ((n (first def))
     1239                    (si (lookup-def n steady-state-index-map))
     1240                    (ni (lookup-def n state-index-map)))
     1241               (if si (begin
     1242                        (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n))
     1243                        (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")")
     1244                                                       (s+ "y_[" ni "]")))))))
     1245           rate-eq-defs)
     1246      (pp indent+ "gsl_vector_free (ys);")
     1247      ))
     1248|#
     1249
     1250    (for-each
     1251     (lambda (def)
     1252       (let ((n (first def)) (b (second def)))
     1253         (if (not (lookup-def n init-eq-defs))
     1254             (pp indent+ ,(expr->string/C++ b n)))))
     1255     reaction-eq-defs)
     1256   
     1257    (for-each
     1258     (lambda (def)
     1259       (pp indent+ ,(expr->string/C++ (second def) (first def))))
     1260     i-eqs)
     1261
     1262    (if v-eq
     1263        (let ((vi (lookup-def 'v state-index-map)))
     1264          (if vi
     1265              (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "y_[" vi "]"))))))
     1266     
     1267    (pp indent  #\})
     1268   
     1269    (pp indent ,nl (,(s+ sysname "::State_::State_") (const State_& s) ": r_(s.r_)"))
     1270    (pp indent  #\{)
     1271    (pp indent+ (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N)))
     1272    (pp indent  #\})
     1273
     1274    (pp indent ,nl (,(s+ sysname "::State_& " sysname "::State_::operator=(const State_& s)")))
     1275
     1276    (pp indent  #\{)
     1277
     1278
     1279    (pp indent+ (,(sprintf #<<EOF
     1280   assert(this != &s); 
     1281     for ( size_t i = 0 ; i < ~A ; ++i )
     1282       y_[i] = s.y_[i];
     1283     r_ = s.r_;
     1284     return *this;
     1285EOF
     1286N))))
     1287
     1288     (pp indent  #\})
     1289
     1290))
     1291
     1292
     1293
     1294(define (gsl-fminimizer sysname N indent)
     1295
     1296  (pp indent (,(sprintf #<<EOF
     1297
     1298  int minimizer_status, minimizer_iterc;
     1299
     1300  const gsl_multimin_fdfminimizer_type *T;
     1301  gsl_multimin_fdfminimizer *S;
     1302
     1303  gsl_vector *ys;
     1304  gsl_multimin_function_fdf F;
     1305
     1306  F.f   = &~A_dynamics;
     1307  F.df  = &~A_dynamics;
     1308  F.fdf = &~A_dynamics;
     1309  F.n   = ~A;
     1310  F.params = (void *)&p;
     1311
     1312  ys = gsl_vector_alloc (~A);
     1313
     1314  for (int i = 0; i < ~A; i++)
     1315  {
     1316     gsl_vector_set (ys, i, 0.0);
     1317  }
     1318
     1319  T = gsl_multimin_fdfminimizer_conjugate_fr;
     1320  S = gsl_multimin_fdfminimizer_alloc (T, ~A);
     1321
     1322  gsl_multimin_fdfminimizer_set (S, &F, ys, 0.01, 1e-4);
     1323
     1324  do
     1325    {
     1326      minimizer_iterc++;
     1327      minimizer_status = gsl_multimin_fdfminimizer_iterate (S);
     1328
     1329      if (minimizer_status)
     1330        break;
     1331
     1332      minimizer_status = gsl_multimin_test_gradient (S->gradient, 1e-3);
     1333
     1334    }
     1335  while (minimizer_status == GSL_CONTINUE && minimizer_iterc < 100);
     1336
     1337  gsl_multimin_fdfminimizer_free (S);
     1338
     1339EOF
     1340sysname sysname sysname N N N N)))
    19661341)
     1342
     1343
     1344
     1345(define (output-accessors+modifiers
     1346         sysname state-index-map const-defs asgn-eq-defs rate-eq-defs
     1347         reaction-eq-defs i-eqs v-eq pool-ions perm-ions indent indent+)
     1348
     1349  (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
     1350  (pp indent  #\{)
     1351
     1352  (for-each
     1353   (lambda (def)
     1354     (let ((n (first def)))
     1355       (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
     1356   const-defs)
     1357
     1358  (pp indent  #\})
     1359
     1360  (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
     1361  (pp indent  #\{)
     1362
     1363  (for-each
     1364   (lambda (def)
     1365     (let ((n (first def)))
     1366       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))))
     1367   const-defs)
     1368
     1369  (pp indent  #\})
     1370
     1371  (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
     1372  (pp indent  #\{)
     1373
     1374  (for-each
     1375   (lambda (def)
     1376     (let ((n (first def)) (i (second def)))
     1377       (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
     1378   state-index-map)
     1379
     1380  (pp indent  #\})
     1381
     1382  (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
     1383  (pp indent  #\{)
     1384
     1385  (for-each
     1386   (lambda (def)
     1387     (let ((n (first def)) (i (second def)))
     1388       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
     1389   state-index-map)
     1390
     1391
     1392  (pp indent  #\})
     1393
     1394)
     1395
     1396
     1397
     1398(define (output-handle
     1399         sysname state-index-map steady-state-index-map
     1400         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     1401         reaction-eq-defs i-eqs pool-ions perm-ions
     1402         synapse-info
     1403         indent indent+)
     1404 
     1405  (pp indent  ,nl (,(sprintf #<<EOF
     1406void ~A::handle(SpikeEvent & e)
     1407  {
     1408    int flag;
     1409    assert(e.get_delay() > 0);
     1410    flag = 0;
     1411EOF
     1412   sysname)))
     1413
     1414   (let ((isyns (lookup-def 'i-synapses synapse-info))
     1415          (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
     1416      (for-each
     1417       (lambda (isyn psc)
     1418         (pp indent+
     1419             (,(sprintf #<<EOF
     1420    if ((flag==0) && (e.get_weight() > P_.~A))
     1421      {
     1422        B_.spike_~A.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
     1423                              abs(e.get_weight()) * e.get_multiplicity());
     1424        flag = 1;
     1425      }
     1426EOF
     1427             (nest-name (second isyn))
     1428             (second psc)))))
     1429       isyns pscs))
     1430
     1431  (pp indent ,nl #\} ,nl)
     1432
     1433  (pp indent  ,nl (,(sprintf #<<EOF
     1434void ~A::handle(CurrentEvent& e)
     1435  {
     1436    assert(e.get_delay() > 0);
     1437
     1438    const double_t c=e.get_current();
     1439    const double_t w=e.get_weight();
     1440
     1441    B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
     1442                        w *c);
     1443  }
     1444
     1445void ~A::handle(DataLoggingRequest& e)
     1446  {
     1447    B_.logger_.handle(e);
     1448  }
     1449EOF
     1450  sysname sysname)))
     1451
     1452)
     1453
     1454
     1455(define (output-recordables
     1456         sysname state-index-map steady-state-index-map
     1457         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     1458         reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
     1459 
     1460  (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
     1461  (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
     1462  (pp indent #\{)
     1463  (for-each
     1464   (lambda (def)
     1465     (let ((n (first def)) (i (second def)))
     1466       (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);"
     1467                              (->string n) sysname sysname (string-upcase (->string n)))))
     1468       ))
     1469   state-index-map)
     1470
     1471  (if (lookup-def 'v state-index-map)
     1472      (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);"
     1473                             sysname sysname (string-upcase (->string 'v))))))
     1474
     1475  (pp indent #\})
     1476  )
     1477
     1478
     1479(define (output-prelude sysname steady-state-index-map indent)
     1480  (pp indent (,(sprintf "#include \"~A.h\"
     1481#include \"exceptions.h\"
     1482#include \"network.h\"
     1483#include \"dict.h\"
     1484#include \"integerdatum.h\"
     1485#include \"doubledatum.h\"
     1486#include \"dictutils.h\"
     1487#include \"numerics.h\"
     1488#include <limits>
     1489
     1490#include \"universal_data_logger_impl.h\"
     1491
     1492#include <iomanip>
     1493#include <iostream>
     1494#include <cstdio>
     1495#include <cmath>
     1496"  sysname)))
     1497  )
     1498
     1499(define (output-header
     1500         sysname method state-index-map steady-state-index-map
     1501         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     1502         reaction-eq-defs i-eqs pool-ions perm-ions
     1503         synapse-info
     1504         indent indent+)
     1505
     1506  (pp indent ("#include \"nest.h\"
     1507#include \"event.h\"
     1508#include \"archiving_node.h\"
     1509#include \"ring_buffer.h\"
     1510#include \"connection.h\"
     1511#include \"universal_data_logger.h\"
     1512#include \"recordables_map.h\"
     1513"))
     1514
     1515  (case method
     1516    ((cvode)
     1517     (pp indent
     1518         ("#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */")
     1519         ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
     1520         ("#include <cvode/cvode_dense.h>       /* prototype for CVDense */")
     1521         ("#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */")
     1522         ("#include <sundials/sundials_types.h> /* definition of type realtype */")
     1523         ("#define Ith(v,i)    NV_Ith_S(v,i-1)       /* Ith numbers components 1..NEQ */")))
     1524    ((gsl)
     1525     (pp indent
     1526         ("#include <gsl/gsl_errno.h>")
     1527         ("#include <gsl/gsl_matrix.h>")
     1528         ("#include <gsl/gsl_sf_exp.h>")
     1529         ("#include <gsl/gsl_odeiv.h>")))
     1530    )
     1531
     1532  (pp indent "namespace nest {" ,nl)
     1533
     1534  (pp indent (extern "\"C\""
     1535                     int ,(s+ sysname "_dynamics")
     1536                     "(double, const double*, double*, void*)" #\;) ,nl)
     1537
     1538  (pp indent
     1539      (,(sprintf "class ~A:" sysname)
     1540       "public Archiving_Node { "))
     1541  (pp indent+
     1542       ("public: ")
     1543       (,(s+ "~" sysname "();"))
     1544       (,(s+ sysname "(const " sysname "&);"))
     1545       (,(s+ sysname "();")))
     1546
     1547  (pp indent (#<<EOF
     1548    using Node::connect_sender;
     1549    using Node::handle;
     1550
     1551    port check_connection(Connection&, port);
     1552   
     1553    void handle(SpikeEvent &);
     1554    void handle(CurrentEvent &);
     1555    void handle(DataLoggingRequest &);
     1556   
     1557    port connect_sender(SpikeEvent &, port);
     1558    port connect_sender(CurrentEvent &, port);
     1559    port connect_sender(DataLoggingRequest &, port);
     1560   
     1561    void get_status(DictionaryDatum &) const;
     1562    void set_status(const DictionaryDatum &);
     1563   
     1564    void init_node_(const Node& proto);
     1565    void init_state_(const Node& proto);
     1566    void init_buffers_();
     1567    void calibrate();
     1568   
     1569    void update(Time const &, const long_t, const long_t);
     1570EOF
     1571))
     1572
     1573  (pp indent (,(sprintf #<<EOF
     1574    // make dynamics function quasi-member
     1575    friend int ~A_dynamics(double, const double*, double*, void*);
     1576
     1577    // The next two classes need to be friends to access the State_ class/member
     1578    friend class RecordablesMap<~A>;
     1579    friend class UniversalDataLogger<~A>;
     1580EOF
     1581sysname sysname sysname)))
     1582
     1583 
     1584  (pp indent+ ,nl "struct Parameters_ { ")
     1585
     1586  (for-each
     1587   (lambda (x) (pp indent+ (,(sprintf "double ~A;" (first x)))))
     1588   const-defs)
     1589
     1590  (pp indent+
     1591      ("Parameters_();")
     1592      ("void get(DictionaryDatum&) const;")
     1593      ("void set(const DictionaryDatum&);"))
     1594
     1595  (pp indent+ "};")
     1596
     1597  (pp indent+ ,nl "struct State_ { ")
     1598
     1599  (pp indent+ ,nl
     1600      "enum StateVecElems {"
     1601      (,(slp ", " (map (lambda (x)
     1602                         (let ((n (string-upcase (->string (first x))))
     1603                               (ni (second x)))
     1604                           (s+ n " = " ni)))
     1605                       state-index-map)))
     1606      "};"
     1607      (,(sprintf "double y_[~A];" (length state-index-map)))
     1608      "int_t     r_;"
     1609      "State_(const Parameters_& p);"
     1610      "State_(const State_& s);"
     1611     
     1612      "State_& operator=(const State_& s);"
     1613
     1614      "void get(DictionaryDatum&) const;"
     1615      "void set(const DictionaryDatum&, const Parameters_&);")
     1616  (pp indent+ "};")
     1617
     1618
     1619  (pp indent+ ,nl #<<EOF
     1620    struct Variables_ {
     1621      int_t    RefractoryCounts_;
     1622      double   U_old_; // for spike-detection
     1623    };
     1624EOF
     1625)
     1626
     1627  (pp indent+ ,nl "struct Buffers_ { ")
     1628
     1629  (pp indent+ ,nl
     1630      (,(sprintf "Buffers_(~A&);" sysname))
     1631      (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
     1632      (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))
     1633
     1634  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
     1635    (for-each
     1636     (lambda (psc)
     1637       (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
     1638     pscs))
     1639
     1640  (pp indent+ ,nl
     1641  ,(case method
     1642    ((cvode)
     1643#<<EOF
     1644  N_Vector y;    //!< current state vector used by CVode
     1645  void *   sys_;  //!< CVode control structure
     1646EOF
     1647))
     1648    ((gsl)
     1649     (pp indent+ ,nl
     1650#<<EOF
     1651  gsl_odeiv_step*    s_;    //!< stepping function
     1652  gsl_odeiv_control* c_;    //!< adaptive stepsize control function
     1653  gsl_odeiv_evolve*  e_;    //!< evolution function
     1654  gsl_odeiv_system   sys_;  //!< struct describing system
     1655EOF
     1656))
     1657
     1658    )
     1659
     1660  (pp indent+ ,nl
     1661#<<EOF
     1662  RingBuffer currents_;
     1663
     1664  double_t step_;           //!< step size in ms
     1665  double   IntegrationStep_;//!< current integration time step, updated by solver
     1666
     1667  /**
     1668   * Input current injected by CurrentEvent.
     1669   * This variable is used to transport the current applied into the
     1670   * _dynamics function computing the derivative of the state vector.
     1671   * It must be a part of Buffers_, since it is initialized once before
     1672   * the first simulation, but not modified before later Simulate calls.
     1673   */
     1674  double_t I_stim_;
     1675EOF
     1676)
     1677  (pp indent+ "};")
     1678
     1679  (pp indent+
     1680      "template <State_::StateVecElems elem>"
     1681      "double_t get_y_elem_() const { return S_.y_[elem]; }"
     1682      "Parameters_ P_;"
     1683      "State_      S_;"
     1684      "Variables_  V_;"
     1685      "Buffers_    B_;"
     1686      (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
     1687      "}; ")
     1688
     1689  (pp indent+ (,(sprintf #<<EOF
     1690  inline
     1691  port ~A::check_connection(Connection& c, port receptor_type)
     1692  {
     1693    SpikeEvent e;
     1694    e.set_sender(*this);
     1695    c.check_event(e);
     1696    return c.get_target()->connect_sender(e, receptor_type);
     1697  }
     1698
     1699  inline
     1700  port ~A::connect_sender(SpikeEvent&, port receptor_type)
     1701  {
     1702    if (receptor_type != 0)
     1703      throw UnknownReceptorType(receptor_type, get_name());
     1704    return 0;
     1705  }
     1706 
     1707  inline
     1708  port ~A::connect_sender(CurrentEvent&, port receptor_type)
     1709  {
     1710    if (receptor_type != 0)
     1711      throw UnknownReceptorType(receptor_type, get_name());
     1712    return 0;
     1713  }
     1714
     1715  inline
     1716  port ~A::connect_sender(DataLoggingRequest& dlr,
     1717                                      port receptor_type)
     1718  {
     1719    if (receptor_type != 0)
     1720      throw UnknownReceptorType(receptor_type, get_name());
     1721    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
     1722  }
     1723
     1724  inline
     1725    void ~A::get_status(DictionaryDatum &d) const
     1726  {
     1727    P_.get(d);
     1728    S_.get(d);
     1729    Archiving_Node::get_status(d);
     1730
     1731    (*d)[names::recordables] = recordablesMap_.get_list();
     1732
     1733    def<double_t>(d, names::t_spike, get_spiketime_ms());
     1734  }
     1735
     1736  inline
     1737    void ~A::set_status(const DictionaryDatum &d)
     1738  {
     1739    Parameters_ ptmp = P_;  // temporary copy in case of errors
     1740    ptmp.set(d);                       // throws if BadProperty
     1741    State_      stmp = S_;  // temporary copy in case of errors
     1742    stmp.set(d, ptmp);                 // throws if BadProperty
     1743
     1744    // We now know that (ptmp, stmp) are consistent. We do not
     1745    // write them back to (P_, S_) before we are also sure that
     1746    // the properties to be set in the parent class are internally
     1747    // consistent.
     1748    Archiving_Node::set_status(d);
     1749
     1750    // if we get here, temporaries contain consistent set of properties
     1751    P_ = ptmp;
     1752    S_ = stmp;
     1753
     1754    calibrate();
     1755  }
     1756EOF
     1757  sysname sysname sysname
     1758  sysname sysname sysname)))
     1759
     1760  (pp indent "}" ,nl)
     1761  )
     1762
     1763(define (output-update
     1764         sysname method state-index-map steady-state-index-map
     1765         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     1766         reaction-eq-defs i-eqs pool-ions perm-ions
     1767         synapse-info
     1768         indent indent+)
     1769 
     1770  (let* ((vi (lookup-def 'v state-index-map))
     1771         (threshold (sprintf #<<EOF
     1772            // (threshold && maximum)
     1773            if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
     1774              {
     1775                S_.r_ = V_.RefractoryCounts_;
     1776               
     1777                set_spiketime(Time::step(origin.get_steps()+lag+1));
     1778               
     1779                SpikeEvent se;
     1780                network()->send(*this, se, lag);
     1781              }
     1782EOF
     1783        vi vi)))
     1784   
     1785    (pp indent  (,(sprintf #<<EOF
     1786void ~A::update(Time const & origin, const long_t from, const long_t to)
     1787EOF
     1788  sysname)))
     1789
     1790    (pp indent  #\{)
     1791
     1792    (pp indent+ (,(sprintf #<<EOF
     1793assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
     1794    assert(from < to);
     1795
     1796    for ( long_t lag = from ; lag < to ; ++lag )
     1797      {
     1798   
     1799        double tt = 0.0 ; //it's all relative!
     1800        V_.U_old_ = S_.y_[~A];
     1801EOF
     1802vi)))
     1803
     1804    (case method
     1805      ((cvode)
     1806       (pp indent+ (,(sprintf #<<EOF
     1807                             
     1808        // adaptive step integration
     1809        while (tt < B_.step_)
     1810        {
     1811          const int status = CVode(B_.sys, B_.step_, y, &tt, CV_NORMAL);
     1812
     1813          if ( status != CV_SUCCESS )
     1814            throw GSLSolverFailure(get_name(), status);
     1815
     1816        }
     1817
     1818EOF
     1819))))
     1820      ((gsl)
     1821       (pp indent+ (,(sprintf #<<EOF
     1822                             
     1823        // adaptive step integration
     1824        while (tt < B_.step_)
     1825        {
     1826          const int status = gsl_odeiv_evolve_apply(B_.e_, B_.c_, B_.s_,
     1827                                 &B_.sys_,              // system of ODE
     1828                                 &tt,                   // from t...
     1829                                  B_.step_,             // ...to t=t+h
     1830                                 &B_.IntegrationStep_ , // integration window (written on!)
     1831                                  S_.y_);               // neuron state
     1832
     1833          if ( status != GSL_SUCCESS )
     1834            throw GSLSolverFailure(get_name(), status);
     1835        }
     1836EOF
     1837))))
     1838      )
     1839
     1840    (let ((isyns (lookup-def 'i-synapses synapse-info))
     1841          (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
     1842
     1843      (for-each (lambda (isyn psc)
     1844                  (let ((gi (lookup-def (nest-name (fourth isyn)) state-index-map)))
     1845                    (if gi
     1846                        (pp indent+ (,(sprintf "      S_.y_[~A] += B_.spike_~A.get_value(lag);"
     1847                                               gi (second psc)))
     1848                        ))))
     1849                isyns pscs))
     1850
     1851    (pp indent+ (,(sprintf #<<EOF
     1852        // sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
     1853        // refractory?
     1854        if (S_.r_)
     1855          {
     1856            --S_.r_;
     1857          }
     1858        else
     1859          {
     1860           ~A
     1861          }
     1862   
     1863        // set new input current
     1864        B_.I_stim_ = B_.currents_.get_value(lag);
     1865
     1866        // log state data
     1867        B_.logger_.record_data(origin.get_steps() + lag);
     1868
     1869      }
     1870EOF
     1871  (if (lookup-def 'V_t const-defs) threshold "")
     1872  )))
     1873       
     1874  (pp indent  #\})
     1875
     1876))
     1877
     1878(define (output-buffers+nodes
     1879         sysname method state-index-map steady-state-index-map
     1880         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     1881         reaction-eq-defs i-eqs pool-ions perm-ions
     1882         synapse-info
     1883         indent indent+)
     1884
     1885  (let ((N (length state-index-map)))
     1886
     1887  (case method
     1888
     1889
     1890  ((cvode)
     1891
     1892     (pp indent ,nl
     1893         (,(sprintf #<<EOF
     1894~A::Buffers_::Buffers_(~A& n)
     1895    : logger_(n),
     1896      s_(0),
     1897      c_(0),
     1898      e_(0)
     1899{
     1900    // Initialization of the remaining members is deferred to
     1901    // init_buffers_().
     1902}
     1903EOF
     1904  sysname sysname))
     1905
     1906  (,(sprintf #<<EOF
     1907~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
     1908    : logger_(n),
     1909      sys_(0)
     1910{
     1911    // Initialization of the remaining members is deferred to
     1912    // init_buffers_().
     1913}
     1914EOF
     1915  sysname sysname))
     1916    ))
     1917
     1918    ((gsl)
     1919
     1920     (pp indent ,nl
     1921         (,(sprintf #<<EOF
     1922~A::Buffers_::Buffers_(~A& n)
     1923    : logger_(n),
     1924      s_(0),
     1925      c_(0),
     1926      e_(0)
     1927{
     1928    // Initialization of the remaining members is deferred to
     1929    // init_buffers_().
     1930}
     1931EOF
     1932  sysname sysname))
     1933
     1934  (,(sprintf #<<EOF
     1935~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
     1936    : logger_(n),
     1937      s_(0),
     1938      c_(0),
     1939      e_(0)
     1940{
     1941    // Initialization of the remaining members is deferred to
     1942    // init_buffers_().
     1943}
     1944EOF
     1945  sysname sysname))
     1946    ))
     1947
     1948
     1949   )
     1950
     1951
     1952  (pp indent  ,nl  (,(sprintf #<<EOF
     1953~A::~A()
     1954    : Archiving_Node(),
     1955      P_(),
     1956      S_(P_),
     1957      B_(*this)
     1958{
     1959    recordablesMap_.create();
     1960}
     1961EOF
     1962  sysname sysname)))
     1963
     1964
     1965  (pp indent  ,nl  (,(sprintf #<<EOF
     1966~A::~A(const ~A& n)
     1967    : Archiving_Node(n),
     1968      P_(n.P_),
     1969      S_(n.S_),
     1970      B_(n.B_, *this)
     1971{
     1972}
     1973EOF
     1974  sysname sysname sysname)))
     1975
     1976  (pp indent (,(s+ sysname "::~" sysname "()")))
     1977
     1978  (case method
     1979
     1980   ((gsl)
     1981    (pp indent #<<EOF
     1982{
     1983    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
     1984    if ( B_.s_ ) gsl_odeiv_step_free(B_.s_);
     1985    if ( B_.c_ ) gsl_odeiv_control_free(B_.c_);
     1986    if ( B_.e_ ) gsl_odeiv_evolve_free(B_.e_);
     1987}
     1988EOF
     1989  ))
     1990
     1991   ((cvode)
     1992    (pp indent #<<EOF
     1993{
     1994
     1995  if ( B_.sys_ )
     1996  {
     1997    /* Free y vector */
     1998    N_VDestroy_Serial(B.y);
     1999
     2000    /* Free integrator memory */
     2001    CVodeFree(&B.sys_);
     2002  }
     2003}
     2004EOF
     2005    ))
     2006   )
     2007  )
     2008
     2009  (pp indent ,nl (,(sprintf #<<EOF
     2010  void ~A::init_node_(const Node& proto)
     2011{
     2012    const ~A& pr = downcast<~A>(proto);
     2013    P_ = pr.P_;
     2014    S_ = pr.S_;
     2015}
     2016EOF
     2017  sysname sysname sysname)))
     2018
     2019  (pp indent  ,nl (,(sprintf #<<EOF
     2020void ~A::init_state_(const Node& proto)
     2021{
     2022    const ~A& pr = downcast<~A>(proto);
     2023    S_ = pr.S_;
     2024}
     2025EOF
     2026  sysname sysname sysname)))
     2027
     2028  (pp indent ,nl (,(sprintf #<<EOF
     2029void ~A::init_buffers_()
     2030{
     2031EOF
     2032 sysname)))
     2033
     2034     (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
     2035      (for-each
     2036       (lambda (psc)
     2037         (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
     2038       pscs))
     2039    (pp indent+ (#<<EOF
     2040    B_.currents_.clear();           
     2041    Archiving_Node::clear_history();
     2042
     2043    B_.logger_.reset();
     2044
     2045    B_.step_ = Time::get_resolution().get_ms();
     2046    B_.IntegrationStep_ = B_.step_;
     2047
     2048    B_.I_stim_ = 0.0;
     2049EOF
     2050)))
     2051
     2052
     2053    (case method
     2054
     2055    ((cvode)
     2056     
     2057    (pp indent+ (,(sprintf #<<EOF
     2058 
     2059    /* Creates serial vector of length N */
     2060    B.y = N_VNewEmpty_Serial(~A);
     2061    if (check_flag((void *)y, "N_VNew_Serial", 0)) abort();
     2062    N_VSetArrayPointer (S_.y_,B.y);
     2063 
     2064    /* Calls CVodeCreate to create the solver memory and specify the
     2065     * Backward Differentiation Formula and the use of a Newton iteration */
     2066    B.sys_ = CVodeCreate(CV_BDF, CV_NEWTON);
     2067    if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) abort();
     2068
     2069   /* Calls CVodeInit to initialize the integrator memory and specify the
     2070    * right hand side function in y'=f(t,y), the initial time, and
     2071    * the initial values. */
     2072   status = CVodeInit(B.sys_, ~A_dynamics, 0.0, y);
     2073   if (check_flag(&status, "CVodeInit", 1)) abort();
     2074
     2075   /* Calls CVodeSetUserData to configure the integrator to pass the
     2076    * params structure to the right-hand function */
     2077   status = CVodeSetUserData(B.sys_, reinterpret_cast<void*>(this));
     2078   if (check_flag(&status, "CVodeSetUserData", 1)) abort();
     2079EOF
     2080    N sysname))))
     2081
     2082    ((gsl)
     2083
     2084    (pp indent+ (,(sprintf #<<EOF
     2085
     2086    static const gsl_odeiv_step_type* T1 = gsl_odeiv_step_rkf45;
     2087 
     2088    if ( B_.s_ == 0 )
     2089      B_.s_ = gsl_odeiv_step_alloc (T1, ~A);
     2090    else
     2091      gsl_odeiv_step_reset(B_.s_);
     2092   
     2093    if ( B_.c_ == 0 ) 
     2094      B_.c_ = gsl_odeiv_control_y_new (1e-3, 0.0);
     2095    else
     2096      gsl_odeiv_control_init(B_.c_, 1e-3, 0.0, 1.0, 0.0);
     2097   
     2098    if ( B_.e_ == 0 ) 
     2099      B_.e_ = gsl_odeiv_evolve_alloc(~A);
     2100    else
     2101      gsl_odeiv_evolve_reset(B_.e_);
     2102 
     2103    B_.sys_.function  = ~A_dynamics;
     2104    B_.sys_.jacobian  = 0;
     2105    B_.sys_.dimension = ~A;
     2106    B_.sys_.params    = reinterpret_cast<void*>(this);
     2107}
     2108EOF
     2109  N N sysname N)))
     2110  )
     2111
     2112
     2113  (let ((iv (lookup-def 'v state-index-map)))
     2114    (if iv
     2115        (pp indent ,nl (,(sprintf #<<EOF
     2116void ~A::calibrate()
     2117{
     2118    B_.logger_.init(); 
     2119    V_.RefractoryCounts_ = 20;
     2120    V_.U_old_ = S_.y_[~A];
     2121}
     2122EOF
     2123  sysname iv)))
     2124  ))
     2125
     2126)
     2127
     2128
     2129
     2130
Note: See TracChangeset for help on using the changeset viewer.