source: project/release/4/nemo/trunk/nemo-nest.scm @ 29940

Last change on this file since 29940 was 29940, checked in by Ivan Raikov, 8 years ago

nemo: adding option --nest-ss-method to allow the selection of steady state solving method

File size: 103.1 KB
Line 
1;;       
2;; An extension for translating NEMO models to NEST code.
3;;
4;; Copyright 2011-2013 Ivan Raikov and the Okinawa Institute of Science and Technology
5;;
6;; This program is free software: you can redistribute it and/or
7;; modify it under the terms of the GNU General Public License as
8;; published by the Free Software Foundation, either version 3 of the
9;; License, or (at your option) any later version.
10;;
11;; This program is distributed in the hope that it will be useful, but
12;; WITHOUT ANY WARRANTY; without even the implied warranty of
13;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14;; General Public License for more details.
15;;
16;; A full copy of the GPL license can be found at
17;; <http://www.gnu.org/licenses/>.
18;;
19
20(module nemo-nest
21
22        (nemo:nest-translator)
23
24        (import scheme chicken utils data-structures ports extras srfi-1 srfi-13 srfi-69)
25       
26        (require-extension lolevel posix matchable strictly-pretty 
27                           varsubst datatype nemo-core nemo-utils nemo-units
28                           nemo-geometry nemo-defaults nemo-constraints
29                           nemo-gate-complex nemo-synapse nemo-currents
30
31                           )
32
33
34(define C++-ops
35  `(+ - * / > < <= >= =))
36
37
38
39(define (nest-name s)
40  (let ((cs (string->list (->string s))))
41    (let loop ((lst (list)) (cs cs))
42      (if (null? cs) (string->symbol (list->string (reverse lst)))
43          (let* ((c (car cs))
44                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
45                           (else #\_))))
46            (loop (cons c1 lst) (cdr cs)))))))
47
48
49(define builtin-consts
50  (append `(params)
51          (map (lambda (x) (nest-name (s+ "M_" (first x)))) nemo:math-constants)))
52
53
54(define builtin-fns
55  `(+ - * / pow neg abs atan asin acos sin cos exp ln
56      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
57      > < <= >= = and or round ceiling floor max min
58      ))
59
60
61(define (rewrite-pow expr)
62  (match expr
63         (('pow x y)  (if (and (integer? y)  (positive? y))
64                          (if (> y 1)  (let ((tmp (gensym "x")))
65                                         `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
66                              x)
67                          (if (and (number? y) (zero? y)) 1.0 expr)))
68         (else expr)))
69
70           
71(define (rhsexpr/C++ expr)
72  (match expr 
73         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/C++ x)) es)))
74         (('pow x y)  (cond ((and (integer? y) (= y 1)) x)
75                            ((and (number? y) (zero? y)) 1.0)
76                            (else expr)))
77         ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (nest-name s))
78                                             (map (lambda (x) (rhsexpr/C++ x)) es)) expr))
79         (id          (if (symbol? id) 
80                          (if (assoc id nemo:math-constants)
81                              (nest-name (s+ "M_" id))
82                              (nest-name id))
83                          id))
84         ))
85
86
87(define (nest-state-name n s)
88  (nest-name (s+ n s)))
89
90
91(define-syntax pp
92  (syntax-rules ()
93    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
94
95
96(define group/C++   (doc:block 2 (doc:text "(") (doc:text ")")))
97(define block/C++   (doc:block 2 (doc:text "{") (doc:text "}")))
98(define (stmt/C++ x) 
99  (match x
100         (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
101         (else  (doc:cons x (doc:text ";")))))
102
103
104(define (ifthen/C++ c e1 e2)
105  (doc:nest 2
106    (doc:connect (doc:group (doc:connect (doc:text "if") c))
107                 (doc:connect (doc:nest 2 e1)
108                              (doc:nest 2 (doc:connect 
109                                           (doc:text "else") 
110                                           e2))))
111    ))
112
113
114(define (letblk/C++ e1 e2)
115  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
116        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
117        (else (doc:connect (doc:group (doc:nest 2 (stmt/C++ e1)))
118                           (doc:group (doc:nest 2 e2))))))
119       
120
121(define (format-op/C++ indent op args)
122  (let ((op1 (doc:text (->string op))))
123    (if (null? args) op1
124        (match args
125               ((x)      (doc:concat (list op1 x)))
126               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
127               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
128               (lst      (let* ((n   (length lst))
129                                (n/2 (inexact->exact (round (/ n 2)))))
130                           (doc:concat 
131                            (intersperse 
132                             (list (format-op/C++ indent op (take lst n/2 )) op1 
133                                   (format-op/C++ indent op (drop lst n/2 )))
134                             (doc:space)))))))))
135
136
137(define (format-fncall/C++ indent op args)
138  (let ((op1 (doc:text (->string op))))
139    (doc:cons op1 (group/C++ ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
140
141
142(define (name-normalize expr)
143  (match expr 
144         (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
145         (('let bs e)
146          `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
147         ((f . es) 
148          (cons (if (member f builtin-fns) f (nest-name f)) (map name-normalize es)))
149         ((? symbol? ) (nest-name expr))
150         ((? atom? ) expr)))
151
152
153(define (canonicalize-expr/C++ expr)
154  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term)))
155    (let* ((expr1 (if-convert expr))
156           (expr2 (subst-convert expr1 subst-empty))
157           (expr3 (let-lift expr2))
158           (expr4 (name-normalize expr3)))
159      expr4)))
160
161
162(define (format-expr/C++ indent expr . rest) 
163  (let-optionals rest ((rv #f))
164   (let ((indent+ (+ 2 indent)))
165    (match expr
166       (('let bindings body)
167        (letblk/C++
168         (fold-right 
169           (lambda (x ax)
170             (letblk/C++
171              (match (second x)
172                     (('if c t e)
173                      (ifthen/C++
174                       (group/C++ (format-expr/C++ indent c))
175                       (block/C++ (format-expr/C++ indent t (first x)))
176                       (block/C++ (format-expr/C++ indent e (first x)))))
177                     (else
178                      (stmt/C++
179                       (format-op/C++ indent+ " = "
180                                         (list (format-expr/C++ indent (first x) )
181                                               (format-expr/C++ indent (second x)))))))
182              ax))
183           (doc:empty) bindings)
184         (match body
185                (('let _ _) (format-expr/C++ indent body rv))
186                (else
187                 (let ((body1 (doc:nest indent (format-expr/C++ indent body))))
188                   (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) body1)))
189                       body1))))))
190       
191       (('if . rest) (error 'format-expr/C++ "invalid if statement " expr))
192
193       ((op . rest) 
194       (let ((op (case op ((ln) 'log) ((abs) 'fabs) (else op))))
195         (let ((fe
196                (if (member op C++-ops)
197                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
198                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
199                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
200                      (case op
201                        ((/) 
202                         (format-op/C++ indent op 
203                                          (map (lambda (x) 
204                                                 (let ((fx (format-expr/C++ indent+ x)))
205                                                   (if (or (symbol? x) (number? x)) fx
206                                                       (if (or mul? plmin?) (group/C++ fx) fx)))) rest)))
207                        ((*) 
208                         (format-op/C++ indent op 
209                                          (map (lambda (x) 
210                                                 (let ((fx (format-expr/C++ indent+ x)))
211                                                   (if (or (symbol? x) (number? x)) fx
212                                                       (if plmin? (group/C++ fx) fx)))) rest)))
213                       
214                        (else
215                         (format-op/C++ indent op 
216                                          (map (lambda (x) 
217                                                 (let ((fx (format-expr/C++ indent+ x))) fx)) rest)))))
218                   
219                    (let ((op (case op ((neg) '-) (else op))))
220                      (format-fncall/C++ indent op (map (lambda (x) (format-expr/C++ indent+ x)) rest))))))
221           (if rv 
222               (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
223               fe))))
224     
225      (else  (let ((fe (doc:text (->string expr))))
226               (if rv 
227                   (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
228                   fe)))))))
229               
230         
231(define (expr->string/C++ x . rest)
232  (let-optionals rest ((rv #f) (width 72))
233    (sdoc->string (doc:format width (format-expr/C++ 2 x rv)))))
234
235
236(define (make-define-fn sysname )
237  (lambda (indent n proc)
238
239    (let ((lst (procedure-data proc))
240          (indent+ (+ 2 indent)))
241
242      (let ((rt       (or (lookup-def 'rt lst) 'double))
243            (formals  (lookup-def 'formals lst))
244            (vars     (lookup-def 'vars lst))
245            (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
246            (body     (lookup-def 'body lst))
247            (rv       (gensym 'rv)))
248
249        (let ((argument-list 
250               (append
251                (if (null? vars) '() (map (lambda (x) (sprintf "double ~A" (nest-name x))) vars))
252                '("const void* params"))))
253          (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) "{" ))
254          (let* ((body0 (rhsexpr/C++ body))
255                 (body1 (canonicalize-expr/C++ (add-params-to-fncall body0 builtin-fns)))
256                 (lbs   (enum-bnds body1 (list))))
257            (pp indent+ (double ,rv ";"))
258            (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs) ";")))
259            (if (not (null? consts)) 
260                (begin (pp indent+ (double ,(slp ", " (map (compose nest-name car) consts)) ";")
261                           (,(sprintf "const ~A::Parameters_" sysname) "& p =  " 
262                            ,(sprintf " *(reinterpret_cast< const ~A::Parameters_ *>(params));" sysname)))
263                       (for-each (lambda (x) (let ((n (car x))) 
264                                               (pp indent+ (,(nest-name n) = ,(sprintf "p.~A;" (nest-name n))))))
265                                 consts)
266                       ))
267            (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
268            (pp indent+ ,(s+ "return " rv ";"))
269            ))
270        (pp indent "}"))
271    )))
272
273
274(define (make-define-fn-header sysname )
275  (lambda (indent n proc)
276
277    (let ((lst (procedure-data proc))
278          (indent+ (+ 2 indent)))
279
280      (let ((rt       (or (lookup-def 'rt lst) 'double))
281            (formals  (lookup-def 'formals lst))
282            (vars     (lookup-def 'vars lst))
283            (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
284            (body     (lookup-def 'body lst))
285            (rv       (gensym 'rv)))
286
287        (let ((argument-list 
288               (append
289                (if (null? vars) '() (map (lambda (x) (sprintf "double ~A" (nest-name x))) vars))
290                '("const void* params"))))
291          (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) ";" ))
292          ))
293      )))
294
295
296
297(define (ith v i) (sprintf "Ith(~A,~A)" v i))
298
299
300(define (output-prelude sysname indent)
301
302  (pp indent ,#<#EOF
303##include "#{sysname}.h"
304##include "exceptions.h"
305##include "network.h"
306##include "dict.h"
307##include "integerdatum.h"
308##include "doubledatum.h"
309##include "dictutils.h"
310##include "numerics.h"
311##include <limits>
312
313##include "universal_data_logger_impl.h"
314
315##include <iomanip>
316##include <iostream>
317##include <cstdio>
318##include <cstring>
319##include <cmath>
320EOF
321))
322
323
324
325(define kinsol-prelude-header
326#<<EOF
327  /**
328   * Exception to be thrown if a KinSol solver does not return KIN_SUCCESS
329   * @ingroup KernelExceptions
330   */
331  class KINSolverFailure: public KernelException
332  {
333  public:
334    /**
335    * @note model should be passed from get_name() to ensure that
336    *             names of copied models are reported correctly. 
337    * @param model name of model causing problem
338    * @param status exit status of the KINSOL solver
339    */
340    KINSolverFailure(const std::string& model,
341                     const int status)
342      : KernelException("KINSolverFailure"),
343        model_(model),
344        status_(status)
345      {}
346    ~KINSolverFailure() throw() {}
347   
348    std::string message()
349    {
350      std::ostringstream msg;
351      msg << "In model " << model_ << ", the KINSOL solver "
352          << "returned with exit status " << status_ << ".\n";
353      return msg.str();
354    }
355
356    private:
357      const std::string model_;
358      const int status_;
359  };
360
361EOF
362)
363
364
365(define cvode-prelude-header
366#<<EOF
367
368  /**
369   * Exception to be thrown if a CVode solver does not return CV_SUCCESS
370   * @ingroup KernelExceptions
371   */
372  class CVodeSolverFailure: public KernelException
373  {
374  public:
375    /**
376    * @note model should be passed from get_name() to ensure that
377    *             names of copied models are reported correctly. 
378    * @param model name of model causing problem
379    * @param status exit status of the CVode solver
380    */
381    CVodeSolverFailure(const std::string& model,
382                     const int status)
383      : KernelException("CVodeSolverFailure"),
384        model_(model),
385        status_(status)
386      {}
387    ~CVodeSolverFailure() throw() {}
388   
389    std::string message()
390    {
391      std::ostringstream msg;
392      msg << "In model " << model_ << ", the CVODE solver "
393          << "returned with exit status " << status_ << ".\n";
394      return msg.str();
395    }
396
397    private:
398      const std::string model_;
399      const int status_;
400  };
401
402
403EOF
404)
405
406
407(define ida-prelude-header
408#<<EOF
409
410  /**
411   * Exception to be thrown if an IDA solver does not return IDA_SUCCESS
412   * @ingroup KernelExceptions
413   */
414  class IDASolverFailure: public KernelException
415  {
416  public:
417    /**
418    * @note model should be passed from get_name() to ensure that
419    *             names of copied models are reported correctly. 
420    * @param model name of model causing problem
421    * @param status exit status of the IDA solver
422    */
423    IDASolverFailure(const std::string& model,
424                     const int status)
425      : KernelException("IDASolverFailure"),
426        model_(model),
427        status_(status)
428      {}
429    ~IDASolverFailure() throw() {}
430   
431    std::string message()
432    {
433      std::ostringstream msg;
434      msg << "In model " << model_ << ", the IDA solver "
435          << "returned with exit status " << status_ << ".\n";
436      return msg.str();
437    }
438
439    private:
440      const std::string model_;
441      const int status_;
442  };
443
444
445EOF
446)
447
448
449(define sundials-prelude 
450#<<EOF
451
452  /* 
453   * Check function return value.
454   *    opt == 0 means SUNDIALS function allocates memory so check if
455   *             returned NULL pointer
456   *    opt == 1 means SUNDIALS function returns a flag so check if
457   *             flag >= 0
458   *    opt == 2 means function allocates memory so check if returned
459   *             NULL pointer 
460   */
461
462  static int check_flag(void *flagvalue, const char *funcname, int opt)
463  {
464    int *errflag;
465
466    /* Check if CVode/IDA function returned NULL pointer - no memory allocated */
467    if (opt == 0 && flagvalue == NULL) 
468    {
469      fprintf(stderr, 
470              "\nCVode/IDA error: %s() failed - returned NULL pointer\n\n",
471              funcname);
472      return(1);
473    }
474
475    /* Check if flag < 0 */
476    else if (opt == 1) 
477    {
478      errflag = (int *) flagvalue;
479      if (*errflag < 0) {
480        fprintf(stderr, 
481                "\nCVode/IDA error: %s() failed with flag = %d\n\n",
482                funcname, *errflag);
483        return(1); }
484    }
485
486  /* Check if function returned NULL pointer - no memory allocated */
487  else if (opt == 2 && flagvalue == NULL) 
488  {
489    fprintf(stderr, 
490            "\nMemory error: %s() failed - returned NULL pointer\n\n",
491            funcname);
492    return(1);
493  }
494
495    return(0);
496  }
497
498
499  void adjust_zero_crossings (N_Vector v, double abstol)
500  {
501    int i;
502    for (i = 0; i < NV_LENGTH_S(v); i++)
503    {
504        if (fabs(NV_Ith_S(v,i)) < abstol) NV_Ith_S(v,i) = 0;
505    }
506    return;
507  }
508
509EOF
510)
511
512(define (fsolve-prelude ss-method)
513  (case ss-method
514    ((kinsol)
515#<<EOF
516
517  int fsolve (KINSysFn f, int N, N_Vector fval, void *user_data, std::string name)
518  {
519      int status;
520      N_Vector u0, sc;
521      void *kmem;
522
523      u0 = N_VNew_Serial(N);
524      N_VConst_Serial(0.0,u0);
525      N_VConst_Serial(0.0,fval);
526
527      sc = N_VNew_Serial(N);
528      N_VConst_Serial(1.0,sc);
529
530      kmem = KINCreate();
531
532      status = KINSetUserData (kmem, user_data);
533      if (check_flag (&status, "KinSetUserData", 1)) throw KINSolverFailure (name, status);
534
535      status = KINInit (kmem, f, u0);
536      if (check_flag (&status, "KinInit", 1)) throw KINSolverFailure (name, status);
537
538      status = KINDense (kmem, N);
539      if (check_flag (&status, "KinDense", 1)) throw KINSolverFailure (name, status);
540
541      status = KINSol (kmem, fval, KIN_NONE, sc, sc);
542      if (check_flag (&status, "KINSol", 1)) throw KINSolverFailure (name, status);
543     
544      N_VDestroy_Serial(sc);
545      N_VDestroy_Serial(u0);
546
547      KINFree (&kmem);
548
549      return 0;
550  }
551
552EOF
553  )
554  (else
555#<<EOF
556
557  int fsolve (int (*fss)(const gsl_vector *, void *user_data, gsl_vector *), int N, gsl_vector *fval, void *user_data, std::string name)
558  {
559     
560      const gsl_multiroot_fsolver_type * T = gsl_multiroot_fsolver_hybrid;
561      gsl_multiroot_fsolver * s = gsl_multiroot_fsolver_alloc (T, N);
562      gsl_multiroot_function f = {fss, N, user_data};
563
564      int status, iter;
565      gsl_vector *x = gsl_vector_alloc (N);
566     
567      for (int i = 0; i < N; i++)
568      {
569         gsl_vector_set (x, i, 0.0);
570      }
571
572      gsl_multiroot_fsolver_set (s, &f, x);
573
574      iter = 0;
575      do
576      {
577         iter++;
578         status = gsl_multiroot_fsolver_iterate (s);
579
580         if ((status == GSL_EBADFUNC) ||
581             (status == GSL_ENOPROG))
582         {
583            throw GSLSolverFailure(name, status);
584         }
585         
586         status =  gsl_multiroot_test_residual (s->f, 1e-7);
587      }
588      while (status == GSL_CONTINUE && iter < 1000);
589
590      for (int i = 0; i < N; i++)
591      {
592         gsl_vector_set (fval, i, gsl_vector_get (s->x, i));
593      }
594
595      gsl_vector_free (x);
596      gsl_multiroot_fsolver_free (s);
597
598      return 0;
599  }
600
601EOF
602   
603))
604)
605
606
607(define (output-event-handle
608         sysname state-index-map 
609         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
610         reaction-eq-defs i-eqs pool-ions perm-ions 
611         synapse-info
612         indent indent+)
613
614  (let ((isyns  (lookup-def 'i-synapses synapse-info))
615        (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
616 
617  (pp indent  ,nl (,#<#EOF
618void #{sysname}::handle(SpikeEvent & e)
619  {
620    int flag;
621    assert(e.get_delay() > 0);
622    flag = 0;
623EOF
624   ))
625
626  (for-each
627   (lambda (isyn psc)
628     (let ((wscale (fourth isyn))
629           (wthreshold (fifth isyn)))
630
631       (pp indent+ (,#<#EOF
632    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)))
633      {
634        B_.spike_#(nest-name (second psc)).add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
635                                                     fabs(e.get_weight()) * e.get_multiplicity());
636        flag = 1;
637      }
638EOF
639       )) 
640   ))
641   isyns pscs)
642
643  (pp indent ,nl #\} ,nl)
644
645  (pp indent (,#<#EOF
646void #{sysname}::handle(CurrentEvent& e)
647  {
648    assert(e.get_delay() > 0);
649
650    const double_t c = e.get_current();
651    const double_t w = e.get_weight();
652
653    B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), w * c);
654  }
655
656void #{sysname}::handle(DataLoggingRequest& e)
657  {
658    B_.logger_.handle(e);
659  }
660EOF
661  ))
662
663))
664
665
666(define (output-header 
667         sysname method ss-method state-index-map steady-state-index-map 
668         defaults external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
669         reaction-eq-defs i-eqs pool-ions perm-ions 
670         synapse-info
671         indent indent+)
672
673  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
674
675    (pp indent
676        ("#include \"nest.h\"")
677        ("#include \"event.h\"")
678        ("#include \"archiving_node.h\"")
679        ("#include \"ring_buffer.h\"")
680        ("#include \"connection.h\"")
681        ("#include \"universal_data_logger.h\"")
682        ("#include \"recordables_map.h\"")
683        )
684
685    (case method
686      ((leapfrog) 
687       (pp indent
688           ("#define Ith(v,i)    (v[i])   /* Ith component in a vector */")
689           ))
690      ((cvode) 
691       (pp indent
692           ("#include <sundials/sundials_types.h> /* definition of type realtype */")
693           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
694           ("#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */")
695           ("#include <cvode/cvode_diag.h>        /* prototype for CVDiag */")
696           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
697           ))
698      ((ida) 
699       (pp indent
700           ("#include <sundials/sundials_types.h> /* definition of type realtype */")
701           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
702           ("#include <ida/ida.h>                 /* prototypes for IDA fcts., consts. */")
703           ("#include <ida/ida_dense.h>")
704           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
705           ))
706      ((gsl) 
707       (pp indent
708           ("#include <gsl/gsl_errno.h>")
709           ("#include <gsl/gsl_matrix.h>")
710           ("#include <gsl/gsl_sf_exp.h>")
711           ("#include <gsl/gsl_odeiv2.h>")
712           ("#define Ith(v,i)    (v[i])")
713           ))
714      )
715
716    (pp indent ,nl)
717
718    (case ss-method
719      ((kinsol)
720       (pp indent
721           ("#include <sundials/sundials_types.h> /* definition of type realtype */")
722           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
723           ("#include <kinsol/kinsol.h>           /* prototypes for KINSOL fcts. */")
724           ("#include <kinsol/kinsol_dense.h>     /* prototype for KINDense */")))
725      (else
726       (pp indent
727           ("#include <gsl/gsl_vector.h>")
728           ("#include <gsl/gsl_multiroots.h>")))
729           
730      )
731       
732    (pp indent ,nl)
733
734    (pp indent "namespace nest {" ,nl)
735 
736    (case method
737      ((cvode) 
738       (pp indent ,cvode-prelude-header))
739      ((ida) 
740       (pp indent ,ida-prelude-header)))
741
742    (if (not (null? steady-state-index-map))
743        (case method
744          ((cvode ida)
745           (pp indent ,kinsol-prelude-header))))
746     
747 
748    (case method
749      ((leapfrog) 
750       (pp indent 
751           ("typedef int (*rhsfn_t)(double t, const double *y, double *ydot, void *user_data);") ,nl
752           (extern "\"C\""
753                   int ,(sprintf "~A_dynamics" sysname)
754                   "(double, const double *, double *, void*)" #\;) ,nl))
755      ((cvode) 
756       (pp indent (extern "\"C\""
757                          int ,(sprintf "~A_dynamics" sysname)
758                          "(double, const N_Vector, N_Vector, void*)" #\;) ,nl
759                   (extern "\"C\""
760                           int ,(sprintf "~A_event" sysname)
761                           "(double, N_Vector, double *, void*)" #\;) ,nl
762                           ))
763      ((ida) 
764       (pp indent (extern "\"C\""
765                          int ,(sprintf "~A_residual" sysname)
766                          "(double, N_Vector, N_Vector, N_Vector, void*)" #\;) ,nl
767                   (extern "\"C\""
768                           int ,(sprintf "~A_event" sysname)
769                           "(double, N_Vector, N_Vector, double *, void*)" #\;) ,nl
770                           ))
771      (else
772       (pp indent (extern "\"C\""
773                          int ,(sprintf  "~A_dynamics" sysname)
774                          "(double, const double*, double*, void*)" #\;) ,nl))
775      )
776
777    (if (not (null? steady-state-index-map))
778        (case ss-method
779          ((kinsol)
780           (pp indent (extern "\"C\""
781                              int ,(sprintf "~A_steadystate" sysname)
782                              "(N_Vector, N_Vector, void*)" #\;) ,nl)
783           )
784          (else
785           (pp indent (extern "\"C\""
786                              int ,(sprintf "~A_steadystate" sysname)
787                              "(const gsl_vector *, void *, gsl_vector *)" #\;) ,nl)
788           )
789          ))
790
791    (pp indent 
792        (,(sprintf "class ~A:" sysname) "public Archiving_Node { "))
793
794    (pp indent+ ( ,#<#EOF
795public: 
796 ~#{sysname} ();
797 #{sysname} (const #{sysname} &);
798 #{sysname} ();
799EOF
800    ))
801
802    (pp indent ( #<<EOF
803    using Node::connect_sender;
804    using Node::handle;
805
806    port check_connection(Connection&, port);
807   
808    void handle(SpikeEvent &);
809    void handle(CurrentEvent &);
810    void handle(DataLoggingRequest &);
811   
812    port connect_sender(SpikeEvent &, port);
813    port connect_sender(CurrentEvent &, port);
814    port connect_sender(DataLoggingRequest &, port);
815   
816    void get_status(DictionaryDatum &) const;
817    void set_status(const DictionaryDatum &);
818   
819    void init_node_(const Node& proto);
820    void init_state_(const Node& proto);
821    void init_buffers_();
822    void calibrate();
823   
824    void update(Time const &, const long_t, const long_t);
825
826EOF
827))
828
829  (pp indent (,#<#EOF
830
831    /**
832     * Minimal spike receptor type.
833     * @note Start with 1 so we can forbid port 0 to avoid accidental
834     *       creation of connections with no receptor type set.
835     */
836    static const port MIN_SPIKE_RECEPTOR = 1;
837
838    /** 
839     * Spike receptors.
840     */
841    enum SpikeSynapseTypes { #{(slp ", " (append
842                                         (if (pair? pscs)
843                                             (cons (sprintf "~A_SPIKE_RECEPTOR=MIN_SPIKE_RECEPTOR" (first (car pscs)))
844                                                   (map (lambda (x) (sprintf "~A_SPIKE_RECEPTOR" (first x))) (cdr pscs))) '())
845                                         (list "SUP_SPIKE_RECEPTOR ")))}
846                             };
847
848
849EOF
850))
851 
852  (case method
853    ((cvode) 
854     (pp indent+ (friend int ,(sprintf "~A_dynamics" sysname)
855                         "(double, const N_Vector, N_Vector, void*)" #\; ) ,nl
856                 (friend int ,(sprintf "~A_event" sysname)
857                         "(double, N_Vector, double *, void*)" #\;) ,nl
858                         ))
859    ((ida) 
860     (pp indent+ (friend int ,(sprintf "~A_residual" sysname)
861                         "(double, N_Vector, N_Vector, N_Vector, void*)" #\; ) ,nl
862                 (friend int ,(sprintf "~A_event" sysname)
863                         "(double, N_Vector, N_Vector, double *, void*)" #\;) ,nl
864                         ))
865    (else
866     (pp indent+ (friend int ,(sprintf "~A_dynamics" sysname)
867                         "(double, const double*, double*, void*)" #\;) ,nl))
868    )
869
870 
871  (if (not (null? steady-state-index-map))
872      (case ss-method
873        ((kinsol)
874         (pp indent (friend  int ,(sprintf "~A_steadystate" sysname)
875                            "(N_Vector, N_Vector, void*)" #\;) ,nl))
876        (else
877         (pp indent (friend  int ,(sprintf "~A_steadystate" sysname)
878                            "(const gsl_vector *, void *, gsl_vector *)" #\;) ,nl))
879        ))
880         
881
882  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
883    (for-each
884     (lambda (psc)
885       (pp indent+ (,(sprintf "int ~A_transients (long_t lag);" (first psc)))))
886     pscs))
887 
888
889  (pp indent (,#<#EOF
890    // The next two classes need to be friends to access the State_ class/member
891    friend class RecordablesMap<#{sysname}>;
892    friend class UniversalDataLogger<#{sysname}>;
893EOF
894   ))
895
896 
897  (pp indent+ ,nl "struct Parameters_ { ")
898
899  (for-each
900   (lambda (x) (pp indent+ (,(sprintf "double ~A;" x))))
901   (append
902    (map (compose nest-name first) const-defs)
903    (map (compose nest-name first) defaults)))
904
905  (pp indent+ 
906      ("Parameters_();")
907      ("void get(DictionaryDatum&) const;")
908      ("void set(const DictionaryDatum&);"))
909
910  (pp indent+ "};")
911
912  (pp indent+ ,nl "struct State_ { ")
913
914  (pp indent+ ,nl 
915      "enum StateVecElems {"
916      (,(slp ", " (map (lambda (x) 
917                         (let ((n (string-upcase (->string (first x))))
918                               (ni (second x)))
919                           (sprintf "~A = ~A" n ni)))
920                       state-index-map)))
921      "};")
922
923  (pp indent+ ,nl
924      (,(sprintf "double y_[~A];" (length state-index-map)))
925
926      "State_(const Parameters_& p);" 
927      "State_(const State_& s);"
928      "State_& operator=(const State_& s);"
929      "void get(DictionaryDatum&) const;"
930      "void set(const DictionaryDatum&, const Parameters_&);")
931
932  (case method
933    ((gsl leapfrog) 
934     (pp indent+ ,nl "int_t r_; /* refractory counter */")))
935
936  (pp indent+ "};")
937
938  (case method
939    ((cvode ida) 
940     (pp indent+ ,nl "struct Variables_ {};"))
941    (else
942     (pp indent+ ,nl "struct Variables_ { int_t RefractoryCounts_; double U_old_; /* for spike-detection */ };"))
943    )
944
945  (pp indent+ ,nl "struct Buffers_ { ")
946
947  (pp indent+ ,nl 
948      (,(sprintf "Buffers_(~A&);" sysname))
949      (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
950      (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))
951
952  (pp indent+ ,nl
953
954      (,(case method
955
956    ((leapfrog)
957#<<EOF
958  rhsfn_t sys_; // pointer to function that computes dv/dt
959  double *u, *v, *dvdt;  // intermediate state vectors used by leapfrog method
960  unsigned int  N;  // size of state vector used by leapfrog method
961EOF
962)
963
964    ((cvode)
965#<<EOF
966  N_Vector y;    //!< current state vector used by CVode
967  void *   sys_;  //!< CVode control structure
968EOF
969)
970    ((ida)
971#<<EOF
972  N_Vector y, y1;    //!< current state vector used by IDA
973  N_Vector yp;    //!< derivatives vector used by IDA
974  void *   sys_;  //!< IDA control structure
975EOF
976)
977
978    ((gsl)
979#<<EOF
980  gsl_odeiv2_step*    s_;    //!< stepping function
981  gsl_odeiv2_control* c_;    //!< adaptive stepsize control function
982  gsl_odeiv2_evolve*  e_;    //!< evolution function
983  gsl_odeiv2_system   sys_;  //!< struct describing system
984  unsigned int N;  // size of state vector used by Jacobian
985  double *u, *jac;  // intermediate state vectors used for Jacobian approximation
986EOF
987)
988
989    )))
990
991
992  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
993    (for-each
994     (lambda (psc)
995       (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
996     pscs))
997
998  (pp indent+ ,nl
999#<<EOF
1000  RingBuffer currents_;
1001
1002  double_t step_;           //!< step size in ms
1003  double   IntegrationStep_;//!< current integration time step, updated by solver
1004
1005  /** 
1006   * Input current injected by CurrentEvent.
1007   * This variable is used to transport the current applied into the
1008   * _dynamics function computing the derivative of the state vector.
1009   * It must be a part of Buffers_, since it is initialized once before
1010   * the first simulation, but not modified before later Simulate calls.
1011   */
1012  double_t I_stim_;
1013EOF
1014)
1015  (pp indent+ "};")
1016
1017  (pp indent+
1018      "template <State_::StateVecElems elem>"
1019      "double_t get_y_elem_() const { return S_.y_[elem]; }"
1020      "Parameters_ P_;"
1021      "State_      S_;"
1022      "Variables_  V_;"
1023      "Buffers_    B_;"
1024      (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
1025      "}; ")
1026
1027  (pp indent+ (,#<#EOF
1028  inline
1029  port #{sysname}::check_connection(Connection& c, port receptor_type)
1030  {
1031    SpikeEvent e;
1032    e.set_sender(*this);
1033    c.check_event(e);
1034    return c.get_target()->connect_sender(e, receptor_type);
1035  }
1036
1037  inline
1038  port #{sysname}::connect_sender(SpikeEvent&, port receptor_type)
1039  {
1040    if ( receptor_type < MIN_SPIKE_RECEPTOR || receptor_type >= SUP_SPIKE_RECEPTOR )
1041    {
1042      if ( receptor_type < 0 || receptor_type >= SUP_SPIKE_RECEPTOR )
1043        throw UnknownReceptorType(receptor_type, get_name());
1044      else
1045        throw IncompatibleReceptorType(receptor_type, get_name(), "SpikeEvent");
1046    }
1047    return receptor_type;
1048  }
1049 
1050  inline
1051  port #{sysname}::connect_sender(CurrentEvent&, port receptor_type)
1052  {
1053    if (receptor_type != 0)
1054      throw UnknownReceptorType(receptor_type, get_name());
1055    return 0;
1056  }
1057
1058  inline
1059  port #{sysname}::connect_sender(DataLoggingRequest& dlr, 
1060                                      port receptor_type)
1061  {
1062    if (receptor_type != 0)
1063      throw UnknownReceptorType(receptor_type, get_name());
1064    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
1065  }
1066
1067  inline
1068    void #{sysname}::get_status(DictionaryDatum &d) const
1069  {
1070    P_.get(d);
1071    S_.get(d);
1072    Archiving_Node::get_status(d);
1073
1074    (*d)[names::recordables] = recordablesMap_.get_list();
1075
1076    def<double_t>(d, names::t_spike, get_spiketime_ms());
1077
1078    DictionaryDatum receptor_dict_ = new Dictionary();
1079
1080    #{(string-concatenate
1081       (map (lambda (psc) 
1082              (sprintf "    (*receptor_dict_)[Name(\"~A\")]  = ~A_SPIKE_RECEPTOR;~%"
1083                       (first psc) (first psc)))
1084            pscs))}
1085               
1086    (*d)[names::receptor_types] = receptor_dict_;
1087  }
1088
1089  inline
1090    void #{sysname}::set_status(const DictionaryDatum &d)
1091  {
1092    Parameters_ ptmp = P_;  // temporary copy in case of errors
1093    ptmp.set(d);                       // throws if BadProperty
1094    State_      stmp = S_;  // temporary copy in case of errors
1095    stmp.set(d, ptmp);                 // throws if BadProperty
1096
1097    // We now know that (ptmp, stmp) are consistent. We do not 
1098    // write them back to (P_, S_) before we are also sure that 
1099    // the properties to be set in the parent class are internally 
1100    // consistent.
1101    Archiving_Node::set_status(d);
1102
1103    // if we get here, temporaries contain consistent set of properties
1104    P_ = ptmp;
1105    S_ = stmp;
1106
1107    calibrate();
1108  }
1109EOF
1110  ))
1111
1112  (pp indent "}" ,nl)
1113  )
1114)
1115
1116
1117
1118(define (output-jac sysname method state-index-map pool-ions i-eqs v-eq indent indent+)
1119
1120;; Diagonal Jacobian approximation: (f(s+.01) - f(s))/.001
1121
1122;; CVODE:
1123;;  typedef int (*CVDlsDenseJacFn)(long int N, realtype t,
1124;;                               N_Vector y, N_Vector fy,
1125;;                               DlsMat Jac, void *user_data,
1126;;                               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
1127
1128;; GSL:
1129;; int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
1130
1131
1132
1133  (case method
1134
1135    ((gsl)
1136     (pp indent ,nl
1137         (,(sprintf "extern \"C\" int ~A_jacobian" sysname)
1138          (,(slp ", " `("double t"
1139                        "const double y[]"
1140                        "double *dfdy"
1141                        "double dfdt[]"
1142                        "void* pnode"
1143                        ))))
1144          #\{
1145
1146          ("// cast the node ptr to an object of the proper type")
1147          ("assert(pnode);")
1148          ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1149          (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1150
1151          ("// state is a reference to the model state ")
1152          (struct ,(sprintf "~A::Buffers_" sysname) "*b;")
1153          ("b = &(vnode.B_);")
1154          (,#<#EOF
1155   for (int i = 0; i < b->N; i++)
1156   {
1157       b->u[i] = y[i] + 0.01;
1158   }
1159   #{sysname}_dynamics(t, b->u, b->jac, pnode);
1160   for (int i = 0; i < b->N; i++)
1161   {
1162       dfdt[i*b->N + i] = (b->jac[i] - dfdy[i]) / .001;
1163   }
1164   return 0;
1165EOF
1166)
1167   #\}
1168   ))
1169))
1170
1171
1172(define (output-cvode-event sysname method imports const-defs state-index-map 
1173           external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
1174           defaults i-eqs v-eq indent indent+)
1175
1176  (let ((vi (lookup-def 'v state-index-map)))
1177   
1178    (pp indent ,nl
1179        (,(sprintf "extern \"C\" int ~A_event" sysname)
1180         (,(slp ", " `("double t"
1181                       "N_Vector y"
1182                       "double *g"
1183                       "void* pnode"
1184                       )))
1185         #\{
1186         ))
1187
1188      (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;")
1189
1190      ,nl
1191
1192      ("// S is shorthand for the type that describes the model state ")
1193      (typedef ,(sprintf "~A::State_" sysname) "S;")
1194           
1195      ,nl
1196           
1197      ("// cast the node ptr to an object of the proper type")
1198      ("assert(pnode);")
1199      ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1200     
1201      ,nl
1202     
1203      ("// params is a reference to the model parameters ")
1204      (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
1205      ("params = &(node.P_);")
1206     
1207      ,nl
1208     
1209      ("// state is a reference to the model state ")
1210      (const struct ,(sprintf "~A::State_" sysname) "*state;")
1211      ("state = &(node.S_);")
1212     
1213      )
1214     
1215      (let ((vt (lookup-def 'V_t defaults)))
1216        (if (and vi vt)
1217            (pp indent+
1218                (,(expr->string/C++ (ith 'y vi) 'v))
1219                (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt))
1220                (,(expr->string/C++ '(- v vt) "g[0]"))
1221                ,nl
1222                ("return 0; "))
1223            (pp indent+
1224                (,(expr->string/C++ -1.0 "g[0]")))
1225            ))
1226
1227     
1228      (pp indent (#\}))
1229     
1230      ))
1231
1232
1233(define (output-ida-event sysname method imports const-defs state-index-map 
1234           external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
1235           defaults i-eqs v-eq indent indent+)
1236
1237  (let ((vi (lookup-def 'v state-index-map)))
1238   
1239    (pp indent ,nl
1240        (,(sprintf "extern \"C\" int ~A_event" sysname)
1241         (,(slp ", " `("double t"
1242                       "N_Vector y"
1243                       "N_Vector yp"
1244                       "double *g"
1245                       "void* pnode"
1246                       )))
1247         #\{
1248         ))
1249
1250      (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;")
1251
1252      ,nl
1253
1254      ("// S is shorthand for the type that describes the model state ")
1255      (typedef ,(sprintf "~A::State_" sysname) "S;")
1256           
1257      ,nl
1258           
1259      ("// cast the node ptr to an object of the proper type")
1260      ("assert(pnode);")
1261      ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1262     
1263      ,nl
1264     
1265      ("// params is a reference to the model parameters ")
1266      (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
1267      ("params = &(node.P_);")
1268     
1269      ,nl
1270     
1271      ("// state is a reference to the model state ")
1272      (const struct ,(sprintf "~A::State_" sysname) "*state;")
1273      ("state = &(node.S_);")
1274     
1275      )
1276     
1277      (let ((vt (lookup-def 'V_t defaults)))
1278        (if (and vi vt)
1279            (pp indent+
1280                (,(expr->string/C++ (ith 'y vi) 'v))
1281                (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt))
1282                (,(expr->string/C++ '(- v vt) "g[0]"))
1283                ,nl
1284                ("return 0; "))
1285            (pp indent+
1286                (,(expr->string/C++ -1.0 "g[0]")))
1287            ))
1288
1289     
1290      (pp indent (#\}))
1291     
1292      ))
1293
1294
1295
1296
1297(define (cvode-solve spike-handle)
1298#<#EOF
1299        int N = NV_LENGTH_S (B_.y);
1300        tout = (current_steps+lag)*h;
1301
1302        // adaptive step integration
1303        while (tt < tout)
1304        {
1305          const int status = CVode (B_.sys_, tout, B_.y, &tt, CV_NORMAL);
1306
1307          switch (status)
1308            {
1309                case CV_SUCCESS:      continue;
1310                case CV_ROOT_RETURN:  #{spike-handle};
1311                default:              throw CVodeSolverFailure (get_name(), 0);
1312            }
1313        }
1314
1315        for (int i = 0; i < N; i++)
1316        {
1317           S_.y_[i] = Ith(B_.y,i);
1318        }
1319EOF
1320)
1321
1322
1323(define (ida-solve spike-handle)
1324#<#EOF
1325        int N = NV_LENGTH_S (B_.y);
1326        tout = (current_steps+lag)*h;
1327
1328        // adaptive step integration
1329        while (tt < tout)
1330        {
1331          const int status = IDASolve (B_.sys_, tout, &tt, B_.y, B_.yp, IDA_NORMAL);
1332
1333          switch (status)
1334            {
1335                case IDA_SUCCESS:      continue;
1336                case IDA_ROOT_RETURN:  #{spike-handle};
1337                case IDA_TSTOP_RETURN: break;
1338                default:               throw IDASolverFailure (get_name(), 0);
1339            }
1340        }
1341
1342        for (int i = 0; i < N; i++)
1343        {
1344           S_.y_[i] = Ith(B_.y,i);
1345        }
1346EOF
1347)
1348
1349
1350(define (gsl-solve vi spike-handle)
1351#<#EOF
1352
1353   V_.U_old_ = S_.y_[#{vi}];
1354
1355   while (tt < h)
1356   {
1357
1358     const int status = gsl_odeiv2_evolve_apply
1359                         (B_.e_, B_.c_, B_.s_, 
1360                          &B_.sys_,              // system of ODE
1361                          &tt,                   // from t...
1362                          h,             // ...to t=t+h
1363                         &B_.IntegrationStep_ , // integration window (written on!)
1364                          S_.y_);               // neuron state
1365
1366          if ( status != GSL_SUCCESS )
1367            throw GSLSolverFailure(get_name(), status);
1368
1369   }
1370
1371   #{spike-handle}
1372EOF
1373)
1374
1375
1376(define (leapfrog-solve vi spike-handle)
1377#<#EOF
1378        V_.U_old_ = S_.y_[#{vi}];
1379
1380        tt = (current_steps+lag)*h;
1381        tout = (current_steps+lag+1)*h;
1382
1383        while (tt < tout)
1384        {
1385           B_.sys_ (tt, S_.y_, B_.dvdt, reinterpret_cast<void*>(this));
1386
1387           for (int i = 0; i < B_.N; i++)
1388           {
1389              B_.u[i] = S_.y_[i] + 0.5 * h * B_.dvdt[i];
1390           }
1391
1392           B_.sys_ (tt+0.5*h, B_.u, B_.dvdt, reinterpret_cast<void*>(this));
1393
1394           for (int i = 0; i < B_.N; i++)
1395           { 
1396              B_.v[i] = S_.y_[i] + h * B_.dvdt[i];
1397           }
1398
1399           tt += h;
1400        }
1401
1402        for (int i = 0; i < B_.N; i++)
1403        {
1404           S_.y_[i] = B_.v[i];
1405        }
1406 
1407        #{spike-handle}
1408EOF
1409)
1410
1411
1412
1413(define (output-synapse-transients sysname method state-index-map 
1414                                   const-defs transient-event-defs event-external-eq-defs 
1415                                   synapse-info imports constraints indent indent+)
1416
1417  (let (
1418        (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
1419                                    (list (nest-name n) v)))
1420                      (lookup-def 'c-units constraints)))
1421
1422        (isyns  (lookup-def 'i-synapses synapse-info))
1423        (pscs   (lookup-def 'post-synaptic-conductances synapse-info))
1424        (psc-transients (map (lambda (lst) (map nest-name lst)) 
1425                             (lookup-def 'psc-transients synapse-info)))
1426        (getstate (case method
1427                    ((cvode ida) (lambda (i) (sprintf "Ith(B_.y,~A)" i)))
1428                    (else (lambda (i) (sprintf "Ith(S_.y_,~A)" i)))))
1429        )
1430
1431
1432    (if (not (= (length event-external-eq-defs) (length pscs)))
1433        (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" 
1434               event-external-eq-defs pscs))
1435     
1436    (for-each (lambda (psc transients)
1437
1438                (let* (
1439                       (ltransient-event-defs
1440                        (filter (lambda (x) (member (first x) transients))
1441                                transient-event-defs))
1442
1443                       (levent-external-eq-def
1444                        (car
1445                         (fold (lambda (def ax)
1446                                 (let* ((n (nest-name (first def)) )
1447                                        (b (second def))
1448                                        (events (let ((fvs (enum-freevars b '() '())))
1449                                                  (filter (lambda (x) (member (first x) fvs)) event-external-eq-defs)))
1450                                        )
1451                                   (append events ax)))
1452                               '() ltransient-event-defs)))
1453
1454                       (lconsts
1455                        (delete-duplicates
1456                         (fold (lambda (def ax)
1457                                 (let* ((n (nest-name (first def)) )
1458                                        (b (second def))
1459                                        (consts (let ((fvs (enum-freevars b '() '())))
1460                                                  (filter (lambda (x) (member (first x) fvs)) const-defs)))
1461                                       )
1462                                   (append consts ax)))
1463                               '() ltransient-event-defs)
1464                         (lambda (x y) (equal? (first x) (first y))))
1465                        ))
1466                 
1467                  (pp indent (,(sprintf "inline int ~A::~A_transients (long_t lag) {" sysname (first psc))))
1468                 
1469                  (if (not (null? lconsts)) 
1470                      (pp indent+ (double ,(slp ", " (map (compose nest-name car) lconsts)) ";")))
1471
1472                  (let ((n (second levent-external-eq-def)))
1473                    (pp indent+
1474                        (,(sprintf "double_t ~A;" n))
1475                        (,(sprintf "~A = B_.spike_~A.get_value(lag);" 
1476                                   n (nest-name (second psc))))
1477                        )
1478                    )
1479                 
1480                  (pp indent+
1481                      (,(sprintf "if (~A > 0.0)" (second levent-external-eq-def)))
1482                      (#\{))
1483                     
1484                  (let* ((n      (nest-name (first levent-external-eq-def)))
1485                         (nu     (lookup-def n c-units))
1486                         (nscale (and nu (nemo:unit-scale nu)))
1487                         (b      (second levent-external-eq-def))
1488                         )
1489                    (pp (+ 2 indent+)
1490                        (,(sprintf "double_t ~A;" n))
1491                        (,(expr->string/C++ (if nscale `(* ,nscale ,b) b) n)))
1492                    )
1493                 
1494                  (for-each (lambda (def)
1495                              (let* ((n  (nest-name (first def)) )
1496                                     (ni (lookup-def n state-index-map))
1497                                     (b  (second def))
1498                                     (consts (let ((fvs (enum-freevars b '() '())))
1499                                               (filter (lambda (x) (member (first x) fvs)) const-defs)))
1500                                     )
1501                               
1502                              (if (not (null? consts)) 
1503                                  (for-each (lambda (x) 
1504                                              (let ((n (first x)))
1505                                                (pp (+ 2 indent+) (,(nest-name n) = ,(sprintf "P_.~A;" (nest-name n))))))
1506                                            consts)
1507                                  )
1508                             
1509                              (pp (+ 2 indent+ )
1510                                  (,(sprintf "double_t ~A;" n))
1511                                  ,(if ni (expr->string/C++ (getstate ni) n) '())
1512                                  ,(expr->string/C++ b n)
1513                                  ,(if ni (expr->string/C++ n  (getstate ni)) '())
1514                                  )
1515                              ))
1516
1517                            ltransient-event-defs)
1518                 
1519                 
1520                  (pp (+ 2 indent+) "return 1;");
1521                  (pp indent+ (#\}))
1522                  (pp indent+ "return 0;");
1523                  (pp indent (#\}))
1524                  ))
1525
1526               pscs psc-transients)
1527
1528    ))
1529
1530
1531(define (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports indent)
1532
1533  (let (
1534        (isyns  (lookup-def 'i-synapses synapse-info))
1535        (pscs   (lookup-def 'post-synaptic-conductances synapse-info))
1536        )
1537
1538    (pp indent (,(sprintf "int events = 0;")))
1539
1540    (if (not (= (length event-external-eq-defs) (length pscs)))
1541        (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" 
1542               event-external-eq-defs pscs))
1543     
1544    (for-each (lambda (psc)
1545                (pp indent
1546                    (,(sprintf "events = events + ~A_transients (lag);" (first psc))))
1547                )
1548              pscs)
1549
1550    (case method
1551      ((cvode)
1552        (pp indent 
1553            ( "/* Reinitializes CVode state if any synaptic events have occurred */")
1554            ("if (events > 0)")
1555            (#\{)
1556            ("  int status = CVodeReInit (B_.sys_, tt, B_.y);")
1557            ("  if (check_flag(&status, \"CVodeReInit\", 1)) throw CVodeSolverFailure (get_name(), status);")
1558            (#\})
1559            ))
1560      ((ida)
1561        (pp indent 
1562            ( "/* Reinitializes IDA state if any synaptic events have occurred */")
1563            ("if (events > 0)")
1564            (#\{)
1565            ("  int status = IDAReInit (B_.sys_, tt, B_.y, B_.yp);")
1566            ("  if (check_flag(&status, \"IDAReInit\", 1)) throw IDASolverFailure (get_name(), status);")
1567            (#\})
1568            ))
1569      ((gsl)
1570        (pp indent 
1571            ( "/* Resets the GSL stepping function if any synaptic events have occurred */")
1572            ("if (events > 0)")
1573            (#\{)
1574            ("  gsl_odeiv2_step_reset (B_.s_);")
1575            (#\})
1576            ))
1577
1578       )
1579    ))
1580
1581
1582
1583(define (spike-handle method defaults vi abstol)
1584  (case method
1585    ((cvode ida) 
1586#<#EOF
1587              {
1588                set_spiketime(Time::ms(tt));
1589                SpikeEvent se;
1590                for (int i = 0; i < N; i++)
1591                {
1592                  S_.y_[i] = Ith(B_.y,i);
1593                }
1594                network()->send(*this, se, lag);
1595                adjust_zero_crossings(B_.y, #{(or abstol 1e-7)});
1596                continue;
1597              }
1598EOF
1599        )
1600
1601    (else
1602     (if (lookup-def 'V_t defaults)
1603         (sprintf 
1604#<<EOF
1605         if ( S_.r_ > 0 )
1606            S_.r_--;
1607         else
1608            if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
1609              {
1610                S_.r_ = V_.RefractoryCounts_;
1611                set_spiketime(Time::step(origin.get_steps()+lag+1));
1612                SpikeEvent se;
1613                network()->send(*this, se, lag);
1614              }
1615EOF
1616        vi vi)  ))
1617
1618    ))
1619
1620
1621
1622(define (output-update
1623         sysname method state-index-map 
1624         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1625         reaction-eq-defs transient-event-defs
1626         i-eqs pool-ions perm-ions synapse-info 
1627         event-external-eq-defs defaults imports
1628         abstol reltol
1629         indent indent+)
1630
1631  (let ((vi (lookup-def 'v state-index-map)))
1632   
1633    (pp indent  (,(sprintf "void ~A::update(Time const & origin, const long_t from, const long_t to)" sysname)
1634                 #\{))
1635
1636    (pp indent+ 
1637#<<EOF
1638    assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
1639    assert(from < to);
1640
1641    double tout;
1642    long_t current_steps = origin.get_steps();
1643
1644    for ( long_t lag = from ; lag < to ; ++lag )
1645      {
1646        double h = B_.step_;   
1647        double tt = 0.0 ;
1648EOF
1649    )
1650
1651    (case method
1652
1653      ((leapfrog)
1654       (begin
1655         (pp indent+ (,(leapfrog-solve vi (spike-handle method defaults vi abstol))))
1656         ))
1657
1658      ((cvode)
1659       (begin
1660         (pp indent+ (,(cvode-solve (spike-handle method defaults vi abstol) )) ,nl)
1661         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
1662         ))
1663
1664      ((ida)
1665       (begin
1666         (pp indent+ (,(ida-solve (spike-handle method defaults vi abstol) )) ,nl)
1667         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
1668         ))
1669
1670      ((gsl)
1671       (begin
1672         (pp indent+ (,(gsl-solve vi (spike-handle method defaults vi abstol))) ,nl)
1673         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
1674         ))
1675      )
1676
1677
1678    (pp (+ 6 indent+)
1679        ("B_.I_stim_ = B_.currents_.get_value(lag);")
1680        ("B_.logger_.record_data(current_steps + lag);"))
1681
1682  (pp indent+  #\})
1683  (pp indent  #\})
1684  ))
1685
1686
1687(define (output-nodes-leapfrog
1688         sysname method state-index-map 
1689         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1690         reaction-eq-defs i-eqs pool-ions perm-ions 
1691         synapse-info defaults abstol reltol maxstep
1692         indent indent+)
1693
1694  (let ((N (length state-index-map))
1695        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
1696        )
1697
1698    (pp indent ,nl  (,#<#EOF
1699#{sysname}::#{sysname}()
1700    : Archiving_Node(), 
1701      P_(), 
1702      S_(P_),
1703      B_(*this)
1704{
1705    recordablesMap_.create();
1706}
1707EOF
1708  ))
1709
1710  (pp indent  ,nl  (,#<#EOF
1711#{sysname}::#{sysname} (const #{sysname}& n)
1712    : Archiving_Node(n), 
1713      P_(n.P_), 
1714      S_(n.S_),
1715      B_(n.B_, *this)
1716{
1717}
1718EOF
1719  ))
1720
1721
1722  (pp indent #<#EOF
1723#{sysname}::~#{sysname} ()
1724{
1725    if ( B_.u != NULL ) free (B_.u);
1726    if ( B_.v != NULL ) free (B_.v);
1727    if ( B_.dvdt != NULL ) free (B_.dvdt);
1728}
1729EOF
1730  )
1731
1732  (pp indent ,nl (,#<#EOF
1733  void #{sysname}::init_node_(const Node& proto)
1734{
1735    const #{sysname}& pr = downcast<#{sysname}>(proto);
1736    P_ = pr.P_;
1737    S_ = State_(P_);
1738}
1739EOF
1740  ))
1741
1742  (pp indent ,nl (,#<#EOF
1743void #{sysname}::init_state_(const Node& proto)
1744{
1745    const #{sysname}& pr = downcast<#{sysname}>(proto);
1746    S_ = State_(pr.P_);
1747}
1748EOF
1749  ))
1750))
1751
1752
1753(define (output-nodes-cvode
1754         sysname method state-index-map 
1755         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1756         reaction-eq-defs i-eqs pool-ions perm-ions 
1757         synapse-info defaults abstol reltol maxstep
1758         indent indent+)
1759
1760  (let (
1761        (N (length state-index-map))
1762        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
1763        )
1764
1765    (pp indent ,nl  (,#<#EOF
1766#{sysname}::#{sysname}()
1767    : Archiving_Node(), 
1768      P_(), 
1769      S_(P_),
1770      B_(*this)
1771{
1772    recordablesMap_.create();
1773}
1774EOF
1775  ))
1776
1777
1778  (pp indent  ,nl  (,#<#EOF
1779#{sysname}::#{sysname}(const #{sysname}& n)
1780    : Archiving_Node(n), 
1781      P_(n.P_), 
1782      S_(n.S_),
1783      B_(n.B_, *this)
1784{
1785}
1786EOF
1787  ))
1788
1789
1790#|
1791      CVodeGetNumRhsEvals(B_.sys_, &nfevals);
1792      CVodeGetNumSteps(B_.sys_, &nsteps);
1793      CVodeGetNumLinSolvSetups(B_.sys_, &nlinsetups);
1794      CVodeGetNumErrTestFails(B_.sys_, &netfails);
1795      CVodeGetCurrentOrder(B_.sys_, &qcur);
1796      CVodeGetNumStabLimOrderReds(B_.sys_, &nlsred);
1797
1798      printf ("CVODE: %d function evaluations, %d steps, %d linear solves\n",
1799              nfevals, nsteps, nlinsetups);
1800      printf ("CVODE: %d error test failures, %d order reductions\n",
1801              netfails, nlsred);
1802      printf ("CVODE: current order is %d\n",
1803              qcur);
1804|#
1805
1806  (pp indent (,#<#EOF
1807#{sysname}::~#{sysname}()
1808{
1809
1810  if ( B_.sys_ )
1811  {
1812    /* Free y vector */
1813    N_VDestroy_Serial(B_.y);
1814
1815    /* Free integrator memory */
1816    if (B_.sys_ != NULL)
1817    {
1818      CVodeFree(&B_.sys_);
1819      B_.sys_ = NULL;
1820    }
1821  }
1822}
1823EOF
1824   ))
1825
1826
1827  (pp indent ,nl (,#<#EOF
1828void #{sysname}::init_node_(const Node& proto)
1829{
1830    const #{sysname}& pr = downcast<#{sysname}>(proto);
1831    P_ = pr.P_;
1832    S_ = State_(P_);
1833}
1834EOF
1835  ))
1836
1837  (pp indent ,nl (,#<#EOF
1838void #{sysname}::init_state_(const Node& proto)
1839{
1840    const #{sysname}& pr = downcast<#{sysname}>(proto);
1841    S_ = State_(pr.P_);
1842}
1843EOF
1844  ))
1845
1846  (pp indent ,nl (,#<#EOF
1847void #{sysname}::init_buffers_()
1848{
1849EOF
1850  ))
1851
1852  (for-each (lambda (psc)
1853              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
1854            pscs)
1855
1856  (pp indent+ (#<<EOF
1857    B_.currents_.clear();           
1858    Archiving_Node::clear_history();
1859
1860    B_.logger_.reset();
1861
1862    B_.step_ = Time::get_resolution().get_ms();
1863    B_.IntegrationStep_ = B_.step_;
1864
1865    B_.I_stim_ = 0.0;
1866EOF
1867))
1868
1869   (pp indent+ (,#<#EOF
1870
1871    int status, N, rootdir;
1872
1873    N = #{N};
1874    // only positive direction (rising edge) of spike events will be detected
1875    rootdir = 1;
1876
1877    /* Creates serial vector of length N */
1878    B_.y = N_VNew_Serial(N);
1879    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw CVodeSolverFailure (get_name(), 0);
1880
1881    for (int i = 0; i < N; i++)
1882    {
1883       Ith(B_.y,i) = S_.y_[i];
1884    }
1885 
1886    /* Calls CVodeCreate to create the solver memory and specify the 
1887     * Backward Differentiation Formula and the use of a Newton iteration */
1888    B_.sys_ = CVodeCreate(CV_BDF, CV_NEWTON);
1889    if (check_flag((void *)B_.sys_, "CVodeCreate", 0)) throw CVodeSolverFailure (get_name(), 0);
1890
1891   /* Calls CVodeInit to initialize the integrator memory and specify the
1892    * right hand side function in y'=f(t,y), the initial time, and
1893    * the initial values. */
1894    status = CVodeInit (B_.sys_, #{sysname}_dynamics, 0.0, B_.y);
1895    if (check_flag(&status, "CVodeInit", 1)) throw CVodeSolverFailure (get_name(), status);
1896
1897EOF
1898    )
1899
1900    ,(if (lookup-def 'V_t defaults) #<#EOF
1901    /* Spike event handler (detects zero-crossing of V-V_t) */
1902    status = CVodeRootInit(B_.sys_, 1, (CVRootFn)#{sysname}_event);
1903    if (check_flag(&status, "CVodeRootInit", 1)) throw CVodeSolverFailure (get_name(), status);
1904
1905    /* Detect only the rising edge of spikes */
1906    status = CVodeSetRootDirection(B_.sys_, &rootdir);
1907    if (check_flag(&status, "CVodeSetRootDirection", 1)) throw CVodeSolverFailure (get_name(), status);
1908EOF
1909     '())
1910
1911     (,#<#EOF
1912    /* Sets the relative and absolute error tolerances of CVode  */
1913    status = CVodeSStolerances (B_.sys_, #{(or abstol 1e-7)}, #{(or reltol 1e-7)});
1914    if (check_flag(&status, "CVodeSStolerances", 1)) throw CVodeSolverFailure (get_name(), status);
1915
1916    /* Turn on CVode stability limit detection (only applicable for order 3 and greater) */
1917    status = CVodeSetStabLimDet (B_.sys_,true);
1918    if (check_flag(&status, "CVodeSetStabLimDet", 1)) throw CVodeSolverFailure (get_name(), status);
1919
1920    /* Sets the maximum order of CVode  */
1921    status = CVodeSetMaxOrd (B_.sys_,5);
1922    if (check_flag(&status, "CVodeSetMaxOrd", 1)) throw CVodeSolverFailure (get_name(), status);
1923
1924    /* Sets maximum step size. */
1925    status = CVodeSetMaxStep (B_.sys_,#{(or maxstep "B_.step_")});
1926    if (check_flag(&status, "CVodeSetMaxStep", 1)) throw CVodeSolverFailure (get_name(), status);
1927
1928    /* Calls CVodeSetUserData to configure the integrator to pass the 
1929     * params structure to the right-hand function */
1930    status = CVodeSetUserData(B_.sys_, reinterpret_cast<void*>(this));
1931    if (check_flag(&status, "CVodeSetUserData", 1)) throw CVodeSolverFailure (get_name(), status);
1932
1933    /* Initializes diagonal linear solver. */
1934    status = CVDiag (B_.sys_);
1935    if (check_flag(&status, "CVDiag", 1)) throw CVodeSolverFailure (get_name(), status);
1936    }
1937
1938EOF
1939    ))
1940
1941  (pp indent (,#<#EOF
1942void #{sysname}::calibrate()
1943{
1944    B_.logger_.init()
1945}
1946
1947EOF
1948   ))
1949 
1950))
1951
1952
1953(define (output-nodes-ida
1954         sysname method state-index-map 
1955         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1956         reaction-eq-defs i-eqs pool-ions perm-ions 
1957         synapse-info defaults abstol reltol maxstep
1958         indent indent+)
1959
1960  (let (
1961        (N (length state-index-map))
1962        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
1963        )
1964
1965    (pp indent ,nl  (,#<#EOF
1966#{sysname}::#{sysname}()
1967    : Archiving_Node(), 
1968      P_(), 
1969      S_(P_),
1970      B_(*this)
1971{
1972    recordablesMap_.create();
1973}
1974EOF
1975  ))
1976
1977
1978  (pp indent  ,nl  (,#<#EOF
1979#{sysname}::#{sysname}(const #{sysname}& n)
1980    : Archiving_Node(n), 
1981      P_(n.P_), 
1982      S_(n.S_),
1983      B_(n.B_, *this)
1984{
1985}
1986EOF
1987  ))
1988
1989  (pp indent (,#<#EOF
1990#{sysname}::~#{sysname}()
1991{
1992
1993  if ( B_.sys_ )
1994  {
1995    /* Free y vector */
1996    N_VDestroy_Serial(B_.y);
1997    N_VDestroy_Serial(B_.yp);
1998
1999    /* Free integrator memory */
2000    if (B_.sys_ != NULL)
2001    {
2002      IDAFree(&B_.sys_);
2003      B_.sys_ = NULL;
2004    }
2005
2006  }
2007}
2008EOF
2009   ))
2010
2011
2012  (pp indent ,nl (,#<#EOF
2013void #{sysname}::init_node_(const Node& proto)
2014{
2015    const #{sysname}& pr = downcast<#{sysname}>(proto);
2016    P_ = pr.P_;
2017    S_ = State_(P_);
2018}
2019EOF
2020  ))
2021
2022  (pp indent ,nl (,#<#EOF
2023void #{sysname}::init_state_(const Node& proto)
2024{
2025    const #{sysname}& pr = downcast<#{sysname}>(proto);
2026    S_ = State_(pr.P_);
2027}
2028EOF
2029  ))
2030
2031  (pp indent ,nl (,#<#EOF
2032void #{sysname}::init_buffers_()
2033{
2034EOF
2035  ))
2036
2037  (for-each (lambda (psc)
2038              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
2039            pscs)
2040
2041  (pp indent+ (#<<EOF
2042    B_.currents_.clear();           
2043    Archiving_Node::clear_history();
2044
2045    B_.logger_.reset();
2046
2047    B_.step_ = Time::get_resolution().get_ms();
2048    B_.IntegrationStep_ = B_.step_;
2049
2050    B_.I_stim_ = 0.0;
2051EOF
2052))
2053
2054   (pp indent+ (,#<#EOF
2055
2056    int status, N, rootdir;
2057
2058    N = #{N};
2059    // only positive direction (rising edge) of spike events will be detected
2060    rootdir = 1;
2061
2062    /* Creates serial vectors of length N */
2063    B_.y = N_VNew_Serial(N);
2064    B_.y1 = N_VNew_Serial(N);
2065    B_.yp = N_VNew_Serial(N);
2066    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw IDASolverFailure (get_name(), 0);
2067
2068    for (int i = 0; i < N; i++)
2069    {
2070       Ith(B_.y,i) = S_.y_[i];
2071    }
2072 
2073    #{sysname}_dynamics (0.0, B_.y, B_.yp, reinterpret_cast<void*>(this));
2074
2075    /* Calls IDACreate to create the solver memory */ 
2076    B_.sys_ = IDACreate();
2077    if (check_flag((void *)B_.sys_, "IDACreate", 0)) throw IDASolverFailure (get_name(), 0);
2078
2079   /* Calls IDAInit to initialize the integrator memory and specify the
2080    * resdual function, the initial time, and the initial values. */
2081    status = IDAInit (B_.sys_, #{sysname}_residual, 0.0, B_.y, B_.yp);
2082
2083    if (check_flag(&status, "IDAInit", 1)) throw IDASolverFailure (get_name(), status);
2084
2085EOF
2086    )
2087
2088    ,(if (lookup-def 'V_t defaults) #<#EOF
2089    /* Spike event handler (detects zero-crossing of V-V_t) */
2090    status = IDARootInit(B_.sys_, 1, (IDARootFn)#{sysname}_event);
2091    if (check_flag(&status, "IDARootInit", 1)) throw IDASolverFailure (get_name(), status);
2092
2093    /* Detect only the rising edge of spikes */
2094    status = IDASetRootDirection(B_.sys_, &rootdir);
2095    if (check_flag(&status, "IDASetRootDirection", 1)) throw IDASolverFailure (get_name(), status);
2096EOF
2097     '())
2098
2099     (,#<#EOF
2100    /* Sets the relative and absolute error tolerances of IDA  */
2101    status = IDASStolerances (B_.sys_, #{(or abstol 1e-7)}, #{(or reltol 1e-7)});
2102    if (check_flag(&status, "IDASStolerances", 1)) throw IDASolverFailure (get_name(), status);
2103
2104    /* Sets the maximum order of IDA  */
2105    status = IDASetMaxOrd (B_.sys_,5);
2106    if (check_flag(&status, "IDASetMaxOrd", 1)) throw IDASolverFailure (get_name(), status);
2107
2108    /* Sets maximum step size. */
2109    status = IDASetMaxStep (B_.sys_,#{(or maxstep "B_.step_")});
2110    if (check_flag(&status, "IDASetMaxStep", 1)) throw IDASolverFailure (get_name(), status);
2111
2112    /* Calls IDASetUserData to configure the integrator to pass the 
2113     * params structure to the right-hand function */
2114    status = IDASetUserData(B_.sys_, reinterpret_cast<void*>(this));
2115    if (check_flag(&status, "IDASetUserData", 1)) throw IDASolverFailure (get_name(), status);
2116
2117    /* Initializes dense linear solver. */
2118    status = IDADense (B_.sys_, N);
2119    if (check_flag(&status, "IDADense", 1)) throw IDASolverFailure (get_name(), status);
2120
2121    status = IDACalcIC(B_.sys_, IDA_Y_INIT, 0.0);
2122    if (check_flag(&status, "IDACalcIC", 1)) throw IDASolverFailure (get_name(), status);
2123
2124    }
2125
2126EOF
2127    ))
2128
2129  (pp indent (,#<#EOF
2130void #{sysname}::calibrate()
2131{
2132    B_.logger_.init()
2133}
2134
2135EOF
2136   ))
2137 
2138))
2139
2140
2141(define (output-nodes-gsl
2142         sysname method state-index-map 
2143         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2144         reaction-eq-defs i-eqs pool-ions perm-ions 
2145         synapse-info defaults abstol reltol maxstep
2146         indent indent+)
2147
2148  (let ((N    (length state-index-map))
2149        (iv   (lookup-def 'v state-index-map))
2150        (pscs (lookup-def 'post-synaptic-conductances synapse-info)))
2151
2152    (pp indent ,nl  (,#<#EOF
2153#{sysname}::#{sysname}()
2154    : Archiving_Node(), 
2155      P_(), 
2156      S_(P_),
2157      B_(*this)
2158{
2159    recordablesMap_.create();
2160}
2161EOF
2162  ))
2163
2164  (pp indent  ,nl  (,#<#EOF
2165#{sysname}::#{sysname}(const #{sysname}& n)
2166    : Archiving_Node(n), 
2167      P_(n.P_), 
2168      S_(n.S_),
2169      B_(n.B_, *this)
2170{
2171}
2172EOF
2173  ))
2174
2175
2176  (pp indent ,#<#EOF
2177#{sysname}::~#{sysname} ()
2178{
2179    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
2180    if ( B_.s_ != NULL) gsl_odeiv2_step_free (B_.s_);
2181    if ( B_.c_ != NULL) gsl_odeiv2_control_free (B_.c_);
2182    if ( B_.e_ != NULL) gsl_odeiv2_evolve_free (B_.e_);
2183    if ( B_.u != NULL) free (B_.u);
2184    if ( B_.jac != NULL) free (B_.jac);
2185}
2186EOF
2187  )
2188
2189  (pp indent ,nl (,#<#EOF
2190  void #{sysname}::init_node_(const Node& proto)
2191{
2192    const #{sysname}& pr = downcast<#{sysname}>(proto);
2193    P_ = pr.P_;
2194    S_ = State_(P_);
2195}
2196EOF
2197  ))
2198
2199  (pp indent ,nl (,#<#EOF
2200void #{sysname}::init_state_(const Node& proto)
2201{
2202    const #{sysname}& pr = downcast<#{sysname}>(proto);
2203    S_ = State_(pr.P_);
2204}
2205EOF
2206  ))
2207
2208  (pp indent ,nl (,#<#EOF
2209void #{sysname}::init_buffers_()
2210{
2211EOF
2212  ))
2213
2214  (for-each (lambda (psc)
2215              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
2216            pscs)
2217
2218  (pp indent+ (#<<EOF
2219    B_.currents_.clear();           
2220    Archiving_Node::clear_history();
2221
2222    B_.logger_.reset();
2223
2224    B_.step_ = Time::get_resolution().get_ms();
2225    B_.IntegrationStep_ = B_.step_;
2226
2227    B_.I_stim_ = 0.0;
2228EOF
2229))
2230
2231    (pp indent+ (,#<#EOF
2232    static const gsl_odeiv2_step_type* T1 = gsl_odeiv2_step_rk2;
2233    B_.N = #{N};
2234 
2235    if ( B_.s_ == 0 )
2236      B_.s_ = gsl_odeiv2_step_alloc (T1, B_.N);
2237    else
2238      gsl_odeiv2_step_reset(B_.s_);
2239   
2240    if ( B_.c_ == 0 ) 
2241      B_.c_ = gsl_odeiv2_control_standard_new (#{(or abstol 1e-7)}, #{(or reltol 1e-7)}, 1.0, 0.0);
2242    else
2243      gsl_odeiv2_control_init(B_.c_, #{(or abstol 1e-7)}, #{(or reltol 1e-7)}, 1.0, 0.0);
2244   
2245    if ( B_.e_ == 0 ) 
2246      B_.e_ = gsl_odeiv2_evolve_alloc(B_.N);
2247    else
2248      gsl_odeiv2_evolve_reset(B_.e_);
2249 
2250    B_.sys_.function  = #{sysname}_dynamics;
2251    B_.sys_.jacobian  = #{sysname}_jacobian;
2252    B_.sys_.dimension = B_.N;
2253    B_.sys_.params    = reinterpret_cast<void*>(this);
2254
2255    B_.u = (double *)malloc(sizeof(double) * B_.N);
2256    assert (B_.u);
2257    B_.jac = (double *)malloc(sizeof(double) * B_.N);
2258    assert (B_.jac);
2259
2260EOF
2261   ))
2262
2263  (pp indent ,nl #\})
2264
2265
2266  (pp indent (,#<#EOF
2267void #{sysname}::calibrate()
2268{
2269    B_.logger_.init()
2270    #{(if iv (sprintf "V_.U_old_ = S_.y_[~A];" iv) "")}
2271    #{(if (lookup-def 't_ref defaults) 
2272          "V_.RefractoryCounts_ = Time(Time::ms(P_.t_ref)).get_steps();" 
2273          "V_.RefractoryCounts_ = 0;")}
2274}
2275EOF
2276))
2277))
2278
2279
2280
2281(define (output-buffers
2282         sysname method state-index-map 
2283         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2284         reaction-eq-defs i-eqs pool-ions perm-ions 
2285         synapse-info defaults
2286         indent indent+)
2287
2288  (let ((N (length state-index-map)))
2289
2290    (case method
2291
2292      ((leapfrog)
2293
2294       (pp indent ,nl
2295           (,#<#EOF
2296#{sysname}::Buffers_::Buffers_(#{sysname}& n)
2297    : logger_(n),
2298      sys_(0),
2299      N(0),
2300      u(0),
2301      v(0),
2302      dvdt(0)
2303{
2304    // Initialization of the remaining members is deferred to
2305    // init_buffers_().
2306}
2307EOF
2308)))
2309
2310   ((cvode ida)
2311
2312    (pp indent ,nl
2313        (,#<#EOF
2314#{sysname}::Buffers_::Buffers_(#{sysname}& n)
2315    : logger_(n),
2316      sys_(0)
2317{
2318    // Initialization of the remaining members is deferred to
2319    // init_buffers_().
2320}
2321EOF
2322        )
2323       (,#<#EOF
2324#{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n)
2325    : logger_(n),
2326      sys_(0)
2327{
2328    // Initialization of the remaining members is deferred to
2329    // init_buffers_().
2330}
2331EOF
2332    )))
2333
2334    ((gsl)
2335
2336     (pp indent ,nl
2337         ( ,#<#EOF
2338#{sysname}::Buffers_::Buffers_(#{sysname}& n)
2339    : logger_(n),
2340      s_(0),
2341      c_(0),
2342      e_(0),
2343      N(0),
2344      u(0),
2345      jac(0)
2346{
2347    // Initialization of the remaining members is deferred to
2348    // init_buffers_().
2349}
2350
2351EOF
2352          )
2353        ( ,#<#EOF
2354#{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n)
2355    : logger_(n),
2356      s_(0),
2357      c_(0),
2358      e_(0),
2359      N(0),
2360      u(0),
2361      jac(0)
2362{
2363    // Initialization of the remaining members is deferred to
2364    // init_buffers_().
2365}
2366EOF
2367    )))
2368   ))
2369)
2370
2371
2372(define (output-accessors+modifiers
2373         sysname imports defaults state-index-map const-defs asgn-eq-defs rate-eq-defs 
2374         reaction-eq-defs i-eqs pool-ions perm-ions constraints
2375         indent indent+)
2376
2377  (let ((c-eqs (lookup-def 'c-eqs constraints))
2378       
2379        (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
2380                                     (list (nest-name n) v)))
2381                      (lookup-def 'c-units constraints)))
2382        )
2383
2384    (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
2385    (pp indent  #\{)
2386   
2387    (for-each
2388     (lambda (def)
2389       (let ((n (first def)))
2390         (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
2391     (append const-defs
2392             defaults))
2393   
2394    (pp indent  #\})
2395   
2396    (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
2397    (pp indent  #\{)
2398   
2399    (for-each
2400     (lambda (def)
2401       (let* ((n (first def))
2402              (nu (lookup-def n c-units))
2403              (scale (and nu (nemo:unit-scale nu)))
2404              )
2405         (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))
2406         (if scale (pp indent+ (,(sprintf "~A = ~A * ~A;" n scale n ))))
2407         ))
2408     (append const-defs
2409             defaults))
2410
2411    (for-each
2412     (lambda (def)
2413       (let ((n (first def))
2414             (b (second def)))
2415         (pp indent+ (,(expr->string/C++ (nest-name n) (nest-name b))))))
2416     defaults)
2417   
2418    (pp indent  #\})
2419   
2420    (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
2421    (pp indent  #\{)
2422   
2423    (for-each
2424     (lambda (def)
2425       (let* (
2426              (n (first def))
2427              (i (second def))
2428              )
2429         (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))
2430         ))
2431     state-index-map)
2432   
2433    (let ((vi (lookup-def 'v state-index-map)))
2434      (if vi
2435          (pp indent+ (,(sprintf "def<double_t>(d, names::V_m, y_[~A]);"  vi) ))
2436          ))
2437   
2438    (pp indent  #\})
2439   
2440    (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
2441    (pp indent  #\{)
2442   
2443    (for-each
2444     (lambda (def)
2445       (let* ((n     (first def)) 
2446              (i     (second def))
2447              (nu    (lookup-def n c-units))
2448              (scale (and nu (nemo:unit-scale nu)))
2449              )
2450         (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))
2451         (if scale (pp indent+ (,(sprintf "y_[~A] = ~A * y_[~A];" i scale i ))))
2452         ))
2453     state-index-map)
2454   
2455    (let ((vi (lookup-def 'v state-index-map)))
2456      (if vi
2457          (pp indent+ (,(sprintf "updateValue<double_t>(d, names::V_m, y_[~A]);"  vi) ))
2458          ))
2459   
2460    (pp indent  #\})
2461   
2462    ))
2463
2464
2465(define (output-init sysname method ss-method state-index-map steady-state-index-map 
2466                     imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2467                     reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints 
2468                     indent indent+)
2469
2470 
2471  (let* ((N (length state-index-map))
2472
2473         (c-eqs (lookup-def 'c-eqs constraints))
2474
2475         (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
2476                                     (list (nest-name n) v)))
2477                       (lookup-def 'c-units constraints)))
2478
2479         (i-eqs 
2480          (map
2481           (lambda (def) (list (first def) 
2482                               (add-params-to-fncall 
2483                                (canonicalize-expr/C++ (second def)) builtin-fns)))
2484           i-eqs))
2485
2486         (init-eqs 
2487          (append
2488             
2489           (map (lambda (def)
2490                  (let ((n (first def))
2491                        (b (second def)))
2492                    (list (nest-name n) (nest-name b))))
2493                external-eq-defs)
2494           
2495           asgn-eq-defs
2496           init-eq-defs
2497
2498           (map (lambda (pool-ion)
2499                  (let ((n (pool-ion-in pool-ion))
2500                        (b (pool-ion-inq pool-ion)))
2501                    (list n b)))
2502                pool-ions)
2503
2504           (map (lambda (pool-ion)
2505                  (let ((n (pool-ion-out pool-ion))
2506                        (b (pool-ion-outq pool-ion)))
2507                    (list n b)))
2508                pool-ions)
2509           ))
2510
2511         (init-dag 
2512          (map (lambda (def)
2513                 (cons (first def) (enum-freevars (second def) '() '())))
2514               init-eqs))
2515
2516         (init-order
2517          (reverse
2518           (topological-sort init-dag 
2519                             (lambda (x y) (string=? (->string x) (->string y))))))
2520
2521         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
2522
2523         )
2524   
2525    (pp indent ,nl
2526        (,(sprintf "~A::State_::State_" sysname) (const Parameters_& p))
2527        ( #\{))
2528
2529    (pp indent+
2530        (double ,(slp ", " (delete-duplicates
2531                            (map ->string 
2532                                 (filter (lambda (x) (not (member x builtin-consts)))
2533                                         (append
2534                                          init-locals
2535                                          init-order
2536                                          (map first external-eq-defs)
2537                                          (map pool-ion-in pool-ions)
2538                                          (map pool-ion-out pool-ions)
2539                                          (map first i-eqs)
2540                                          (map first state-index-map) 
2541                                          (map first const-defs)
2542                                          (map first reaction-eq-defs)
2543                                          )))
2544                                 string=?))
2545                      ";")
2546        ("const Parameters_ *params = &p;")
2547        )
2548
2549    (case method
2550      ((gsl leapfrog) 
2551       (pp indent+ "r_ = 0;" ,nl)))
2552
2553    (pp indent+ ,(sprintf "memset(y_,0,~A*sizeof(double));" N))
2554
2555    (for-each
2556     (lambda (x) 
2557       (let* ((n  (first x)))
2558         (pp indent+ ,(expr->string/C++ (sprintf "p.~A" n) n))))
2559     const-defs)
2560
2561    (let ((vi (lookup-def 'v state-index-map))
2562          (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
2563      (if (and vi vrest) 
2564          (pp indent+ ,(expr->string/C++  vrest 'v ))))
2565
2566    (for-each (lambda (n)
2567                (let ((b  (lookup-def n init-eqs)))
2568                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
2569              init-order)
2570   
2571    (for-each (lambda (def)
2572                (let* ((n  (first def)) 
2573                       (ni (lookup-def n state-index-map)))
2574                  (if ni (pp indent+ ,(expr->string/C++ n  (sprintf "y_[~A]" ni))))))
2575              init-eq-defs)
2576
2577    (for-each
2578     (lambda (eq) 
2579       (match-let (((op left right)  eq))
2580                  (pp indent+ ,(sprintf "if (!((~A) ~A (~A))) { throw BadProperty (~S); }; " 
2581                                        (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ left))) op
2582                                        (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ right)))
2583                                        (sprintf "Constraint ~A is not satisfied." eq)
2584                                        ))
2585                  ))
2586     c-eqs)
2587
2588
2589    (if (not (null? steady-state-index-map))
2590
2591        (case ss-method
2592
2593          ((kinsol)
2594
2595           (let ((N (length steady-state-index-map))
2596                 (ssvect (gensym 'ssvect)))
2597
2598             (pp indent+
2599                 ,(sprintf "N_Vector ~A;" ssvect)
2600                 ,(expr->string/C++ (sprintf "N_VNew_Serial(~A)" N) ssvect)
2601                 (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ,N ,ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
2602             
2603             (for-each
2604              (lambda (def)
2605                (let* ((n (first def)) 
2606                       (si (lookup-def n steady-state-index-map))
2607                       (ni (lookup-def n state-index-map)))
2608                  (if si
2609                        (pp indent+
2610                            ,(expr->string/C++ (ith ssvect si)
2611                                               (sprintf "y_[~A]" ni))
2612                            ,(expr->string/C++ (ith ssvect si)
2613                                               n))
2614                        )
2615                  ))
2616              rate-eq-defs)
2617             
2618             (pp indent+ (,(sprintf "N_VDestroy_Serial (~A);" ssvect)))
2619
2620             ))
2621
2622          (else
2623
2624           (let ((N (length steady-state-index-map))
2625                 (ssvect (gensym 'ssvect))
2626                 )
2627
2628             (pp indent+
2629                 ,(sprintf "gsl_vector *~A;" ssvect)
2630                 ,(expr->string/C++ (sprintf "gsl_vector_alloc (~A)" N) ssvect)
2631                 (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ,N ,ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
2632             
2633             (for-each
2634              (lambda (def)
2635                (let* ((n (first def)) 
2636                       (si (lookup-def n steady-state-index-map))
2637                       (ni (lookup-def n state-index-map)))
2638                  (if si
2639                        (pp indent+
2640                            ,(expr->string/C++ (sprintf "gsl_vector_get (~A,~A)" ssvect si)
2641                                               (sprintf "y_[~A]" ni))
2642                            ,(expr->string/C++ (sprintf "gsl_vector_get (~A,~A)" ssvect si)
2643                                               n))
2644                        )
2645                  ))
2646              rate-eq-defs)
2647             
2648             (pp indent+ (,(sprintf "gsl_vector_free (~A);" ssvect)))
2649
2650             ))
2651
2652          ))
2653   
2654    (for-each
2655     (lambda (def)
2656       (let ((n (first def)) (b (second def)))
2657         (if (not (lookup-def n init-eq-defs))
2658             (pp indent+ ,(expr->string/C++ b n)))))
2659     reaction-eq-defs)
2660   
2661    (for-each
2662     (lambda (def) 
2663       (pp indent+ ,(expr->string/C++ (second def) (first def))))
2664     i-eqs)
2665
2666    (let ((vi (lookup-def 'v state-index-map)))
2667      (if vi (pp indent+ ,(expr->string/C++ 'v (sprintf "y_[~A]" vi)))))
2668     
2669    (pp indent  #\})
2670
2671    (pp indent ,nl
2672        (,(sprintf "~A::State_::State_" sysname) (const State_& s)  )
2673        ( #\{)
2674       
2675        ,(case method
2676          ((gsl leapfrog) 
2677           `("r_ = s.r_;"))
2678          (else ""))
2679       
2680        (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N))
2681        ( #\})
2682        ,nl
2683        )
2684
2685    (pp indent ,nl
2686        (,(sprintf "~A::State_& ~A::State_::operator=(const State_& s)" sysname sysname))
2687        ( #\{)
2688        )
2689
2690    (case method
2691      ((gsl leapfrog) 
2692       (pp indent+ "r_ = s.r_;" ,nl)))
2693
2694    (pp indent+ ( ,#<#EOF
2695     assert(this != &s)
2696     for ( size_t i = 0 ; i < #{N} ; ++i )
2697       y_[i] = s.y_[i];
2698EOF
2699))
2700
2701        (pp indent+ "return *this;")
2702
2703        (pp indent  #\})
2704))
2705
2706
2707(define (output-parameters sysname imports const-defs external-eq-defs defaults indent indent+)
2708
2709  (let* ((parameter-defs
2710          (map
2711           (lambda (def) (list (first def) 
2712                               (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
2713           const-defs))
2714
2715         (parameter-locals  (find-locals (map second parameter-defs)))
2716
2717         (const-exprs
2718            (map (lambda (def)
2719                   (let* ((n  (first def)) (b (second def)))
2720                     (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
2721                 const-defs) )
2722             
2723          (default-exprs 
2724            (map (lambda (def)
2725                   (let ((n (first def))
2726                         (b (second def)))
2727                     (expr->string/C++ (nest-name b) (nest-name n))))
2728                 defaults))
2729           )
2730
2731      (pp indent ,nl (,(sprintf "~A::Parameters_::Parameters_" sysname) () ":"))
2732     
2733      (if (not (null? parameter-locals)) 
2734          (pp indent+ (double ,(slp ", " parameter-locals) ";")))
2735
2736      (pp indent+ ,(slp ",\n" const-exprs))
2737     
2738      (pp indent  #\{ ,(slp "\n" default-exprs) #\})
2739     
2740      ))
2741
2742
2743(define (output-recordables
2744         sysname state-index-map 
2745         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2746         reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
2747 
2748  (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
2749  (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
2750  (pp indent #\{)
2751
2752  (for-each
2753   (lambda (def)
2754     (let ((n (first def)) (i (second def)))
2755       (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);" 
2756                              (->string n) sysname sysname (string-upcase (->string n)))))
2757       ))
2758   state-index-map)
2759
2760  (if (lookup-def 'v state-index-map)
2761      (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);" 
2762                             sysname sysname (string-upcase (->string 'v))))))
2763
2764  (pp indent #\})
2765  )
2766
2767
2768
2769(define (output-dy sysname method imports const-defs state-index-map 
2770                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
2771                   pool-ions i-eqs v-eq constraints indent indent+)
2772
2773  (let* (
2774         (c-units (map (lambda (x) 
2775                         (let ((n (first x)) (v (second x)))
2776                           (list (nest-name n) v)))
2777                       (lookup-def 'c-units constraints)))
2778
2779         (i-eqs 
2780          (map
2781           (lambda (def) (list (first def) 
2782                               (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
2783           i-eqs))
2784
2785         (v-eq
2786          (and v-eq
2787               (list (first v-eq) 
2788                     (add-params-to-fncall (canonicalize-expr/C++ (second v-eq)) builtin-fns))))
2789
2790         (eqs 
2791          (append
2792
2793           external-eq-defs
2794           asgn-eq-defs
2795           reaction-eq-defs
2796                 
2797           (map (lambda (pool-ion)
2798                  (let ((n (pool-ion-in pool-ion))
2799                        (b (pool-ion-inq pool-ion)))
2800                    (list n b)))
2801                pool-ions)
2802
2803           (map (lambda (pool-ion)
2804                  (let ((n (pool-ion-out pool-ion))
2805                        (b (pool-ion-outq pool-ion)))
2806                    (list n b)))
2807                pool-ions)
2808
2809           i-eqs
2810          ))
2811         
2812         (eq-dag 
2813          (map (lambda (def)
2814                 (cons (first def) (enum-freevars (second def) '() '())))
2815               eqs))
2816
2817         (eq-order
2818          (reverse
2819           (topological-sort eq-dag 
2820                             (lambda (x y) (string=? (->string x) (->string y))))))
2821
2822         (eq-locals  (find-locals 
2823                      (map second
2824                           (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
2825                               (append i-eqs rate-eq-defs eqs)))))
2826         )
2827
2828
2829    (case method
2830      ((leapfrog)
2831       (pp indent ,nl (
2832                       ,(sprintf "extern \"C\" int ~A_dynamics" sysname)
2833                       (,(slp ", " `("double t"
2834                                     "const double y[]"
2835                                     "double f[]"
2836                                     "void* pnode"
2837                                     )))
2838                       #\{
2839                       )))
2840      ((cvode ida)
2841       (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname)
2842                       (,(slp ", " `("double t"
2843                                     "N_Vector y"
2844                                     "N_Vector f"
2845                                     "void* pnode"
2846                                     )))
2847                       #\{
2848                       )))
2849      ((gsl)
2850       (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname)
2851                       (,(slp ", " `("double t"
2852                                     "const double y[]"
2853                                     "double f[]"
2854                                     "void* pnode"
2855                                     )))
2856                       #\{
2857                       )))
2858      )
2859
2860
2861
2862    (pp indent+ (double ,(slp ", " (delete-duplicates 
2863                                    (map (compose ->string nest-name)
2864                                         (filter (lambda (x) 
2865                                                   (not (member x builtin-consts)))
2866                                                 (append
2867                                                  eq-locals
2868                                                  eq-order
2869                                                  (map first i-eqs)
2870                                                  (map first external-eq-defs)
2871                                                  (map first state-index-map)
2872                                                  (map first const-defs)
2873                                                  )))
2874                                    string=?))
2875                        ";"))
2876
2877    (pp indent+ ,nl
2878
2879        ("// S is shorthand for the type that describes the model state ")
2880        (typedef ,(sprintf "~A::State_" sysname) "S;")
2881       
2882        ,nl
2883       
2884        ("// cast the node ptr to an object of the proper type")
2885        ("assert(pnode);")
2886        ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2887        (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2888       
2889        ,nl
2890
2891        ("// params is a reference to the model parameters ")
2892        (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
2893        ("params = &(node.P_);")
2894       
2895        ,nl
2896
2897        ("// state is a reference to the model state ")
2898        (struct ,(sprintf  "~A::State_" sysname) "*state;")
2899        ("state = &(vnode.S_);")
2900       
2901        ,nl
2902
2903        ("// y[] must be the state vector supplied by the integrator, ")
2904        ("// not the state vector in the node, node.S_.y[]. ")
2905
2906        ,nl
2907
2908        )
2909
2910    (for-each (lambda (def)
2911                (let ((n (first def)) )
2912                  (pp indent+
2913                      ,(expr->string/C++ (sprintf "params->~A" n) n))))
2914              const-defs)
2915
2916    (let ((vi (lookup-def 'v state-index-map)))
2917      (if vi (pp indent+ ,(expr->string/C++ (ith 'y vi) 'v))))
2918
2919    (for-each (lambda (def)
2920                (let* ((n (first def))
2921                       (ni (lookup-def n state-index-map)))
2922                  (pp indent+ ,(expr->string/C++ (ith 'y ni) (nest-name n)))))
2923              rate-eq-defs)
2924
2925    (for-each (lambda (n)
2926                (let ((b  (lookup-def n eqs)))
2927                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
2928              eq-order)
2929
2930    (for-each (lambda (def)
2931                (let* ((n (first def))
2932                       (b (second def))
2933                       (fv (ith 'f (lookup-def n state-index-map)))
2934                       )
2935                  (pp indent+ ,(s+ "// rate equation " n)
2936                      ,(expr->string/C++ b fv))))
2937              rate-eq-defs)
2938   
2939    (if v-eq
2940        (let ((vi (lookup-def 'v state-index-map)))
2941          (if vi
2942              (pp indent+ ,(expr->string/C++ (second v-eq) (ith 'f vi)))
2943              )))
2944
2945    (case method
2946      ((leapfrog)
2947        (pp indent+ ,nl ("return 0;")))
2948      ((cvode ida)
2949        (pp indent+ ,nl ("return 0;")))
2950      ((gsl)
2951        (pp indent+ ,nl ("return GSL_SUCCESS;")))
2952      )
2953
2954    (pp indent  #\})
2955
2956    ))
2957
2958
2959(define (output-residual sysname method imports const-defs state-index-map 
2960                         external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
2961                         pool-ions i-eqs v-eq constraints indent indent+)
2962  (case method 
2963      ((ida)
2964       (begin
2965         (pp indent ,nl (,(sprintf "extern \"C\" int ~A_residual" sysname)
2966                         (,(slp ", " `("double t"
2967                                       "N_Vector y"
2968                                       "N_Vector yp"
2969                                       "N_Vector f"
2970                                       "void* pnode"
2971                                       )))
2972                         #\{
2973                         ))
2974
2975         (pp indent+ ,nl
2976             (int status #\; )
2977
2978             ,nl
2979
2980             ("// cast the node ptr to an object of the proper type")
2981             ("assert(pnode);")
2982             ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2983             (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2984             ("N_Vector y1 = vnode.B_.y1;")
2985             
2986             ,nl
2987
2988             (,(sprintf "status = ~A_dynamics (t, y, y1, pnode);" sysname))
2989
2990             ,nl)
2991             
2992         (for-each (lambda (def)
2993                     (let* ((n (first def))
2994                            (i (lookup-def n state-index-map))
2995                            (fv (ith 'f i))
2996                            (y1v (ith 'y1 i))
2997                            (ypv (ith 'yp i))
2998                            )
2999                       (pp indent+
3000                           ,(expr->string/C++ `(- ,y1v ,ypv) fv))))
3001                   rate-eq-defs)
3002
3003         (if v-eq
3004             (let ((vi (lookup-def 'v state-index-map)))
3005               (if vi
3006                   (let ((v1 (ith 'y1 vi))
3007                         (vp (ith 'yp vi)))
3008                     (pp indent+ ,(expr->string/C++ `(- ,v1 ,vp) (ith 'f vi))))
3009                   )))
3010
3011
3012         (pp indent
3013             ("return status; ")
3014             #\})
3015         ))
3016      ))
3017
3018
3019(define (output-steadystate sysname method ss-method state-index-map steady-state-index-map 
3020                            imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3021                            reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints 
3022                            indent indent+)
3023
3024 
3025  (let* (
3026         (N (length steady-state-index-map))
3027
3028         (c-eqs (lookup-def 'c-eqs constraints))
3029
3030         (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
3031                                     (list (nest-name n) v)))
3032                       (lookup-def 'c-units constraints)))
3033
3034         (i-eqs 
3035          (map
3036           (lambda (def) (list (first def) 
3037                               (add-params-to-fncall 
3038                                (canonicalize-expr/C++ (second def)) builtin-fns)))
3039           i-eqs))
3040
3041         (init-eqs 
3042          (append
3043             
3044           (map (lambda (def)
3045                  (let ((n (first def))
3046                        (b (second def)))
3047                    (list (nest-name n) (nest-name b))))
3048                external-eq-defs)
3049           
3050           asgn-eq-defs
3051           init-eq-defs
3052
3053           (map (lambda (pool-ion)
3054                  (let ((n (pool-ion-in pool-ion))
3055                        (b (pool-ion-inq pool-ion)))
3056                    (list n b)))
3057                pool-ions)
3058
3059           (map (lambda (pool-ion)
3060                  (let ((n (pool-ion-out pool-ion))
3061                        (b (pool-ion-outq pool-ion)))
3062                    (list n b)))
3063                pool-ions)
3064           ))
3065
3066         (init-dag 
3067          (map (lambda (def)
3068                 (cons (first def) (enum-freevars (second def) '() '())))
3069               init-eqs))
3070
3071         (init-order
3072          (reverse
3073           (topological-sort init-dag 
3074                             (lambda (x y) (string=? (->string x) (->string y))))))
3075
3076         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
3077
3078         )
3079
3080    (case ss-method
3081      ((kinsol)
3082       (pp indent ,nl
3083           (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(N_Vector u, N_Vector f, void* pnode)" ,nl)
3084           ( #\{)))
3085      (else
3086       (pp indent ,nl
3087           (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(const gsl_vector *u, void *pnode, gsl_vector *f)" ,nl)
3088           ( #\{)))
3089      )
3090
3091    (pp indent+
3092        (double ,(slp ", " (delete-duplicates
3093                            (map ->string 
3094                                 (filter (lambda (x) (not (member x builtin-consts)))
3095                                         (append
3096                                          init-locals
3097                                          init-order
3098                                          (map first external-eq-defs)
3099                                          (map pool-ion-in pool-ions)
3100                                          (map pool-ion-out pool-ions)
3101                                          (map first i-eqs)
3102                                          (map first steady-state-index-map) 
3103                                          (map first const-defs)
3104                                          (map first reaction-eq-defs)
3105                                          )))
3106                                 string=?))
3107                      ";")
3108        )
3109
3110    (pp indent+ ,nl
3111
3112        ("// params is a reference to the model parameters ")
3113        (const struct ,(sprintf "~A::Parameters_" sysname) 
3114               "*params = "
3115               "(" ,(sprintf "struct ~A::Parameters_ *" sysname) ")" "pnode;")
3116       
3117        ,nl
3118        )
3119
3120    (for-each
3121     (lambda (x) 
3122       (let* ((n  (first x)))
3123         (pp indent+ ,(expr->string/C++ (sprintf "params->~A" n) n))))
3124     const-defs)
3125
3126    (case ss-method
3127      ((kinsol)
3128       (for-each
3129        (lambda (def)
3130          (let* ((n   (first def)) 
3131                 (ni  (lookup-def n steady-state-index-map)))
3132            (if ni (pp indent+ ,(expr->string/C++ (ith 'u ni) n)))
3133            ))
3134        rate-eq-defs))
3135      (else
3136       (for-each
3137        (lambda (def)
3138          (let* ((n   (first def)) 
3139                 (ni  (lookup-def n steady-state-index-map)))
3140            (if ni (pp indent+ ,(expr->string/C++ (sprintf "gsl_vector_get (u, ~A)" ni) n)))
3141            ))
3142        rate-eq-defs))
3143      )
3144       
3145
3146    (let ((vi (lookup-def 'v state-index-map)))
3147      (if vi (pp indent+ ,(expr->string/C++ 0. 'v))))
3148   
3149    (for-each (lambda (def) 
3150                (pp indent+ ,(expr->string/C++ 0. (first def))))
3151              i-eqs)
3152
3153    (for-each (lambda (n)
3154                (let ((b  (lookup-def n init-eqs)))
3155                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
3156              init-order)
3157
3158    (case ss-method
3159      ((kinsol)
3160       (for-each
3161        (lambda (def)
3162          (let* ((n   (first def)) 
3163                 (ni  (lookup-def n steady-state-index-map))
3164                 (b   (second def))
3165                 (lb  (find-locals (list b))))
3166            (if ni 
3167                (begin
3168                  (if (not (null? lb))
3169                      (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
3170                  (pp indent+
3171                      (,(expr->string/C++ b (ith 'f ni))))))
3172            ))
3173        rate-eq-defs))
3174      (else
3175       (for-each
3176        (lambda (def)
3177          (let* ((n   (first def)) 
3178                 (ni  (lookup-def n steady-state-index-map))
3179                 (b   (second def))
3180                 (lb  (find-locals (list b))))
3181            (if ni 
3182                (let ((tmp (gensym 't)))
3183                  (pp indent+ (double ,tmp #\;))
3184                  (if (not (null? lb))
3185                      (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
3186                  (pp indent+
3187                      ,(expr->string/C++ b tmp)
3188                      (,(sprintf "gsl_vector_set (f,~A,~A);" ni tmp)))
3189                  ))
3190            ))
3191        rate-eq-defs))
3192      )
3193
3194    (pp indent ,nl
3195        ( "return 0;" )
3196        ( #\}))
3197
3198    ))
3199
3200
3201(define (nemo:nest-translator sys . rest)
3202
3203  (define (cid x)  (second x))
3204  (define (cn x)   (first x))
3205
3206  (let-optionals rest ((dirname ".")  (method #f) (ss-method #f) (abstol #f) (reltol #f) (maxstep #f))
3207
3208    (let ((method (or method 'gsl)))
3209
3210      (if (not (member method '(gsl cvode ida leapfrog)))
3211          (nemo:error 'nemo:nest-translator ": unknown method " method))
3212
3213  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
3214    (let ((imports  ((dis 'imports)  sys))
3215          (exports  ((dis 'exports)  sys)))
3216      (let* ((indent      0)
3217             (indent+     (+ 2 indent ))
3218
3219             (sysname     (nest-name ((dis 'sysname) sys)))
3220             (prefix      (->string sysname))
3221             (deps*       ((dis 'depgraph*)  sys))
3222             (consts      ((dis 'consts)     sys))
3223             (asgns       ((dis 'asgns)      sys))
3224             (states      ((dis 'states)     sys))
3225             (reactions   ((dis 'reactions)  sys))
3226             (defuns      ((dis 'defuns)     sys))
3227             (components  ((dis 'components) sys))
3228             
3229             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
3230             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
3231
3232             (const-defs       (filter-map
3233                                (lambda (nv)
3234                                  (and (not (member (first nv) builtin-consts))
3235                                       (let ((v1 (canonicalize-expr/C++ (second nv))))
3236                                         (list (nest-name (first nv)) v1))))
3237                                consts))
3238             
3239             (defaults             (nemo:defaults-query sys))
3240
3241             (geometry             (nemo:geometry-query sys))
3242
3243             (gate-complex-info    (nemo:gate-complex-query sys))
3244             (perm-ions       (map (match-lambda ((comp i e erev val) `(,comp ,(nest-name i) ,(nest-name e) ,erev)))
3245                                   (lookup-def 'perm-ions gate-complex-info)))
3246             (acc-ions        (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out)))))
3247                                   (lookup-def 'acc-ions gate-complex-info)))
3248             (epools          (lookup-def 'pool-ions gate-complex-info))
3249             (pool-ions       (pool-ion-name-map nest-name  epools))
3250
3251             (comprc         (any (match-lambda ((name 'membrane-tau id) (list name id)) (else #f)) components))
3252             (compcap        (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
3253             (mrc            (or (and comprc (car ((dis 'component-exports) sys (cid comprc))))
3254                                 (lookup-def 'membrane-tau defaults)
3255                                 (lookup-def 'tau_m defaults)
3256                                 (and compcap (car ((dis 'component-exports) sys (cid compcap))))
3257                                 (lookup-def 'membrane-capacitance defaults)
3258                                 (lookup-def 'C_m defaults)
3259                                 ))
3260
3261             (soma-geometry  (lookup-def 'soma geometry))
3262             (marea          (and soma-geometry (third soma-geometry)))
3263
3264             (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
3265             (synapse-info         (nemo:post-synaptic-conductance-query sys))
3266
3267
3268             (pscs           (lookup-def 'post-synaptic-conductances synapse-info))
3269
3270             (i-syns         (lookup-def 'i-synapses synapse-info))
3271               
3272             (i-gates        (lookup-def 'i-gates gate-complex-info))
3273
3274             (i-defs         (nemo:ionic-current-definitions
3275                              gate-complexes i-gates i-syns pscs marea
3276                              (lambda (x) (state-power sys x))
3277                              (lambda (x) ((dis 'component-exports) sys x))
3278                              (lambda (x) ((dis 'component-subcomps) sys x))
3279                              nest-name rhsexpr/C++ canonicalize-expr/C++
3280                              builtin-fns))
3281
3282             (i-eqs          (lookup-def 'i-eqs i-defs))
3283             (i-names        (lookup-def 'i-names i-defs))
3284
3285             (constraints    (nemo:constraint-definitions 
3286                              gate-complexes i-gates i-syns pscs marea imports
3287                              (lambda (x) (state-power sys x))
3288                              (lambda (x) (quantity-unit sys x))
3289                              (lambda (x) ((dis 'component-exports) sys x))
3290                              (lambda (x) ((dis 'component-subcomps) sys x))
3291                              nest-name))
3292
3293             (external-eq-defs   (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++
3294                                                        namespace-filter: (lambda (x) (not (equal? x 'event)))))
3295
3296             (event-external-eq-defs   (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++
3297                                                              namespace-filter: (lambda (x) (equal? x 'event))))
3298
3299             (asgn-eq-defs       (poset->asgn-eq-defs* poset sys nest-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
3300             
3301             (rate-eq-defs       (reverse (poset->rate-eq-defs* poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)))
3302             
3303             (reaction-eq-defs   (poset->reaction-eq-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++))
3304
3305             (transient-event-defs  (poset->transient-event-defs poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)) 
3306             
3307             (init-eq-defs       (poset->init-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
3308             
3309             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
3310                                      (poset->state-conserve-eq-defs poset sys nest-name nest-state-name)))
3311             
3312             (imports-sans-v (filter (lambda (x) (not (equal? 'v (first x)))) imports))
3313
3314             (v-eq    (and (not (null? i-names))
3315                           (let ((istim "(node.B_.I_stim_)" )) 
3316                             (cond
3317
3318                              ((and mrc marea)
3319                               (list 'v (rhsexpr/C++ 
3320                                         `(/ (+ (* ,istim (/ 100. ,marea)) 
3321                                                (* -1e3 ,(sum i-names))) ,mrc))))
3322                              (marea
3323                               (list 'v (rhsexpr/C++ 
3324                                         `(+ (* ,istim (/ 100. ,marea))
3325                                             (* -1e-3 ,(sum i-names))))))
3326                              (mrc
3327                               (list 'v (rhsexpr/C++ `(/ (+ ,istim (* -1e-3 ,(sum i-names))) ,mrc))))
3328                             
3329                              (else
3330                               (list 'v (rhsexpr/C++ `(+ ,istim (* -1e-3 ,(sum i-names))))))
3331                              ))
3332                           ))
3333
3334             
3335             (state-index-map  (let ((acc (fold (lambda (def ax)
3336                                                  (let ((st-name (first def)))
3337                                                    (list (+ 1 (first ax)) 
3338                                                          (cons `(,st-name ,(first ax)) (second ax)))))
3339                                                (list 0 (list)) 
3340                                                (if (lookup-def 'v imports)
3341                                                    (cons (list 'v) rate-eq-defs)
3342                                                    rate-eq-defs)
3343                                                )))
3344                                 (second acc)))
3345             
3346             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
3347                                                         (let ((st-name (first def)))
3348                                                           (if (not (alist-ref st-name init-eq-defs))
3349                                                               (list (+ 1 (first ax)) 
3350                                                                     (cons `(,st-name ,(first ax)) (second ax)))
3351                                                               ax)))
3352                                                       (list 0 (list)) 
3353                                                       rate-eq-defs)))
3354                                        (second acc)))
3355             
3356             (dfenv 
3357              (map (lambda (x) (let ((n (first x))) (list n (nest-name (sprintf "d_~A" n )))))
3358                   defuns))
3359
3360             )
3361       
3362       
3363       
3364        (for-each
3365         (lambda (a)
3366           (let ((acc-ion   (car a)))
3367             (if (assoc acc-ion perm-ions)
3368                 (nemo:error 'nemo:nest-translator 
3369                             ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
3370         acc-ions)
3371
3372
3373        (let ((cpp-output  (open-output-file (make-output-fname dirname prefix ".cpp")))
3374              (hpp-output  (open-output-file (make-output-fname dirname prefix ".h"))))
3375         
3376          ;; include files and other prelude definitions
3377          (with-output-to-port cpp-output
3378            (lambda ()
3379              (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
3380              (output-prelude sysname indent)
3381              ))
3382         
3383          (with-output-to-port cpp-output
3384            (lambda () (pp indent ,nl "namespace nest {" ,nl)))
3385         
3386          (case method
3387            ((cvode ida) 
3388             (with-output-to-port cpp-output
3389               (lambda ()
3390                 (pp indent ,sundials-prelude)))))
3391
3392
3393          (if (not (null? steady-state-index-map))
3394              (with-output-to-port cpp-output
3395                (lambda ()
3396                  (pp indent ,(fsolve-prelude ss-method)))
3397                ))
3398
3399          ;; user-defined functions
3400          (let ((define-fn  (make-define-fn sysname))
3401                (define-fn-header (make-define-fn-header sysname)))
3402           
3403            (for-each (lambda (fndef) 
3404                        (if (not (member (car fndef) builtin-fns))
3405                            (with-output-to-port cpp-output
3406                              (lambda ()
3407                                (apply define-fn-header (cons indent fndef))
3408                                (pp indent ,nl)))
3409                            ))
3410                      defuns)
3411
3412            (for-each (lambda (fndef) 
3413                        (if (not (member (car fndef) builtin-fns))
3414                            (with-output-to-port cpp-output
3415                              (lambda ()
3416                                (apply define-fn (cons indent fndef))
3417                                (pp indent ,nl)))
3418                            ))
3419                      defuns)
3420            )
3421         
3422         
3423          ;; derivative function
3424          (with-output-to-port cpp-output
3425            (lambda ()
3426              (output-dy sysname method imports-sans-v const-defs state-index-map 
3427                         external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3428                         pool-ions i-eqs v-eq constraints
3429                         indent indent+)
3430              ))
3431
3432          ;; residual function (used by IDA only)
3433          (with-output-to-port cpp-output
3434            (lambda ()
3435              (output-residual sysname method imports-sans-v const-defs state-index-map 
3436                               external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3437                               pool-ions i-eqs v-eq constraints
3438                               indent indent+)
3439              ))
3440
3441          ;; steady state function
3442          (if (not (null? steady-state-index-map))
3443              (with-output-to-port cpp-output
3444                (lambda ()
3445                  (output-steadystate sysname method ss-method state-index-map steady-state-index-map 
3446                                      imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3447                                      reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+)
3448                  (pp indent ,nl)
3449                  )))
3450
3451          ;; Jacobian function
3452          (with-output-to-port cpp-output
3453            (lambda ()
3454              (output-jac  sysname method state-index-map pool-ions i-eqs v-eq indent indent+)
3455              ))
3456
3457          ;;  system recordables
3458          (with-output-to-port cpp-output
3459            (lambda ()
3460              (output-recordables
3461               sysname state-index-map 
3462               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3463               reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
3464              (pp indent ,nl)
3465              ))
3466         
3467          ;; system parameters
3468          (with-output-to-port cpp-output
3469            (lambda ()
3470              (output-parameters sysname imports-sans-v const-defs 
3471                                 external-eq-defs defaults indent indent+)
3472              ))
3473         
3474          ;; initial values function
3475          (with-output-to-port cpp-output
3476            (lambda ()
3477              (output-init sysname method ss-method state-index-map steady-state-index-map 
3478                           imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3479                           reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+)
3480              (pp indent ,nl)
3481              ))
3482         
3483         
3484          ;; accessors and modifiers for states and parameters
3485          (with-output-to-port cpp-output
3486            (lambda ()
3487              (output-accessors+modifiers
3488               sysname imports-sans-v defaults state-index-map 
3489               const-defs asgn-eq-defs rate-eq-defs 
3490               reaction-eq-defs i-eqs pool-ions perm-ions 
3491               constraints
3492               indent indent+)
3493              (pp indent ,nl)
3494              ))
3495         
3496         
3497          ;; buffer and node initialization
3498          (with-output-to-port cpp-output
3499            (lambda ()
3500              (output-buffers
3501               sysname method state-index-map 
3502               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3503               reaction-eq-defs i-eqs pool-ions perm-ions 
3504               synapse-info defaults
3505               indent indent+)
3506              (pp indent ,nl)
3507              ))
3508         
3509
3510          (with-output-to-port cpp-output
3511            (lambda ()
3512              ((case method
3513                 ((leapfrog) 
3514                  output-nodes-leapfrog)
3515                 ((gsl)     
3516                  output-nodes-gsl)
3517                 ((cvode)
3518                  output-nodes-cvode)
3519                 ((ida)
3520                  output-nodes-ida)
3521                 )
3522               sysname method state-index-map 
3523               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3524               reaction-eq-defs i-eqs pool-ions perm-ions 
3525               synapse-info defaults abstol reltol maxstep
3526               indent indent+)
3527              (pp indent ,nl)
3528              ))
3529         
3530          ;; spike threshold detect (cvode and ida methods only)
3531          (case method
3532            ((cvode)
3533              (with-output-to-port cpp-output
3534                (lambda ()
3535                  (output-cvode-event
3536                   sysname method imports const-defs state-index-map 
3537                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3538                   defaults i-eqs v-eq indent indent+)
3539                  (pp indent ,nl)
3540                  )))
3541            ((ida)
3542              (with-output-to-port cpp-output
3543                (lambda ()
3544                  (output-ida-event
3545                   sysname method imports const-defs state-index-map 
3546                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3547                   defaults i-eqs v-eq indent indent+)
3548                  (pp indent ,nl)
3549                  )))
3550            )
3551
3552          ;; synaptic event transients
3553
3554          (with-output-to-port cpp-output
3555            (lambda ()
3556              (output-synapse-transients
3557               sysname method state-index-map
3558               const-defs transient-event-defs event-external-eq-defs 
3559               synapse-info imports constraints
3560               indent indent+)
3561              (pp indent ,nl)
3562              ))
3563
3564          ;; state update
3565          (with-output-to-port cpp-output
3566            (lambda ()
3567              (output-update
3568               sysname method state-index-map 
3569               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3570               reaction-eq-defs transient-event-defs i-eqs
3571               pool-ions perm-ions synapse-info 
3572               event-external-eq-defs
3573               defaults imports abstol reltol
3574               indent indent+)
3575              (pp indent ,nl)
3576              ))
3577         
3578          ;; spike handler
3579          (with-output-to-port cpp-output
3580            (lambda ()
3581              (output-event-handle
3582               sysname state-index-map 
3583               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3584               reaction-eq-defs i-eqs pool-ions perm-ions 
3585               synapse-info
3586               indent indent+)
3587              (pp indent ,nl)
3588              ))
3589         
3590          (with-output-to-port cpp-output
3591            (lambda () (pp indent "}" ,nl)))
3592
3593             (with-output-to-port hpp-output
3594               (lambda ()
3595                 (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
3596                 (output-header
3597                  sysname method ss-method state-index-map steady-state-index-map defaults
3598                  external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3599                  reaction-eq-defs i-eqs pool-ions perm-ions 
3600                  synapse-info
3601                  indent indent+)
3602                 (pp indent ,nl)
3603                 ))
3604             
3605             (close-output-port cpp-output)
3606             (close-output-port hpp-output)
3607             
3608             ))
3609      ))
3610  ))
3611
3612))
Note: See TracBrowser for help on using the repository browser.