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

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

nemo: implemented gsl steady state solver in nest backend

File size: 102.2 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 method)
513  (case method
514    ((cvode ida)
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 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           ("#include <kinsol/kinsol.h>           /* prototypes for KINSOL fcts. */")
697           ("#include <kinsol/kinsol_dense.h>     /* prototype for KINDense*/")
698           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
699           ))
700      ((ida) 
701       (pp indent
702           ("#include <sundials/sundials_types.h> /* definition of type realtype */")
703           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
704           ("#include <ida/ida.h>                 /* prototypes for IDA fcts., consts. */")
705           ("#include <ida/ida_dense.h>")
706           ("#include <kinsol/kinsol.h>           /* prototypes for KINSOL fcts. */")
707           ("#include <kinsol/kinsol_dense.h>     /* prototype for KINDense */")
708           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
709           ))
710      ((gsl) 
711       (pp indent
712           ("#include <gsl/gsl_errno.h>")
713           ("#include <gsl/gsl_matrix.h>")
714           ("#include <gsl/gsl_sf_exp.h>")
715           ("#include <gsl/gsl_odeiv2.h>")
716           ("#include <gsl/gsl_vector.h>")
717           ("#include <gsl/gsl_multiroots.h>")
718           ("#define Ith(v,i)    (v[i])")
719           ))
720      )
721
722    (pp indent ,nl)
723
724    (pp indent "namespace nest {" ,nl)
725 
726    (case method
727      ((cvode) 
728       (pp indent ,cvode-prelude-header))
729      ((ida) 
730       (pp indent ,ida-prelude-header)))
731
732    (if (not (null? steady-state-index-map))
733        (case method
734          ((cvode ida)
735           (pp indent ,kinsol-prelude-header))))
736     
737 
738    (case method
739      ((leapfrog) 
740       (pp indent 
741           ("typedef int (*rhsfn_t)(double t, const double *y, double *ydot, void *user_data);") ,nl
742           (extern "\"C\""
743                   int ,(sprintf "~A_dynamics" sysname)
744                   "(double, const double *, double *, void*)" #\;) ,nl))
745      ((cvode) 
746       (pp indent (extern "\"C\""
747                          int ,(sprintf "~A_dynamics" sysname)
748                          "(double, const N_Vector, N_Vector, void*)" #\;) ,nl
749                   (extern "\"C\""
750                           int ,(sprintf "~A_event" sysname)
751                           "(double, N_Vector, double *, void*)" #\;) ,nl
752                           ))
753      ((ida) 
754       (pp indent (extern "\"C\""
755                          int ,(sprintf "~A_residual" sysname)
756                          "(double, N_Vector, N_Vector, N_Vector, void*)" #\;) ,nl
757                   (extern "\"C\""
758                           int ,(sprintf "~A_event" sysname)
759                           "(double, N_Vector, N_Vector, double *, void*)" #\;) ,nl
760                           ))
761      (else
762       (pp indent (extern "\"C\""
763                          int ,(sprintf  "~A_dynamics" sysname)
764                          "(double, const double*, double*, void*)" #\;) ,nl))
765      )
766
767    (if (not (null? steady-state-index-map))
768        (case method
769          ((cvode ida)
770           (pp indent (extern "\"C\""
771                              int ,(sprintf "~A_steadystate" sysname)
772                              "(N_Vector, N_Vector, void*)" #\;) ,nl)
773           )
774          (else
775           (pp indent (extern "\"C\""
776                              int ,(sprintf "~A_steadystate" sysname)
777                              "(const gsl_vector *, void *, gsl_vector *)" #\;) ,nl)
778           )
779          ))
780
781    (pp indent 
782        (,(sprintf "class ~A:" sysname) "public Archiving_Node { "))
783
784    (pp indent+ ( ,#<#EOF
785public: 
786 ~#{sysname} ();
787 #{sysname} (const #{sysname} &);
788 #{sysname} ();
789EOF
790    ))
791
792    (pp indent ( #<<EOF
793    using Node::connect_sender;
794    using Node::handle;
795
796    port check_connection(Connection&, port);
797   
798    void handle(SpikeEvent &);
799    void handle(CurrentEvent &);
800    void handle(DataLoggingRequest &);
801   
802    port connect_sender(SpikeEvent &, port);
803    port connect_sender(CurrentEvent &, port);
804    port connect_sender(DataLoggingRequest &, port);
805   
806    void get_status(DictionaryDatum &) const;
807    void set_status(const DictionaryDatum &);
808   
809    void init_node_(const Node& proto);
810    void init_state_(const Node& proto);
811    void init_buffers_();
812    void calibrate();
813   
814    void update(Time const &, const long_t, const long_t);
815
816EOF
817))
818
819  (pp indent (,#<#EOF
820
821    /**
822     * Minimal spike receptor type.
823     * @note Start with 1 so we can forbid port 0 to avoid accidental
824     *       creation of connections with no receptor type set.
825     */
826    static const port MIN_SPIKE_RECEPTOR = 1;
827
828    /** 
829     * Spike receptors.
830     */
831    enum SpikeSynapseTypes { #{(slp ", " (append
832                                         (if (pair? pscs)
833                                             (cons (sprintf "~A_SPIKE_RECEPTOR=MIN_SPIKE_RECEPTOR" (first (car pscs)))
834                                                   (map (lambda (x) (sprintf "~A_SPIKE_RECEPTOR" (first x))) (cdr pscs))) '())
835                                         (list "SUP_SPIKE_RECEPTOR ")))}
836                             };
837
838
839EOF
840))
841 
842  (case method
843    ((cvode) 
844     (pp indent+ (friend int ,(sprintf "~A_dynamics" sysname)
845                         "(double, const N_Vector, N_Vector, void*)" #\; ) ,nl
846                 (friend int ,(sprintf "~A_event" sysname)
847                         "(double, N_Vector, double *, void*)" #\;) ,nl
848                         ))
849    ((ida) 
850     (pp indent+ (friend int ,(sprintf "~A_residual" sysname)
851                         "(double, N_Vector, N_Vector, N_Vector, void*)" #\; ) ,nl
852                 (friend int ,(sprintf "~A_event" sysname)
853                         "(double, N_Vector, N_Vector, double *, void*)" #\;) ,nl
854                         ))
855    (else
856     (pp indent+ (friend int ,(sprintf "~A_dynamics" sysname)
857                         "(double, const double*, double*, void*)" #\;) ,nl))
858    )
859
860 
861  (if (not (null? steady-state-index-map))
862      (case method
863        ((cvode ida)
864         (pp indent (friend  int ,(sprintf "~A_steadystate" sysname)
865                            "(N_Vector, N_Vector, void*)" #\;) ,nl))
866        (else
867         (pp indent (friend  int ,(sprintf "~A_steadystate" sysname)
868                            "(const gsl_vector *, void *, gsl_vector *)" #\;) ,nl))
869        ))
870         
871
872  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
873    (for-each
874     (lambda (psc)
875       (pp indent+ (,(sprintf "int ~A_transients (long_t lag);" (first psc)))))
876     pscs))
877 
878
879  (pp indent (,#<#EOF
880    // The next two classes need to be friends to access the State_ class/member
881    friend class RecordablesMap<#{sysname}>;
882    friend class UniversalDataLogger<#{sysname}>;
883EOF
884   ))
885
886 
887  (pp indent+ ,nl "struct Parameters_ { ")
888
889  (for-each
890   (lambda (x) (pp indent+ (,(sprintf "double ~A;" x))))
891   (append
892    (map (compose nest-name first) const-defs)
893    (map (compose nest-name first) defaults)))
894
895  (pp indent+ 
896      ("Parameters_();")
897      ("void get(DictionaryDatum&) const;")
898      ("void set(const DictionaryDatum&);"))
899
900  (pp indent+ "};")
901
902  (pp indent+ ,nl "struct State_ { ")
903
904  (pp indent+ ,nl 
905      "enum StateVecElems {"
906      (,(slp ", " (map (lambda (x) 
907                         (let ((n (string-upcase (->string (first x))))
908                               (ni (second x)))
909                           (sprintf "~A = ~A" n ni)))
910                       state-index-map)))
911      "};")
912
913  (pp indent+ ,nl
914      (,(sprintf "double y_[~A];" (length state-index-map)))
915
916      "State_(const Parameters_& p);" 
917      "State_(const State_& s);"
918      "State_& operator=(const State_& s);"
919      "void get(DictionaryDatum&) const;"
920      "void set(const DictionaryDatum&, const Parameters_&);")
921
922  (case method
923    ((gsl leapfrog) 
924     (pp indent+ ,nl "int_t r_; /* refractory counter */")))
925
926  (pp indent+ "};")
927
928  (case method
929    ((cvode ida) 
930     (pp indent+ ,nl "struct Variables_ {};"))
931    (else
932     (pp indent+ ,nl "struct Variables_ { int_t RefractoryCounts_; double U_old_; /* for spike-detection */ };"))
933    )
934
935  (pp indent+ ,nl "struct Buffers_ { ")
936
937  (pp indent+ ,nl 
938      (,(sprintf "Buffers_(~A&);" sysname))
939      (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
940      (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))
941
942  (pp indent+ ,nl
943
944      (,(case method
945
946    ((leapfrog)
947#<<EOF
948  rhsfn_t sys_; // pointer to function that computes dv/dt
949  double *u, *v, *dvdt;  // intermediate state vectors used by leapfrog method
950  unsigned int  N;  // size of state vector used by leapfrog method
951EOF
952)
953
954    ((cvode)
955#<<EOF
956  N_Vector y;    //!< current state vector used by CVode
957  void *   sys_;  //!< CVode control structure
958EOF
959)
960    ((ida)
961#<<EOF
962  N_Vector y, y1;    //!< current state vector used by IDA
963  N_Vector yp;    //!< derivatives vector used by IDA
964  void *   sys_;  //!< IDA control structure
965EOF
966)
967
968    ((gsl)
969#<<EOF
970  gsl_odeiv2_step*    s_;    //!< stepping function
971  gsl_odeiv2_control* c_;    //!< adaptive stepsize control function
972  gsl_odeiv2_evolve*  e_;    //!< evolution function
973  gsl_odeiv2_system   sys_;  //!< struct describing system
974  unsigned int N;  // size of state vector used by Jacobian
975  double *u, *jac;  // intermediate state vectors used for Jacobian approximation
976EOF
977)
978
979    )))
980
981
982  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
983    (for-each
984     (lambda (psc)
985       (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
986     pscs))
987
988  (pp indent+ ,nl
989#<<EOF
990  RingBuffer currents_;
991
992  double_t step_;           //!< step size in ms
993  double   IntegrationStep_;//!< current integration time step, updated by solver
994
995  /** 
996   * Input current injected by CurrentEvent.
997   * This variable is used to transport the current applied into the
998   * _dynamics function computing the derivative of the state vector.
999   * It must be a part of Buffers_, since it is initialized once before
1000   * the first simulation, but not modified before later Simulate calls.
1001   */
1002  double_t I_stim_;
1003EOF
1004)
1005  (pp indent+ "};")
1006
1007  (pp indent+
1008      "template <State_::StateVecElems elem>"
1009      "double_t get_y_elem_() const { return S_.y_[elem]; }"
1010      "Parameters_ P_;"
1011      "State_      S_;"
1012      "Variables_  V_;"
1013      "Buffers_    B_;"
1014      (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
1015      "}; ")
1016
1017  (pp indent+ (,#<#EOF
1018  inline
1019  port #{sysname}::check_connection(Connection& c, port receptor_type)
1020  {
1021    SpikeEvent e;
1022    e.set_sender(*this);
1023    c.check_event(e);
1024    return c.get_target()->connect_sender(e, receptor_type);
1025  }
1026
1027  inline
1028  port #{sysname}::connect_sender(SpikeEvent&, port receptor_type)
1029  {
1030    if ( receptor_type < MIN_SPIKE_RECEPTOR || receptor_type >= SUP_SPIKE_RECEPTOR )
1031    {
1032      if ( receptor_type < 0 || receptor_type >= SUP_SPIKE_RECEPTOR )
1033        throw UnknownReceptorType(receptor_type, get_name());
1034      else
1035        throw IncompatibleReceptorType(receptor_type, get_name(), "SpikeEvent");
1036    }
1037    return receptor_type;
1038  }
1039 
1040  inline
1041  port #{sysname}::connect_sender(CurrentEvent&, port receptor_type)
1042  {
1043    if (receptor_type != 0)
1044      throw UnknownReceptorType(receptor_type, get_name());
1045    return 0;
1046  }
1047
1048  inline
1049  port #{sysname}::connect_sender(DataLoggingRequest& dlr, 
1050                                      port receptor_type)
1051  {
1052    if (receptor_type != 0)
1053      throw UnknownReceptorType(receptor_type, get_name());
1054    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
1055  }
1056
1057  inline
1058    void #{sysname}::get_status(DictionaryDatum &d) const
1059  {
1060    P_.get(d);
1061    S_.get(d);
1062    Archiving_Node::get_status(d);
1063
1064    (*d)[names::recordables] = recordablesMap_.get_list();
1065
1066    def<double_t>(d, names::t_spike, get_spiketime_ms());
1067
1068    DictionaryDatum receptor_dict_ = new Dictionary();
1069
1070    #{(string-concatenate
1071       (map (lambda (psc) 
1072              (sprintf "    (*receptor_dict_)[Name(\"~A\")]  = ~A_SPIKE_RECEPTOR;~%"
1073                       (first psc) (first psc)))
1074            pscs))}
1075               
1076    (*d)[names::receptor_types] = receptor_dict_;
1077  }
1078
1079  inline
1080    void #{sysname}::set_status(const DictionaryDatum &d)
1081  {
1082    Parameters_ ptmp = P_;  // temporary copy in case of errors
1083    ptmp.set(d);                       // throws if BadProperty
1084    State_      stmp = S_;  // temporary copy in case of errors
1085    stmp.set(d, ptmp);                 // throws if BadProperty
1086
1087    // We now know that (ptmp, stmp) are consistent. We do not 
1088    // write them back to (P_, S_) before we are also sure that 
1089    // the properties to be set in the parent class are internally 
1090    // consistent.
1091    Archiving_Node::set_status(d);
1092
1093    // if we get here, temporaries contain consistent set of properties
1094    P_ = ptmp;
1095    S_ = stmp;
1096
1097    calibrate();
1098  }
1099EOF
1100  ))
1101
1102  (pp indent "}" ,nl)
1103  )
1104)
1105
1106
1107
1108(define (output-jac sysname method state-index-map pool-ions i-eqs v-eq indent indent+)
1109
1110;; Diagonal Jacobian approximation: (f(s+.01) - f(s))/.001
1111
1112;; CVODE:
1113;;  typedef int (*CVDlsDenseJacFn)(long int N, realtype t,
1114;;                               N_Vector y, N_Vector fy,
1115;;                               DlsMat Jac, void *user_data,
1116;;                               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
1117
1118;; GSL:
1119;; int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
1120
1121
1122
1123  (case method
1124
1125    ((gsl)
1126     (pp indent ,nl
1127         (,(sprintf "extern \"C\" int ~A_jacobian" sysname)
1128          (,(slp ", " `("double t"
1129                        "const double y[]"
1130                        "double *dfdy"
1131                        "double dfdt[]"
1132                        "void* pnode"
1133                        ))))
1134          #\{
1135
1136          ("// cast the node ptr to an object of the proper type")
1137          ("assert(pnode);")
1138          ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1139          (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1140
1141          ("// state is a reference to the model state ")
1142          (struct ,(sprintf "~A::Buffers_" sysname) "*b;")
1143          ("b = &(vnode.B_);")
1144          (,#<#EOF
1145   for (int i = 0; i < b->N; i++)
1146   {
1147       b->u[i] = y[i] + 0.01;
1148   }
1149   #{sysname}_dynamics(t, b->u, b->jac, pnode);
1150   for (int i = 0; i < b->N; i++)
1151   {
1152       dfdt[i*b->N + i] = (b->jac[i] - dfdy[i]) / .001;
1153   }
1154   return 0;
1155EOF
1156)
1157   #\}
1158   ))
1159))
1160
1161
1162(define (output-cvode-event sysname method imports const-defs state-index-map 
1163           external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
1164           defaults i-eqs v-eq indent indent+)
1165
1166  (let ((vi (lookup-def 'v state-index-map)))
1167   
1168    (pp indent ,nl
1169        (,(sprintf "extern \"C\" int ~A_event" sysname)
1170         (,(slp ", " `("double t"
1171                       "N_Vector y"
1172                       "double *g"
1173                       "void* pnode"
1174                       )))
1175         #\{
1176         ))
1177
1178      (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;")
1179
1180      ,nl
1181
1182      ("// S is shorthand for the type that describes the model state ")
1183      (typedef ,(sprintf "~A::State_" sysname) "S;")
1184           
1185      ,nl
1186           
1187      ("// cast the node ptr to an object of the proper type")
1188      ("assert(pnode);")
1189      ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1190     
1191      ,nl
1192     
1193      ("// params is a reference to the model parameters ")
1194      (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
1195      ("params = &(node.P_);")
1196     
1197      ,nl
1198     
1199      ("// state is a reference to the model state ")
1200      (const struct ,(sprintf "~A::State_" sysname) "*state;")
1201      ("state = &(node.S_);")
1202     
1203      )
1204     
1205      (let ((vt (lookup-def 'V_t defaults)))
1206        (if (and vi vt)
1207            (pp indent+
1208                (,(expr->string/C++ (ith 'y vi) 'v))
1209                (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt))
1210                (,(expr->string/C++ '(- v vt) "g[0]"))
1211                ,nl
1212                ("return 0; "))
1213            (pp indent+
1214                (,(expr->string/C++ -1.0 "g[0]")))
1215            ))
1216
1217     
1218      (pp indent (#\}))
1219     
1220      ))
1221
1222
1223(define (output-ida-event sysname method imports const-defs state-index-map 
1224           external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
1225           defaults i-eqs v-eq indent indent+)
1226
1227  (let ((vi (lookup-def 'v state-index-map)))
1228   
1229    (pp indent ,nl
1230        (,(sprintf "extern \"C\" int ~A_event" sysname)
1231         (,(slp ", " `("double t"
1232                       "N_Vector y"
1233                       "N_Vector yp"
1234                       "double *g"
1235                       "void* pnode"
1236                       )))
1237         #\{
1238         ))
1239
1240      (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;")
1241
1242      ,nl
1243
1244      ("// S is shorthand for the type that describes the model state ")
1245      (typedef ,(sprintf "~A::State_" sysname) "S;")
1246           
1247      ,nl
1248           
1249      ("// cast the node ptr to an object of the proper type")
1250      ("assert(pnode);")
1251      ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
1252     
1253      ,nl
1254     
1255      ("// params is a reference to the model parameters ")
1256      (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
1257      ("params = &(node.P_);")
1258     
1259      ,nl
1260     
1261      ("// state is a reference to the model state ")
1262      (const struct ,(sprintf "~A::State_" sysname) "*state;")
1263      ("state = &(node.S_);")
1264     
1265      )
1266     
1267      (let ((vt (lookup-def 'V_t defaults)))
1268        (if (and vi vt)
1269            (pp indent+
1270                (,(expr->string/C++ (ith 'y vi) 'v))
1271                (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt))
1272                (,(expr->string/C++ '(- v vt) "g[0]"))
1273                ,nl
1274                ("return 0; "))
1275            (pp indent+
1276                (,(expr->string/C++ -1.0 "g[0]")))
1277            ))
1278
1279     
1280      (pp indent (#\}))
1281     
1282      ))
1283
1284
1285
1286
1287(define (cvode-solve spike-handle)
1288#<#EOF
1289        int N = NV_LENGTH_S (B_.y);
1290        tout = (current_steps+lag)*h;
1291
1292        // adaptive step integration
1293        while (tt < tout)
1294        {
1295          const int status = CVode (B_.sys_, tout, B_.y, &tt, CV_NORMAL);
1296
1297          switch (status)
1298            {
1299                case CV_SUCCESS:      continue;
1300                case CV_ROOT_RETURN:  #{spike-handle};
1301                default:              throw CVodeSolverFailure (get_name(), 0);
1302            }
1303        }
1304
1305        for (int i = 0; i < N; i++)
1306        {
1307           S_.y_[i] = Ith(B_.y,i);
1308        }
1309EOF
1310)
1311
1312
1313(define (ida-solve spike-handle)
1314#<#EOF
1315        int N = NV_LENGTH_S (B_.y);
1316        tout = (current_steps+lag)*h;
1317
1318        // adaptive step integration
1319        while (tt < tout)
1320        {
1321          const int status = IDASolve (B_.sys_, tout, &tt, B_.y, B_.yp, IDA_NORMAL);
1322
1323          switch (status)
1324            {
1325                case IDA_SUCCESS:      continue;
1326                case IDA_ROOT_RETURN:  #{spike-handle};
1327                case IDA_TSTOP_RETURN: break;
1328                default:               throw IDASolverFailure (get_name(), 0);
1329            }
1330        }
1331
1332        for (int i = 0; i < N; i++)
1333        {
1334           S_.y_[i] = Ith(B_.y,i);
1335        }
1336EOF
1337)
1338
1339
1340(define (gsl-solve vi spike-handle)
1341#<#EOF
1342
1343   V_.U_old_ = S_.y_[#{vi}];
1344
1345   while (tt < h)
1346   {
1347
1348     const int status = gsl_odeiv2_evolve_apply
1349                         (B_.e_, B_.c_, B_.s_, 
1350                          &B_.sys_,              // system of ODE
1351                          &tt,                   // from t...
1352                          h,             // ...to t=t+h
1353                         &B_.IntegrationStep_ , // integration window (written on!)
1354                          S_.y_);               // neuron state
1355
1356          if ( status != GSL_SUCCESS )
1357            throw GSLSolverFailure(get_name(), status);
1358
1359   }
1360
1361   #{spike-handle}
1362EOF
1363)
1364
1365
1366(define (leapfrog-solve vi spike-handle)
1367#<#EOF
1368        V_.U_old_ = S_.y_[#{vi}];
1369
1370        tt = (current_steps+lag)*h;
1371        tout = (current_steps+lag+1)*h;
1372
1373        while (tt < tout)
1374        {
1375           B_.sys_ (tt, S_.y_, B_.dvdt, reinterpret_cast<void*>(this));
1376
1377           for (int i = 0; i < B_.N; i++)
1378           {
1379              B_.u[i] = S_.y_[i] + 0.5 * h * B_.dvdt[i];
1380           }
1381
1382           B_.sys_ (tt+0.5*h, B_.u, B_.dvdt, reinterpret_cast<void*>(this));
1383
1384           for (int i = 0; i < B_.N; i++)
1385           { 
1386              B_.v[i] = S_.y_[i] + h * B_.dvdt[i];
1387           }
1388
1389           tt += h;
1390        }
1391
1392        for (int i = 0; i < B_.N; i++)
1393        {
1394           S_.y_[i] = B_.v[i];
1395        }
1396 
1397        #{spike-handle}
1398EOF
1399)
1400
1401
1402
1403(define (output-synapse-transients sysname method state-index-map 
1404                                   const-defs transient-event-defs event-external-eq-defs 
1405                                   synapse-info imports constraints indent indent+)
1406
1407  (let (
1408        (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
1409                                    (list (nest-name n) v)))
1410                      (lookup-def 'c-units constraints)))
1411
1412        (isyns  (lookup-def 'i-synapses synapse-info))
1413        (pscs   (lookup-def 'post-synaptic-conductances synapse-info))
1414        (psc-transients (map (lambda (lst) (map nest-name lst)) 
1415                             (lookup-def 'psc-transients synapse-info)))
1416        (getstate (case method
1417                    ((cvode ida) (lambda (i) (sprintf "Ith(B_.y,~A)" i)))
1418                    (else (lambda (i) (sprintf "Ith(S_.y_,~A)" i)))))
1419        )
1420
1421
1422    (if (not (= (length event-external-eq-defs) (length pscs)))
1423        (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" 
1424               event-external-eq-defs pscs))
1425     
1426    (for-each (lambda (psc transients)
1427
1428                (let* (
1429                       (ltransient-event-defs
1430                        (filter (lambda (x) (member (first x) transients))
1431                                transient-event-defs))
1432
1433                       (levent-external-eq-def
1434                        (car
1435                         (fold (lambda (def ax)
1436                                 (let* ((n (nest-name (first def)) )
1437                                        (b (second def))
1438                                        (events (let ((fvs (enum-freevars b '() '())))
1439                                                  (filter (lambda (x) (member (first x) fvs)) event-external-eq-defs)))
1440                                        )
1441                                   (append events ax)))
1442                               '() ltransient-event-defs)))
1443
1444                       (lconsts
1445                        (delete-duplicates
1446                         (fold (lambda (def ax)
1447                                 (let* ((n (nest-name (first def)) )
1448                                        (b (second def))
1449                                        (consts (let ((fvs (enum-freevars b '() '())))
1450                                                  (filter (lambda (x) (member (first x) fvs)) const-defs)))
1451                                       )
1452                                   (append consts ax)))
1453                               '() ltransient-event-defs)
1454                         (lambda (x y) (equal? (first x) (first y))))
1455                        ))
1456                 
1457                  (pp indent (,(sprintf "inline int ~A::~A_transients (long_t lag) {" sysname (first psc))))
1458                 
1459                  (if (not (null? lconsts)) 
1460                      (pp indent+ (double ,(slp ", " (map (compose nest-name car) lconsts)) ";")))
1461
1462                  (let ((n (second levent-external-eq-def)))
1463                    (pp indent+
1464                        (,(sprintf "double_t ~A;" n))
1465                        (,(sprintf "~A = B_.spike_~A.get_value(lag);" 
1466                                   n (nest-name (second psc))))
1467                        )
1468                    )
1469                 
1470                  (pp indent+
1471                      (,(sprintf "if (~A > 0.0)" (second levent-external-eq-def)))
1472                      (#\{))
1473                     
1474                  (let* ((n      (nest-name (first levent-external-eq-def)))
1475                         (nu     (lookup-def n c-units))
1476                         (nscale (and nu (nemo:unit-scale nu)))
1477                         (b      (second levent-external-eq-def))
1478                         )
1479                    (pp (+ 2 indent+)
1480                        (,(sprintf "double_t ~A;" n))
1481                        (,(expr->string/C++ (if nscale `(* ,nscale ,b) b) n)))
1482                    )
1483                 
1484                  (for-each (lambda (def)
1485                              (let* ((n  (nest-name (first def)) )
1486                                     (ni (lookup-def n state-index-map))
1487                                     (b  (second def))
1488                                     (consts (let ((fvs (enum-freevars b '() '())))
1489                                               (filter (lambda (x) (member (first x) fvs)) const-defs)))
1490                                     )
1491                               
1492                              (if (not (null? consts)) 
1493                                  (for-each (lambda (x) 
1494                                              (let ((n (first x)))
1495                                                (pp (+ 2 indent+) (,(nest-name n) = ,(sprintf "P_.~A;" (nest-name n))))))
1496                                            consts)
1497                                  )
1498                             
1499                              (pp (+ 2 indent+ )
1500                                  (,(sprintf "double_t ~A;" n))
1501                                  ,(if ni (expr->string/C++ (getstate ni) n) '())
1502                                  ,(expr->string/C++ b n)
1503                                  ,(if ni (expr->string/C++ n  (getstate ni)) '())
1504                                  )
1505                              ))
1506
1507                            ltransient-event-defs)
1508                 
1509                 
1510                  (pp (+ 2 indent+) "return 1;");
1511                  (pp indent+ (#\}))
1512                  (pp indent+ "return 0;");
1513                  (pp indent (#\}))
1514                  ))
1515
1516               pscs psc-transients)
1517
1518    ))
1519
1520
1521(define (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports indent)
1522
1523  (let (
1524        (isyns  (lookup-def 'i-synapses synapse-info))
1525        (pscs   (lookup-def 'post-synaptic-conductances synapse-info))
1526        )
1527
1528    (pp indent (,(sprintf "int events = 0;")))
1529
1530    (if (not (= (length event-external-eq-defs) (length pscs)))
1531        (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" 
1532               event-external-eq-defs pscs))
1533     
1534    (for-each (lambda (psc)
1535                (pp indent
1536                    (,(sprintf "events = events + ~A_transients (lag);" (first psc))))
1537                )
1538              pscs)
1539
1540    (case method
1541      ((cvode)
1542        (pp indent 
1543            ( "/* Reinitializes CVode state if any synaptic events have occurred */")
1544            ("if (events > 0)")
1545            (#\{)
1546            ("  int status = CVodeReInit (B_.sys_, tt, B_.y);")
1547            ("  if (check_flag(&status, \"CVodeReInit\", 1)) throw CVodeSolverFailure (get_name(), status);")
1548            (#\})
1549            ))
1550      ((ida)
1551        (pp indent 
1552            ( "/* Reinitializes IDA state if any synaptic events have occurred */")
1553            ("if (events > 0)")
1554            (#\{)
1555            ("  int status = IDAReInit (B_.sys_, tt, B_.y, B_.yp);")
1556            ("  if (check_flag(&status, \"IDAReInit\", 1)) throw IDASolverFailure (get_name(), status);")
1557            (#\})
1558            ))
1559      ((gsl)
1560        (pp indent 
1561            ( "/* Resets the GSL stepping function if any synaptic events have occurred */")
1562            ("if (events > 0)")
1563            (#\{)
1564            ("  gsl_odeiv2_step_reset (B_.s_);")
1565            (#\})
1566            ))
1567
1568       )
1569    ))
1570
1571
1572
1573(define (spike-handle method defaults vi abstol)
1574  (case method
1575    ((cvode ida) 
1576#<#EOF
1577              {
1578                set_spiketime(Time::ms(tt));
1579                SpikeEvent se;
1580                for (int i = 0; i < N; i++)
1581                {
1582                  S_.y_[i] = Ith(B_.y,i);
1583                }
1584                network()->send(*this, se, lag);
1585                adjust_zero_crossings(B_.y, #{(or abstol 1e-7)});
1586                continue;
1587              }
1588EOF
1589        )
1590
1591    (else
1592     (if (lookup-def 'V_t defaults)
1593         (sprintf 
1594#<<EOF
1595         if ( S_.r_ > 0 )
1596            S_.r_--;
1597         else
1598            if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
1599              {
1600                S_.r_ = V_.RefractoryCounts_;
1601                set_spiketime(Time::step(origin.get_steps()+lag+1));
1602                SpikeEvent se;
1603                network()->send(*this, se, lag);
1604              }
1605EOF
1606        vi vi)  ))
1607
1608    ))
1609
1610
1611
1612(define (output-update
1613         sysname method state-index-map 
1614         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1615         reaction-eq-defs transient-event-defs
1616         i-eqs pool-ions perm-ions synapse-info 
1617         event-external-eq-defs defaults imports
1618         abstol reltol
1619         indent indent+)
1620
1621  (let ((vi (lookup-def 'v state-index-map)))
1622   
1623    (pp indent  (,(sprintf "void ~A::update(Time const & origin, const long_t from, const long_t to)" sysname)
1624                 #\{))
1625
1626    (pp indent+ 
1627#<<EOF
1628    assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
1629    assert(from < to);
1630
1631    double tout;
1632    long_t current_steps = origin.get_steps();
1633
1634    for ( long_t lag = from ; lag < to ; ++lag )
1635      {
1636        double h = B_.step_;   
1637        double tt = 0.0 ;
1638EOF
1639    )
1640
1641    (case method
1642
1643      ((leapfrog)
1644       (begin
1645         (pp indent+ (,(leapfrog-solve vi (spike-handle method defaults vi abstol))))
1646         ))
1647
1648      ((cvode)
1649       (begin
1650         (pp indent+ (,(cvode-solve (spike-handle method defaults vi abstol) )) ,nl)
1651         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
1652         ))
1653
1654      ((ida)
1655       (begin
1656         (pp indent+ (,(ida-solve (spike-handle method defaults vi abstol) )) ,nl)
1657         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
1658         ))
1659
1660      ((gsl)
1661       (begin
1662         (pp indent+ (,(gsl-solve vi (spike-handle method defaults vi abstol))) ,nl)
1663         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
1664         ))
1665      )
1666
1667
1668    (pp (+ 6 indent+)
1669        ("B_.I_stim_ = B_.currents_.get_value(lag);")
1670        ("B_.logger_.record_data(current_steps + lag);"))
1671
1672  (pp indent+  #\})
1673  (pp indent  #\})
1674  ))
1675
1676
1677(define (output-nodes-leapfrog
1678         sysname method state-index-map 
1679         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1680         reaction-eq-defs i-eqs pool-ions perm-ions 
1681         synapse-info defaults abstol reltol maxstep
1682         indent indent+)
1683
1684  (let ((N (length state-index-map))
1685        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
1686        )
1687
1688    (pp indent ,nl  (,#<#EOF
1689#{sysname}::#{sysname}()
1690    : Archiving_Node(), 
1691      P_(), 
1692      S_(P_),
1693      B_(*this)
1694{
1695    recordablesMap_.create();
1696}
1697EOF
1698  ))
1699
1700  (pp indent  ,nl  (,#<#EOF
1701#{sysname}::#{sysname} (const #{sysname}& n)
1702    : Archiving_Node(n), 
1703      P_(n.P_), 
1704      S_(n.S_),
1705      B_(n.B_, *this)
1706{
1707}
1708EOF
1709  ))
1710
1711
1712  (pp indent #<#EOF
1713#{sysname}::~#{sysname} ()
1714{
1715    if ( B_.u != NULL ) free (B_.u);
1716    if ( B_.v != NULL ) free (B_.v);
1717    if ( B_.dvdt != NULL ) free (B_.dvdt);
1718}
1719EOF
1720  )
1721
1722  (pp indent ,nl (,#<#EOF
1723  void #{sysname}::init_node_(const Node& proto)
1724{
1725    const #{sysname}& pr = downcast<#{sysname}>(proto);
1726    P_ = pr.P_;
1727    S_ = State_(P_);
1728}
1729EOF
1730  ))
1731
1732  (pp indent ,nl (,#<#EOF
1733void #{sysname}::init_state_(const Node& proto)
1734{
1735    const #{sysname}& pr = downcast<#{sysname}>(proto);
1736    S_ = State_(pr.P_);
1737}
1738EOF
1739  ))
1740))
1741
1742
1743(define (output-nodes-cvode
1744         sysname method state-index-map 
1745         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1746         reaction-eq-defs i-eqs pool-ions perm-ions 
1747         synapse-info defaults abstol reltol maxstep
1748         indent indent+)
1749
1750  (let (
1751        (N (length state-index-map))
1752        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
1753        )
1754
1755    (pp indent ,nl  (,#<#EOF
1756#{sysname}::#{sysname}()
1757    : Archiving_Node(), 
1758      P_(), 
1759      S_(P_),
1760      B_(*this)
1761{
1762    recordablesMap_.create();
1763}
1764EOF
1765  ))
1766
1767
1768  (pp indent  ,nl  (,#<#EOF
1769#{sysname}::#{sysname}(const #{sysname}& n)
1770    : Archiving_Node(n), 
1771      P_(n.P_), 
1772      S_(n.S_),
1773      B_(n.B_, *this)
1774{
1775}
1776EOF
1777  ))
1778
1779  (pp indent (,#<#EOF
1780#{sysname}::~#{sysname}()
1781{
1782
1783  if ( B_.sys_ )
1784  {
1785    /* Free y vector */
1786    N_VDestroy_Serial(B_.y);
1787
1788    /* Free integrator memory */
1789    if (B_.sys_ != NULL)
1790    {
1791      CVodeFree(&B_.sys_);
1792      B_.sys_ = NULL;
1793    }
1794  }
1795}
1796EOF
1797   ))
1798
1799
1800  (pp indent ,nl (,#<#EOF
1801void #{sysname}::init_node_(const Node& proto)
1802{
1803    const #{sysname}& pr = downcast<#{sysname}>(proto);
1804    P_ = pr.P_;
1805    S_ = State_(P_);
1806}
1807EOF
1808  ))
1809
1810  (pp indent ,nl (,#<#EOF
1811void #{sysname}::init_state_(const Node& proto)
1812{
1813    const #{sysname}& pr = downcast<#{sysname}>(proto);
1814    S_ = State_(pr.P_);
1815}
1816EOF
1817  ))
1818
1819  (pp indent ,nl (,#<#EOF
1820void #{sysname}::init_buffers_()
1821{
1822EOF
1823  ))
1824
1825  (for-each (lambda (psc)
1826              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
1827            pscs)
1828
1829  (pp indent+ (#<<EOF
1830    B_.currents_.clear();           
1831    Archiving_Node::clear_history();
1832
1833    B_.logger_.reset();
1834
1835    B_.step_ = Time::get_resolution().get_ms();
1836    B_.IntegrationStep_ = B_.step_;
1837
1838    B_.I_stim_ = 0.0;
1839EOF
1840))
1841
1842   (pp indent+ (,#<#EOF
1843
1844    int status, N, rootdir;
1845
1846    N = #{N};
1847    // only positive direction (rising edge) of spike events will be detected
1848    rootdir = 1;
1849
1850    /* Creates serial vector of length N */
1851    B_.y = N_VNew_Serial(N);
1852    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw CVodeSolverFailure (get_name(), 0);
1853
1854    for (int i = 0; i < N; i++)
1855    {
1856       Ith(B_.y,i) = S_.y_[i];
1857    }
1858 
1859    /* Calls CVodeCreate to create the solver memory and specify the 
1860     * Backward Differentiation Formula and the use of a Newton iteration */
1861    B_.sys_ = CVodeCreate(CV_BDF, CV_NEWTON);
1862    if (check_flag((void *)B_.sys_, "CVodeCreate", 0)) throw CVodeSolverFailure (get_name(), 0);
1863
1864   /* Calls CVodeInit to initialize the integrator memory and specify the
1865    * right hand side function in y'=f(t,y), the initial time, and
1866    * the initial values. */
1867    status = CVodeInit (B_.sys_, #{sysname}_dynamics, 0.0, B_.y);
1868    if (check_flag(&status, "CVodeInit", 1)) throw CVodeSolverFailure (get_name(), status);
1869
1870EOF
1871    )
1872
1873    ,(if (lookup-def 'V_t defaults) #<#EOF
1874    /* Spike event handler (detects zero-crossing of V-V_t) */
1875    status = CVodeRootInit(B_.sys_, 1, (CVRootFn)#{sysname}_event);
1876    if (check_flag(&status, "CVodeRootInit", 1)) throw CVodeSolverFailure (get_name(), status);
1877
1878    /* Detect only the rising edge of spikes */
1879    status = CVodeSetRootDirection(B_.sys_, &rootdir);
1880    if (check_flag(&status, "CVodeSetRootDirection", 1)) throw CVodeSolverFailure (get_name(), status);
1881EOF
1882     '())
1883
1884     (,#<#EOF
1885    /* Sets the relative and absolute error tolerances of CVode  */
1886    status = CVodeSStolerances (B_.sys_, #{(or abstol 1e-7)}, #{(or reltol 1e-7)});
1887    if (check_flag(&status, "CVodeSStolerances", 1)) throw CVodeSolverFailure (get_name(), status);
1888
1889    /* Turn on CVode stability limit detection (only applicable for order 3 and greater) */
1890    status = CVodeSetStabLimDet (B_.sys_,true);
1891    if (check_flag(&status, "CVodeSetStabLimDet", 1)) throw CVodeSolverFailure (get_name(), status);
1892
1893    /* Sets the maximum order of CVode  */
1894    status = CVodeSetMaxOrd (B_.sys_,5);
1895    if (check_flag(&status, "CVodeSetMaxOrd", 1)) throw CVodeSolverFailure (get_name(), status);
1896
1897    /* Sets maximum step size. */
1898    status = CVodeSetMaxStep (B_.sys_,#{(or maxstep "B_.step_")});
1899    if (check_flag(&status, "CVodeSetMaxStep", 1)) throw CVodeSolverFailure (get_name(), status);
1900
1901    /* Calls CVodeSetUserData to configure the integrator to pass the 
1902     * params structure to the right-hand function */
1903    status = CVodeSetUserData(B_.sys_, reinterpret_cast<void*>(this));
1904    if (check_flag(&status, "CVodeSetUserData", 1)) throw CVodeSolverFailure (get_name(), status);
1905
1906    /* Initializes diagonal linear solver. */
1907    status = CVDiag (B_.sys_);
1908    if (check_flag(&status, "CVDiag", 1)) throw CVodeSolverFailure (get_name(), status);
1909    }
1910
1911EOF
1912    ))
1913
1914  (pp indent (,#<#EOF
1915void #{sysname}::calibrate()
1916{
1917    B_.logger_.init()
1918}
1919
1920EOF
1921   ))
1922 
1923))
1924
1925
1926(define (output-nodes-ida
1927         sysname method state-index-map 
1928         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1929         reaction-eq-defs i-eqs pool-ions perm-ions 
1930         synapse-info defaults abstol reltol maxstep
1931         indent indent+)
1932
1933  (let (
1934        (N (length state-index-map))
1935        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
1936        )
1937
1938    (pp indent ,nl  (,#<#EOF
1939#{sysname}::#{sysname}()
1940    : Archiving_Node(), 
1941      P_(), 
1942      S_(P_),
1943      B_(*this)
1944{
1945    recordablesMap_.create();
1946}
1947EOF
1948  ))
1949
1950
1951  (pp indent  ,nl  (,#<#EOF
1952#{sysname}::#{sysname}(const #{sysname}& n)
1953    : Archiving_Node(n), 
1954      P_(n.P_), 
1955      S_(n.S_),
1956      B_(n.B_, *this)
1957{
1958}
1959EOF
1960  ))
1961
1962  (pp indent (,#<#EOF
1963#{sysname}::~#{sysname}()
1964{
1965
1966  if ( B_.sys_ )
1967  {
1968    /* Free y vector */
1969    N_VDestroy_Serial(B_.y);
1970    N_VDestroy_Serial(B_.yp);
1971
1972    /* Free integrator memory */
1973    if (B_.sys_ != NULL)
1974    {
1975      IDAFree(&B_.sys_);
1976      B_.sys_ = NULL;
1977    }
1978
1979  }
1980}
1981EOF
1982   ))
1983
1984
1985  (pp indent ,nl (,#<#EOF
1986void #{sysname}::init_node_(const Node& proto)
1987{
1988    const #{sysname}& pr = downcast<#{sysname}>(proto);
1989    P_ = pr.P_;
1990    S_ = State_(P_);
1991}
1992EOF
1993  ))
1994
1995  (pp indent ,nl (,#<#EOF
1996void #{sysname}::init_state_(const Node& proto)
1997{
1998    const #{sysname}& pr = downcast<#{sysname}>(proto);
1999    S_ = State_(pr.P_);
2000}
2001EOF
2002  ))
2003
2004  (pp indent ,nl (,#<#EOF
2005void #{sysname}::init_buffers_()
2006{
2007EOF
2008  ))
2009
2010  (for-each (lambda (psc)
2011              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
2012            pscs)
2013
2014  (pp indent+ (#<<EOF
2015    B_.currents_.clear();           
2016    Archiving_Node::clear_history();
2017
2018    B_.logger_.reset();
2019
2020    B_.step_ = Time::get_resolution().get_ms();
2021    B_.IntegrationStep_ = B_.step_;
2022
2023    B_.I_stim_ = 0.0;
2024EOF
2025))
2026
2027   (pp indent+ (,#<#EOF
2028
2029    int status, N, rootdir;
2030
2031    N = #{N};
2032    // only positive direction (rising edge) of spike events will be detected
2033    rootdir = 1;
2034
2035    /* Creates serial vectors of length N */
2036    B_.y = N_VNew_Serial(N);
2037    B_.y1 = N_VNew_Serial(N);
2038    B_.yp = N_VNew_Serial(N);
2039    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw IDASolverFailure (get_name(), 0);
2040
2041    for (int i = 0; i < N; i++)
2042    {
2043       Ith(B_.y,i) = S_.y_[i];
2044    }
2045 
2046    #{sysname}_dynamics (0.0, B_.y, B_.yp, reinterpret_cast<void*>(this));
2047
2048    /* Calls IDACreate to create the solver memory */ 
2049    B_.sys_ = IDACreate();
2050    if (check_flag((void *)B_.sys_, "IDACreate", 0)) throw IDASolverFailure (get_name(), 0);
2051
2052   /* Calls IDAInit to initialize the integrator memory and specify the
2053    * resdual function, the initial time, and the initial values. */
2054    status = IDAInit (B_.sys_, #{sysname}_residual, 0.0, B_.y, B_.yp);
2055
2056    if (check_flag(&status, "IDAInit", 1)) throw IDASolverFailure (get_name(), status);
2057
2058EOF
2059    )
2060
2061    ,(if (lookup-def 'V_t defaults) #<#EOF
2062    /* Spike event handler (detects zero-crossing of V-V_t) */
2063    status = IDARootInit(B_.sys_, 1, (IDARootFn)#{sysname}_event);
2064    if (check_flag(&status, "IDARootInit", 1)) throw IDASolverFailure (get_name(), status);
2065
2066    /* Detect only the rising edge of spikes */
2067    status = IDASetRootDirection(B_.sys_, &rootdir);
2068    if (check_flag(&status, "IDASetRootDirection", 1)) throw IDASolverFailure (get_name(), status);
2069EOF
2070     '())
2071
2072     (,#<#EOF
2073    /* Sets the relative and absolute error tolerances of IDA  */
2074    status = IDASStolerances (B_.sys_, #{(or abstol 1e-7)}, #{(or reltol 1e-7)});
2075    if (check_flag(&status, "IDASStolerances", 1)) throw IDASolverFailure (get_name(), status);
2076
2077    /* Sets the maximum order of IDA  */
2078    status = IDASetMaxOrd (B_.sys_,5);
2079    if (check_flag(&status, "IDASetMaxOrd", 1)) throw IDASolverFailure (get_name(), status);
2080
2081    /* Sets maximum step size. */
2082    status = IDASetMaxStep (B_.sys_,#{(or maxstep "B_.step_")});
2083    if (check_flag(&status, "IDASetMaxStep", 1)) throw IDASolverFailure (get_name(), status);
2084
2085    /* Calls IDASetUserData to configure the integrator to pass the 
2086     * params structure to the right-hand function */
2087    status = IDASetUserData(B_.sys_, reinterpret_cast<void*>(this));
2088    if (check_flag(&status, "IDASetUserData", 1)) throw IDASolverFailure (get_name(), status);
2089
2090    /* Initializes dense linear solver. */
2091    status = IDADense (B_.sys_, N);
2092    if (check_flag(&status, "IDADense", 1)) throw IDASolverFailure (get_name(), status);
2093
2094    status = IDACalcIC(B_.sys_, IDA_Y_INIT, 0.0);
2095    if (check_flag(&status, "IDACalcIC", 1)) throw IDASolverFailure (get_name(), status);
2096
2097    }
2098
2099EOF
2100    ))
2101
2102  (pp indent (,#<#EOF
2103void #{sysname}::calibrate()
2104{
2105    B_.logger_.init()
2106}
2107
2108EOF
2109   ))
2110 
2111))
2112
2113
2114(define (output-nodes-gsl
2115         sysname method state-index-map 
2116         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2117         reaction-eq-defs i-eqs pool-ions perm-ions 
2118         synapse-info defaults abstol reltol maxstep
2119         indent indent+)
2120
2121  (let ((N    (length state-index-map))
2122        (iv   (lookup-def 'v state-index-map))
2123        (pscs (lookup-def 'post-synaptic-conductances synapse-info)))
2124
2125    (pp indent ,nl  (,#<#EOF
2126#{sysname}::#{sysname}()
2127    : Archiving_Node(), 
2128      P_(), 
2129      S_(P_),
2130      B_(*this)
2131{
2132    recordablesMap_.create();
2133}
2134EOF
2135  ))
2136
2137  (pp indent  ,nl  (,#<#EOF
2138#{sysname}::#{sysname}(const #{sysname}& n)
2139    : Archiving_Node(n), 
2140      P_(n.P_), 
2141      S_(n.S_),
2142      B_(n.B_, *this)
2143{
2144}
2145EOF
2146  ))
2147
2148
2149  (pp indent ,#<#EOF
2150#{sysname}::~#{sysname} ()
2151{
2152    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
2153    if ( B_.s_ != NULL) gsl_odeiv2_step_free (B_.s_);
2154    if ( B_.c_ != NULL) gsl_odeiv2_control_free (B_.c_);
2155    if ( B_.e_ != NULL) gsl_odeiv2_evolve_free (B_.e_);
2156    if ( B_.u != NULL) free (B_.u);
2157    if ( B_.jac != NULL) free (B_.jac);
2158}
2159EOF
2160  )
2161
2162  (pp indent ,nl (,#<#EOF
2163  void #{sysname}::init_node_(const Node& proto)
2164{
2165    const #{sysname}& pr = downcast<#{sysname}>(proto);
2166    P_ = pr.P_;
2167    S_ = State_(P_);
2168}
2169EOF
2170  ))
2171
2172  (pp indent ,nl (,#<#EOF
2173void #{sysname}::init_state_(const Node& proto)
2174{
2175    const #{sysname}& pr = downcast<#{sysname}>(proto);
2176    S_ = State_(pr.P_);
2177}
2178EOF
2179  ))
2180
2181  (pp indent ,nl (,#<#EOF
2182void #{sysname}::init_buffers_()
2183{
2184EOF
2185  ))
2186
2187  (for-each (lambda (psc)
2188              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
2189            pscs)
2190
2191  (pp indent+ (#<<EOF
2192    B_.currents_.clear();           
2193    Archiving_Node::clear_history();
2194
2195    B_.logger_.reset();
2196
2197    B_.step_ = Time::get_resolution().get_ms();
2198    B_.IntegrationStep_ = B_.step_;
2199
2200    B_.I_stim_ = 0.0;
2201EOF
2202))
2203
2204    (pp indent+ (,#<#EOF
2205    static const gsl_odeiv2_step_type* T1 = gsl_odeiv2_step_rk2;
2206    B_.N = #{N};
2207 
2208    if ( B_.s_ == 0 )
2209      B_.s_ = gsl_odeiv2_step_alloc (T1, B_.N);
2210    else
2211      gsl_odeiv2_step_reset(B_.s_);
2212   
2213    if ( B_.c_ == 0 ) 
2214      B_.c_ = gsl_odeiv2_control_standard_new (#{(or abstol 1e-7)}, #{(or reltol 1e-7)}, 1.0, 0.0);
2215    else
2216      gsl_odeiv2_control_init(B_.c_, #{(or abstol 1e-7)}, #{(or reltol 1e-7)}, 1.0, 0.0);
2217   
2218    if ( B_.e_ == 0 ) 
2219      B_.e_ = gsl_odeiv2_evolve_alloc(B_.N);
2220    else
2221      gsl_odeiv2_evolve_reset(B_.e_);
2222 
2223    B_.sys_.function  = #{sysname}_dynamics;
2224    B_.sys_.jacobian  = #{sysname}_jacobian;
2225    B_.sys_.dimension = B_.N;
2226    B_.sys_.params    = reinterpret_cast<void*>(this);
2227
2228    B_.u = (double *)malloc(sizeof(double) * B_.N);
2229    assert (B_.u);
2230    B_.jac = (double *)malloc(sizeof(double) * B_.N);
2231    assert (B_.jac);
2232
2233EOF
2234   ))
2235
2236  (pp indent ,nl #\})
2237
2238
2239  (pp indent (,#<#EOF
2240void #{sysname}::calibrate()
2241{
2242    B_.logger_.init()
2243    #{(if iv (sprintf "V_.U_old_ = S_.y_[~A];" iv) "")}
2244    #{(if (lookup-def 't_ref defaults) 
2245          "V_.RefractoryCounts_ = Time(Time::ms(P_.t_ref)).get_steps();" 
2246          "V_.RefractoryCounts_ = 0;")}
2247}
2248EOF
2249))
2250))
2251
2252
2253
2254(define (output-buffers
2255         sysname method state-index-map 
2256         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2257         reaction-eq-defs i-eqs pool-ions perm-ions 
2258         synapse-info defaults
2259         indent indent+)
2260
2261  (let ((N (length state-index-map)))
2262
2263    (case method
2264
2265      ((leapfrog)
2266
2267       (pp indent ,nl
2268           (,#<#EOF
2269#{sysname}::Buffers_::Buffers_(#{sysname}& n)
2270    : logger_(n),
2271      sys_(0),
2272      N(0),
2273      u(0),
2274      v(0),
2275      dvdt(0)
2276{
2277    // Initialization of the remaining members is deferred to
2278    // init_buffers_().
2279}
2280EOF
2281)))
2282
2283   ((cvode ida)
2284
2285    (pp indent ,nl
2286        (,#<#EOF
2287#{sysname}::Buffers_::Buffers_(#{sysname}& n)
2288    : logger_(n),
2289      sys_(0)
2290{
2291    // Initialization of the remaining members is deferred to
2292    // init_buffers_().
2293}
2294EOF
2295        )
2296       (,#<#EOF
2297#{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n)
2298    : logger_(n),
2299      sys_(0)
2300{
2301    // Initialization of the remaining members is deferred to
2302    // init_buffers_().
2303}
2304EOF
2305    )))
2306
2307    ((gsl)
2308
2309     (pp indent ,nl
2310         ( ,#<#EOF
2311#{sysname}::Buffers_::Buffers_(#{sysname}& n)
2312    : logger_(n),
2313      s_(0),
2314      c_(0),
2315      e_(0),
2316      N(0),
2317      u(0),
2318      jac(0)
2319{
2320    // Initialization of the remaining members is deferred to
2321    // init_buffers_().
2322}
2323
2324EOF
2325          )
2326        ( ,#<#EOF
2327#{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n)
2328    : logger_(n),
2329      s_(0),
2330      c_(0),
2331      e_(0),
2332      N(0),
2333      u(0),
2334      jac(0)
2335{
2336    // Initialization of the remaining members is deferred to
2337    // init_buffers_().
2338}
2339EOF
2340    )))
2341   ))
2342)
2343
2344
2345(define (output-accessors+modifiers
2346         sysname imports defaults state-index-map const-defs asgn-eq-defs rate-eq-defs 
2347         reaction-eq-defs i-eqs pool-ions perm-ions constraints
2348         indent indent+)
2349
2350  (let ((c-eqs (lookup-def 'c-eqs constraints))
2351       
2352        (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
2353                                     (list (nest-name n) v)))
2354                      (lookup-def 'c-units constraints)))
2355        )
2356
2357    (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
2358    (pp indent  #\{)
2359   
2360    (for-each
2361     (lambda (def)
2362       (let ((n (first def)))
2363         (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
2364     (append const-defs
2365             defaults))
2366   
2367    (pp indent  #\})
2368   
2369    (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
2370    (pp indent  #\{)
2371   
2372    (for-each
2373     (lambda (def)
2374       (let* ((n (first def))
2375              (nu (lookup-def n c-units))
2376              (scale (and nu (nemo:unit-scale nu)))
2377              )
2378         (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))
2379         (if scale (pp indent+ (,(sprintf "~A = ~A * ~A;" n scale n ))))
2380         ))
2381     (append const-defs
2382             defaults))
2383
2384    (for-each
2385     (lambda (def)
2386       (let ((n (first def))
2387             (b (second def)))
2388         (pp indent+ (,(expr->string/C++ (nest-name n) (nest-name b))))))
2389     defaults)
2390   
2391    (pp indent  #\})
2392   
2393    (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
2394    (pp indent  #\{)
2395   
2396    (for-each
2397     (lambda (def)
2398       (let* (
2399              (n (first def))
2400              (i (second def))
2401              )
2402         (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))
2403         ))
2404     state-index-map)
2405   
2406    (let ((vi (lookup-def 'v state-index-map)))
2407      (if vi
2408          (pp indent+ (,(sprintf "def<double_t>(d, names::V_m, y_[~A]);"  vi) ))
2409          ))
2410   
2411    (pp indent  #\})
2412   
2413    (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
2414    (pp indent  #\{)
2415   
2416    (for-each
2417     (lambda (def)
2418       (let* ((n     (first def)) 
2419              (i     (second def))
2420              (nu    (lookup-def n c-units))
2421              (scale (and nu (nemo:unit-scale nu)))
2422              )
2423         (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))
2424         (if scale (pp indent+ (,(sprintf "y_[~A] = ~A * y_[~A];" i scale i ))))
2425         ))
2426     state-index-map)
2427   
2428    (let ((vi (lookup-def 'v state-index-map)))
2429      (if vi
2430          (pp indent+ (,(sprintf "updateValue<double_t>(d, names::V_m, y_[~A]);"  vi) ))
2431          ))
2432   
2433    (pp indent  #\})
2434   
2435    ))
2436
2437
2438(define (output-init sysname method state-index-map steady-state-index-map 
2439                     imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2440                     reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints 
2441                     indent indent+)
2442
2443 
2444  (let* ((N (length state-index-map))
2445
2446         (c-eqs (lookup-def 'c-eqs constraints))
2447
2448         (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
2449                                     (list (nest-name n) v)))
2450                       (lookup-def 'c-units constraints)))
2451
2452         (i-eqs 
2453          (map
2454           (lambda (def) (list (first def) 
2455                               (add-params-to-fncall 
2456                                (canonicalize-expr/C++ (second def)) builtin-fns)))
2457           i-eqs))
2458
2459         (init-eqs 
2460          (append
2461             
2462           (map (lambda (def)
2463                  (let ((n (first def))
2464                        (b (second def)))
2465                    (list (nest-name n) (nest-name b))))
2466                external-eq-defs)
2467           
2468           asgn-eq-defs
2469           init-eq-defs
2470
2471           (map (lambda (pool-ion)
2472                  (let ((n (pool-ion-in pool-ion))
2473                        (b (pool-ion-inq pool-ion)))
2474                    (list n b)))
2475                pool-ions)
2476
2477           (map (lambda (pool-ion)
2478                  (let ((n (pool-ion-out pool-ion))
2479                        (b (pool-ion-outq pool-ion)))
2480                    (list n b)))
2481                pool-ions)
2482           ))
2483
2484         (init-dag 
2485          (map (lambda (def)
2486                 (cons (first def) (enum-freevars (second def) '() '())))
2487               init-eqs))
2488
2489         (init-order
2490          (reverse
2491           (topological-sort init-dag 
2492                             (lambda (x y) (string=? (->string x) (->string y))))))
2493
2494         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
2495
2496         )
2497   
2498    (pp indent ,nl
2499        (,(sprintf "~A::State_::State_" sysname) (const Parameters_& p))
2500        ( #\{))
2501
2502    (pp indent+
2503        (double ,(slp ", " (delete-duplicates
2504                            (map ->string 
2505                                 (filter (lambda (x) (not (member x builtin-consts)))
2506                                         (append
2507                                          init-locals
2508                                          init-order
2509                                          (map first external-eq-defs)
2510                                          (map pool-ion-in pool-ions)
2511                                          (map pool-ion-out pool-ions)
2512                                          (map first i-eqs)
2513                                          (map first state-index-map) 
2514                                          (map first const-defs)
2515                                          (map first reaction-eq-defs)
2516                                          )))
2517                                 string=?))
2518                      ";")
2519        ("const Parameters_ *params = &p;")
2520        )
2521
2522    (case method
2523      ((gsl leapfrog) 
2524       (pp indent+ "r_ = 0;" ,nl)))
2525
2526    (pp indent+ ,(sprintf "memset(y_,0,~A*sizeof(double));" N))
2527
2528    (for-each
2529     (lambda (x) 
2530       (let* ((n  (first x)))
2531         (pp indent+ ,(expr->string/C++ (sprintf "p.~A" n) n))))
2532     const-defs)
2533
2534    (let ((vi (lookup-def 'v state-index-map))
2535          (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
2536      (if (and vi vrest) 
2537          (pp indent+ ,(expr->string/C++  vrest 'v ))))
2538
2539    (for-each (lambda (n)
2540                (let ((b  (lookup-def n init-eqs)))
2541                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
2542              init-order)
2543   
2544    (for-each (lambda (def)
2545                (let* ((n  (first def)) 
2546                       (ni (lookup-def n state-index-map)))
2547                  (if ni (pp indent+ ,(expr->string/C++ n  (sprintf "y_[~A]" ni))))))
2548              init-eq-defs)
2549
2550    (for-each
2551     (lambda (eq) 
2552       (match-let (((op left right)  eq))
2553                  (pp indent+ ,(sprintf "if (!((~A) ~A (~A))) { throw BadProperty (~S); }; " 
2554                                        (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ left))) op
2555                                        (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ right)))
2556                                        (sprintf "Constraint ~A is not satisfied." eq)
2557                                        ))
2558                  ))
2559     c-eqs)
2560
2561
2562    (if (not (null? steady-state-index-map))
2563
2564        (case method
2565
2566          ((cvode ida)
2567
2568           (let ((N (length steady-state-index-map))
2569                 (ssvect (gensym 'ssvect)))
2570
2571             (pp indent+
2572                 ,(sprintf "N_Vector ~A;" ssvect)
2573                 ,(expr->string/C++ (sprintf "N_VNew_Serial(~A)" N) ssvect)
2574                 (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ,N ,ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
2575             
2576             (for-each
2577              (lambda (def)
2578                (let* ((n (first def)) 
2579                       (si (lookup-def n steady-state-index-map))
2580                       (ni (lookup-def n state-index-map)))
2581                  (if si
2582                        (pp indent+
2583                            ,(expr->string/C++ (ith ssvect si)
2584                                               (sprintf "y_[~A]" ni))
2585                            ,(expr->string/C++ (ith ssvect si)
2586                                               n))
2587                        )
2588                  ))
2589              rate-eq-defs)
2590             
2591             (pp indent+ (,(sprintf "N_VDestroy_Serial (~A);" ssvect)))
2592
2593             ))
2594
2595          (else
2596
2597           (let ((N (length steady-state-index-map))
2598                 (ssvect (gensym 'ssvect))
2599                 )
2600
2601             (pp indent+
2602                 ,(sprintf "gsl_vector *~A;" ssvect)
2603                 ,(expr->string/C++ (sprintf "gsl_vector_alloc (~A)" N) ssvect)
2604                 (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ,N ,ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
2605             
2606             (for-each
2607              (lambda (def)
2608                (let* ((n (first def)) 
2609                       (si (lookup-def n steady-state-index-map))
2610                       (ni (lookup-def n state-index-map)))
2611                  (if si
2612                        (pp indent+
2613                            ,(expr->string/C++ (sprintf "gsl_vector_get (~A,~A)" ssvect si)
2614                                               (sprintf "y_[~A]" ni))
2615                            ,(expr->string/C++ (sprintf "gsl_vector_get (~A,~A)" ssvect si)
2616                                               n))
2617                        )
2618                  ))
2619              rate-eq-defs)
2620             
2621             (pp indent+ (,(sprintf "gsl_vector_free (~A);" ssvect)))
2622
2623             ))
2624
2625          ))
2626   
2627    (for-each
2628     (lambda (def)
2629       (let ((n (first def)) (b (second def)))
2630         (if (not (lookup-def n init-eq-defs))
2631             (pp indent+ ,(expr->string/C++ b n)))))
2632     reaction-eq-defs)
2633   
2634    (for-each
2635     (lambda (def) 
2636       (pp indent+ ,(expr->string/C++ (second def) (first def))))
2637     i-eqs)
2638
2639    (let ((vi (lookup-def 'v state-index-map)))
2640      (if vi (pp indent+ ,(expr->string/C++ 'v (sprintf "y_[~A]" vi)))))
2641     
2642    (pp indent  #\})
2643
2644    (pp indent ,nl
2645        (,(sprintf "~A::State_::State_" sysname) (const State_& s)  )
2646        ( #\{)
2647       
2648        ,(case method
2649          ((gsl leapfrog) 
2650           `("r_ = s.r_;"))
2651          (else ""))
2652       
2653        (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N))
2654        ( #\})
2655        ,nl
2656        )
2657
2658    (pp indent ,nl
2659        (,(sprintf "~A::State_& ~A::State_::operator=(const State_& s)" sysname sysname))
2660        ( #\{)
2661        )
2662
2663    (case method
2664      ((gsl leapfrog) 
2665       (pp indent+ "r_ = s.r_;" ,nl)))
2666
2667    (pp indent+ ( ,#<#EOF
2668     assert(this != &s)
2669     for ( size_t i = 0 ; i < #{N} ; ++i )
2670       y_[i] = s.y_[i];
2671EOF
2672))
2673
2674        (pp indent+ "return *this;")
2675
2676        (pp indent  #\})
2677))
2678
2679
2680(define (output-parameters sysname imports const-defs external-eq-defs defaults indent indent+)
2681
2682  (let* ((parameter-defs
2683          (map
2684           (lambda (def) (list (first def) 
2685                               (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
2686           const-defs))
2687
2688         (parameter-locals  (find-locals (map second parameter-defs)))
2689
2690         (const-exprs
2691            (map (lambda (def)
2692                   (let* ((n  (first def)) (b (second def)))
2693                     (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
2694                 const-defs) )
2695             
2696          (default-exprs 
2697            (map (lambda (def)
2698                   (let ((n (first def))
2699                         (b (second def)))
2700                     (expr->string/C++ (nest-name b) (nest-name n))))
2701                 defaults))
2702           )
2703
2704      (pp indent ,nl (,(sprintf "~A::Parameters_::Parameters_" sysname) () ":"))
2705     
2706      (if (not (null? parameter-locals)) 
2707          (pp indent+ (double ,(slp ", " parameter-locals) ";")))
2708
2709      (pp indent+ ,(slp ",\n" const-exprs))
2710     
2711      (pp indent  #\{ ,(slp "\n" default-exprs) #\})
2712     
2713      ))
2714
2715
2716(define (output-recordables
2717         sysname state-index-map 
2718         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2719         reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
2720 
2721  (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
2722  (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
2723  (pp indent #\{)
2724
2725  (for-each
2726   (lambda (def)
2727     (let ((n (first def)) (i (second def)))
2728       (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);" 
2729                              (->string n) sysname sysname (string-upcase (->string n)))))
2730       ))
2731   state-index-map)
2732
2733  (if (lookup-def 'v state-index-map)
2734      (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);" 
2735                             sysname sysname (string-upcase (->string 'v))))))
2736
2737  (pp indent #\})
2738  )
2739
2740
2741
2742(define (output-dy sysname method imports const-defs state-index-map 
2743                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
2744                   pool-ions i-eqs v-eq constraints indent indent+)
2745
2746  (let* (
2747         (c-units (map (lambda (x) 
2748                         (let ((n (first x)) (v (second x)))
2749                           (list (nest-name n) v)))
2750                       (lookup-def 'c-units constraints)))
2751
2752         (i-eqs 
2753          (map
2754           (lambda (def) (list (first def) 
2755                               (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
2756           i-eqs))
2757
2758         (v-eq
2759          (and v-eq
2760               (list (first v-eq) 
2761                     (add-params-to-fncall (canonicalize-expr/C++ (second v-eq)) builtin-fns))))
2762
2763         (eqs 
2764          (append
2765
2766           external-eq-defs
2767           asgn-eq-defs
2768           reaction-eq-defs
2769                 
2770           (map (lambda (pool-ion)
2771                  (let ((n (pool-ion-in pool-ion))
2772                        (b (pool-ion-inq pool-ion)))
2773                    (list n b)))
2774                pool-ions)
2775
2776           (map (lambda (pool-ion)
2777                  (let ((n (pool-ion-out pool-ion))
2778                        (b (pool-ion-outq pool-ion)))
2779                    (list n b)))
2780                pool-ions)
2781
2782           i-eqs
2783          ))
2784         
2785         (eq-dag 
2786          (map (lambda (def)
2787                 (cons (first def) (enum-freevars (second def) '() '())))
2788               eqs))
2789
2790         (eq-order
2791          (reverse
2792           (topological-sort eq-dag 
2793                             (lambda (x y) (string=? (->string x) (->string y))))))
2794
2795         (eq-locals  (find-locals 
2796                      (map second
2797                           (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
2798                               (append i-eqs rate-eq-defs eqs)))))
2799         )
2800
2801
2802    (case method
2803      ((leapfrog)
2804       (pp indent ,nl (
2805                       ,(sprintf "extern \"C\" int ~A_dynamics" sysname)
2806                       (,(slp ", " `("double t"
2807                                     "const double y[]"
2808                                     "double f[]"
2809                                     "void* pnode"
2810                                     )))
2811                       #\{
2812                       )))
2813      ((cvode ida)
2814       (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname)
2815                       (,(slp ", " `("double t"
2816                                     "N_Vector y"
2817                                     "N_Vector f"
2818                                     "void* pnode"
2819                                     )))
2820                       #\{
2821                       )))
2822      ((gsl)
2823       (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname)
2824                       (,(slp ", " `("double t"
2825                                     "const double y[]"
2826                                     "double f[]"
2827                                     "void* pnode"
2828                                     )))
2829                       #\{
2830                       )))
2831      )
2832
2833
2834
2835    (pp indent+ (double ,(slp ", " (delete-duplicates 
2836                                    (map (compose ->string nest-name)
2837                                         (filter (lambda (x) 
2838                                                   (not (member x builtin-consts)))
2839                                                 (append
2840                                                  eq-locals
2841                                                  eq-order
2842                                                  (map first i-eqs)
2843                                                  (map first external-eq-defs)
2844                                                  (map first state-index-map)
2845                                                  (map first const-defs)
2846                                                  )))
2847                                    string=?))
2848                        ";"))
2849
2850    (pp indent+ ,nl
2851
2852        ("// S is shorthand for the type that describes the model state ")
2853        (typedef ,(sprintf "~A::State_" sysname) "S;")
2854       
2855        ,nl
2856       
2857        ("// cast the node ptr to an object of the proper type")
2858        ("assert(pnode);")
2859        ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2860        (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2861       
2862        ,nl
2863
2864        ("// params is a reference to the model parameters ")
2865        (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
2866        ("params = &(node.P_);")
2867       
2868        ,nl
2869
2870        ("// state is a reference to the model state ")
2871        (struct ,(sprintf  "~A::State_" sysname) "*state;")
2872        ("state = &(vnode.S_);")
2873       
2874        ,nl
2875
2876        ("// y[] must be the state vector supplied by the integrator, ")
2877        ("// not the state vector in the node, node.S_.y[]. ")
2878
2879        ,nl
2880
2881        )
2882
2883    (for-each (lambda (def)
2884                (let ((n (first def)) )
2885                  (pp indent+
2886                      ,(expr->string/C++ (sprintf "params->~A" n) n))))
2887              const-defs)
2888
2889    (let ((vi (lookup-def 'v state-index-map)))
2890      (if vi (pp indent+ ,(expr->string/C++ (ith 'y vi) 'v))))
2891
2892    (for-each (lambda (def)
2893                (let* ((n (first def))
2894                       (ni (lookup-def n state-index-map)))
2895                  (pp indent+ ,(expr->string/C++ (ith 'y ni) (nest-name n)))))
2896              rate-eq-defs)
2897
2898    (for-each (lambda (n)
2899                (let ((b  (lookup-def n eqs)))
2900                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
2901              eq-order)
2902
2903    (for-each (lambda (def)
2904                (let* ((n (first def))
2905                       (b (second def))
2906                       (fv (ith 'f (lookup-def n state-index-map)))
2907                       )
2908                  (pp indent+ ,(s+ "// rate equation " n)
2909                      ,(expr->string/C++ b fv))))
2910              rate-eq-defs)
2911   
2912    (if v-eq
2913        (let ((vi (lookup-def 'v state-index-map)))
2914          (if vi
2915              (pp indent+ ,(expr->string/C++ (second v-eq) (ith 'f vi)))
2916              )))
2917
2918    (case method
2919      ((leapfrog)
2920        (pp indent+ ,nl ("return 0;")))
2921      ((cvode ida)
2922        (pp indent+ ,nl ("return 0;")))
2923      ((gsl)
2924        (pp indent+ ,nl ("return GSL_SUCCESS;")))
2925      )
2926
2927    (pp indent  #\})
2928
2929    ))
2930
2931
2932(define (output-residual sysname method imports const-defs state-index-map 
2933                         external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
2934                         pool-ions i-eqs v-eq constraints indent indent+)
2935  (case method 
2936      ((ida)
2937       (begin
2938         (pp indent ,nl (,(sprintf "extern \"C\" int ~A_residual" sysname)
2939                         (,(slp ", " `("double t"
2940                                       "N_Vector y"
2941                                       "N_Vector yp"
2942                                       "N_Vector f"
2943                                       "void* pnode"
2944                                       )))
2945                         #\{
2946                         ))
2947
2948         (pp indent+ ,nl
2949             (int status #\; )
2950
2951             ,nl
2952
2953             ("// cast the node ptr to an object of the proper type")
2954             ("assert(pnode);")
2955             ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2956             (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
2957             ("N_Vector y1 = vnode.B_.y1;")
2958             
2959             ,nl
2960
2961             (,(sprintf "status = ~A_dynamics (t, y, y1, pnode);" sysname))
2962
2963             ,nl)
2964             
2965         (for-each (lambda (def)
2966                     (let* ((n (first def))
2967                            (i (lookup-def n state-index-map))
2968                            (fv (ith 'f i))
2969                            (y1v (ith 'y1 i))
2970                            (ypv (ith 'yp i))
2971                            )
2972                       (pp indent+
2973                           ,(expr->string/C++ `(- ,y1v ,ypv) fv))))
2974                   rate-eq-defs)
2975
2976         (if v-eq
2977             (let ((vi (lookup-def 'v state-index-map)))
2978               (if vi
2979                   (let ((v1 (ith 'y1 vi))
2980                         (vp (ith 'yp vi)))
2981                     (pp indent+ ,(expr->string/C++ `(- ,v1 ,vp) (ith 'f vi))))
2982                   )))
2983
2984
2985         (pp indent
2986             ("return status; ")
2987             #\})
2988         ))
2989      ))
2990
2991
2992(define (output-steadystate sysname method state-index-map steady-state-index-map 
2993                             imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2994                             reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints 
2995                             indent indent+)
2996
2997 
2998  (let* (
2999         (N (length steady-state-index-map))
3000
3001         (c-eqs (lookup-def 'c-eqs constraints))
3002
3003         (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
3004                                     (list (nest-name n) v)))
3005                       (lookup-def 'c-units constraints)))
3006
3007         (i-eqs 
3008          (map
3009           (lambda (def) (list (first def) 
3010                               (add-params-to-fncall 
3011                                (canonicalize-expr/C++ (second def)) builtin-fns)))
3012           i-eqs))
3013
3014         (init-eqs 
3015          (append
3016             
3017           (map (lambda (def)
3018                  (let ((n (first def))
3019                        (b (second def)))
3020                    (list (nest-name n) (nest-name b))))
3021                external-eq-defs)
3022           
3023           asgn-eq-defs
3024           init-eq-defs
3025
3026           (map (lambda (pool-ion)
3027                  (let ((n (pool-ion-in pool-ion))
3028                        (b (pool-ion-inq pool-ion)))
3029                    (list n b)))
3030                pool-ions)
3031
3032           (map (lambda (pool-ion)
3033                  (let ((n (pool-ion-out pool-ion))
3034                        (b (pool-ion-outq pool-ion)))
3035                    (list n b)))
3036                pool-ions)
3037           ))
3038
3039         (init-dag 
3040          (map (lambda (def)
3041                 (cons (first def) (enum-freevars (second def) '() '())))
3042               init-eqs))
3043
3044         (init-order
3045          (reverse
3046           (topological-sort init-dag 
3047                             (lambda (x y) (string=? (->string x) (->string y))))))
3048
3049         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
3050
3051         )
3052
3053    (case method
3054      ((ida cvode)
3055       (pp indent ,nl
3056           (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(N_Vector u, N_Vector f, void* pnode)" ,nl)
3057           ( #\{)))
3058      (else
3059       (pp indent ,nl
3060           (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(const gsl_vector *u, void *pnode, gsl_vector *f)" ,nl)
3061           ( #\{)))
3062      )
3063
3064    (pp indent+
3065        (double ,(slp ", " (delete-duplicates
3066                            (map ->string 
3067                                 (filter (lambda (x) (not (member x builtin-consts)))
3068                                         (append
3069                                          init-locals
3070                                          init-order
3071                                          (map first external-eq-defs)
3072                                          (map pool-ion-in pool-ions)
3073                                          (map pool-ion-out pool-ions)
3074                                          (map first i-eqs)
3075                                          (map first steady-state-index-map) 
3076                                          (map first const-defs)
3077                                          (map first reaction-eq-defs)
3078                                          )))
3079                                 string=?))
3080                      ";")
3081        )
3082
3083    (pp indent+ ,nl
3084
3085        ("// params is a reference to the model parameters ")
3086        (const struct ,(sprintf "~A::Parameters_" sysname) 
3087               "*params = "
3088               "(" ,(sprintf "struct ~A::Parameters_ *" sysname) ")" "pnode;")
3089       
3090        ,nl
3091        )
3092
3093    (for-each
3094     (lambda (x) 
3095       (let* ((n  (first x)))
3096         (pp indent+ ,(expr->string/C++ (sprintf "params->~A" n) n))))
3097     const-defs)
3098
3099    (case method
3100      ((ida cvode)
3101       (for-each
3102        (lambda (def)
3103          (let* ((n   (first def)) 
3104                 (ni  (lookup-def n steady-state-index-map)))
3105            (if ni (pp indent+ ,(expr->string/C++ (ith 'u ni) n)))
3106            ))
3107        rate-eq-defs))
3108      (else
3109       (for-each
3110        (lambda (def)
3111          (let* ((n   (first def)) 
3112                 (ni  (lookup-def n steady-state-index-map)))
3113            (if ni (pp indent+ ,(expr->string/C++ (sprintf "gsl_vector_get (u, ~A)" ni) n)))
3114            ))
3115        rate-eq-defs))
3116      )
3117       
3118
3119    (let ((vi (lookup-def 'v state-index-map)))
3120      (if vi (pp indent+ ,(expr->string/C++ 0. 'v))))
3121   
3122    (for-each (lambda (def) 
3123                (pp indent+ ,(expr->string/C++ 0. (first def))))
3124              i-eqs)
3125
3126    (for-each (lambda (n)
3127                (let ((b  (lookup-def n init-eqs)))
3128                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
3129              init-order)
3130
3131    (case method
3132      ((ida cvode)
3133       (for-each
3134        (lambda (def)
3135          (let* ((n   (first def)) 
3136                 (ni  (lookup-def n steady-state-index-map))
3137                 (b   (second def))
3138                 (lb  (find-locals (list b))))
3139            (if ni 
3140                (begin
3141                  (if (not (null? lb))
3142                      (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
3143                  (pp indent+
3144                      (,(expr->string/C++ b (ith 'f ni))))))
3145            ))
3146        rate-eq-defs))
3147      (else
3148       (for-each
3149        (lambda (def)
3150          (let* ((n   (first def)) 
3151                 (ni  (lookup-def n steady-state-index-map))
3152                 (b   (second def))
3153                 (lb  (find-locals (list b))))
3154            (if ni 
3155                (let ((tmp (gensym 't)))
3156                  (pp indent+ (double ,tmp #\;))
3157                  (if (not (null? lb))
3158                      (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
3159                  (pp indent+
3160                      ,(expr->string/C++ b tmp)
3161                      (,(sprintf "gsl_vector_set (f,~A,~A);" ni tmp)))
3162                  ))
3163            ))
3164        rate-eq-defs))
3165      )
3166
3167    (pp indent ,nl
3168        ( "return 0;" )
3169        ( #\}))
3170
3171    ))
3172
3173
3174(define (nemo:nest-translator sys . rest)
3175
3176  (define (cid x)  (second x))
3177  (define (cn x)   (first x))
3178
3179  (let-optionals rest ((dirname ".")  (method #f) (abstol #f) (reltol #f) (maxstep #f))
3180
3181    (let ((method (or method 'gsl)))
3182
3183      (if (not (member method '(gsl cvode ida leapfrog)))
3184          (nemo:error 'nemo:nest-translator ": unknown method " method))
3185
3186  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
3187    (let ((imports  ((dis 'imports)  sys))
3188          (exports  ((dis 'exports)  sys)))
3189      (let* ((indent      0)
3190             (indent+     (+ 2 indent ))
3191
3192             (sysname     (nest-name ((dis 'sysname) sys)))
3193             (prefix      (->string sysname))
3194             (deps*       ((dis 'depgraph*)  sys))
3195             (consts      ((dis 'consts)     sys))
3196             (asgns       ((dis 'asgns)      sys))
3197             (states      ((dis 'states)     sys))
3198             (reactions   ((dis 'reactions)  sys))
3199             (defuns      ((dis 'defuns)     sys))
3200             (components  ((dis 'components) sys))
3201             
3202             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
3203             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
3204
3205             (const-defs       (filter-map
3206                                (lambda (nv)
3207                                  (and (not (member (first nv) builtin-consts))
3208                                       (let ((v1 (canonicalize-expr/C++ (second nv))))
3209                                         (list (nest-name (first nv)) v1))))
3210                                consts))
3211             
3212             (defaults             (nemo:defaults-query sys))
3213
3214             (geometry             (nemo:geometry-query sys))
3215
3216             (gate-complex-info    (nemo:gate-complex-query sys))
3217             (perm-ions       (map (match-lambda ((comp i e erev val) `(,comp ,(nest-name i) ,(nest-name e) ,erev)))
3218                                   (lookup-def 'perm-ions gate-complex-info)))
3219             (acc-ions        (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out)))))
3220                                   (lookup-def 'acc-ions gate-complex-info)))
3221             (epools          (lookup-def 'pool-ions gate-complex-info))
3222             (pool-ions       (pool-ion-name-map nest-name  epools))
3223
3224             (comprc         (any (match-lambda ((name 'membrane-tau id) (list name id)) (else #f)) components))
3225             (compcap        (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
3226             (mrc            (or (and comprc (car ((dis 'component-exports) sys (cid comprc))))
3227                                 (lookup-def 'membrane-tau defaults)
3228                                 (lookup-def 'tau_m defaults)
3229                                 (and compcap (car ((dis 'component-exports) sys (cid compcap))))
3230                                 (lookup-def 'membrane-capacitance defaults)
3231                                 (lookup-def 'C_m defaults)
3232                                 ))
3233
3234             (soma-geometry  (lookup-def 'soma geometry))
3235             (marea          (and soma-geometry (third soma-geometry)))
3236
3237             (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
3238             (synapse-info         (nemo:post-synaptic-conductance-query sys))
3239
3240
3241             (pscs           (lookup-def 'post-synaptic-conductances synapse-info))
3242
3243             (i-syns         (lookup-def 'i-synapses synapse-info))
3244               
3245             (i-gates        (lookup-def 'i-gates gate-complex-info))
3246
3247             (i-defs         (nemo:ionic-current-definitions
3248                              gate-complexes i-gates i-syns pscs marea
3249                              (lambda (x) (state-power sys x))
3250                              (lambda (x) ((dis 'component-exports) sys x))
3251                              (lambda (x) ((dis 'component-subcomps) sys x))
3252                              nest-name rhsexpr/C++ canonicalize-expr/C++
3253                              builtin-fns))
3254
3255             (i-eqs          (lookup-def 'i-eqs i-defs))
3256             (i-names        (lookup-def 'i-names i-defs))
3257
3258             (constraints    (nemo:constraint-definitions 
3259                              gate-complexes i-gates i-syns pscs marea imports
3260                              (lambda (x) (state-power sys x))
3261                              (lambda (x) (quantity-unit sys x))
3262                              (lambda (x) ((dis 'component-exports) sys x))
3263                              (lambda (x) ((dis 'component-subcomps) sys x))
3264                              nest-name))
3265
3266             (external-eq-defs   (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++
3267                                                        namespace-filter: (lambda (x) (not (equal? x 'event)))))
3268
3269             (event-external-eq-defs   (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++
3270                                                              namespace-filter: (lambda (x) (equal? x 'event))))
3271
3272             (asgn-eq-defs       (poset->asgn-eq-defs* poset sys nest-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
3273             
3274             (rate-eq-defs       (reverse (poset->rate-eq-defs* poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)))
3275             
3276             (reaction-eq-defs   (poset->reaction-eq-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++))
3277
3278             (transient-event-defs  (poset->transient-event-defs poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)) 
3279             
3280             (init-eq-defs       (poset->init-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
3281             
3282             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
3283                                      (poset->state-conserve-eq-defs poset sys nest-name nest-state-name)))
3284             
3285             (imports-sans-v (filter (lambda (x) (not (equal? 'v (first x)))) imports))
3286
3287             (v-eq    (and (not (null? i-names))
3288                           (let ((istim "(node.B_.I_stim_)" )) 
3289                             (cond
3290
3291                              ((and mrc marea)
3292                               (list 'v (rhsexpr/C++ 
3293                                         `(/ (+ (* ,istim (/ 100. ,marea)) 
3294                                                (* -1e3 ,(sum i-names))) ,mrc))))
3295                              (marea
3296                               (list 'v (rhsexpr/C++ 
3297                                         `(+ (* ,istim (/ 100. ,marea))
3298                                             (* -1e-3 ,(sum i-names))))))
3299                              (mrc
3300                               (list 'v (rhsexpr/C++ `(/ (+ ,istim (* -1e-3 ,(sum i-names))) ,mrc))))
3301                             
3302                              (else
3303                               (list 'v (rhsexpr/C++ `(+ ,istim (* -1e-3 ,(sum i-names))))))
3304                              ))
3305                           ))
3306
3307             
3308             (state-index-map  (let ((acc (fold (lambda (def ax)
3309                                                  (let ((st-name (first def)))
3310                                                    (list (+ 1 (first ax)) 
3311                                                          (cons `(,st-name ,(first ax)) (second ax)))))
3312                                                (list 0 (list)) 
3313                                                (if (lookup-def 'v imports)
3314                                                    (cons (list 'v) rate-eq-defs)
3315                                                    rate-eq-defs)
3316                                                )))
3317                                 (second acc)))
3318             
3319             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
3320                                                         (let ((st-name (first def)))
3321                                                           (if (not (alist-ref st-name init-eq-defs))
3322                                                               (list (+ 1 (first ax)) 
3323                                                                     (cons `(,st-name ,(first ax)) (second ax)))
3324                                                               ax)))
3325                                                       (list 0 (list)) 
3326                                                       rate-eq-defs)))
3327                                        (second acc)))
3328             
3329             (dfenv 
3330              (map (lambda (x) (let ((n (first x))) (list n (nest-name (sprintf "d_~A" n )))))
3331                   defuns))
3332
3333             )
3334       
3335       
3336       
3337        (for-each
3338         (lambda (a)
3339           (let ((acc-ion   (car a)))
3340             (if (assoc acc-ion perm-ions)
3341                 (nemo:error 'nemo:nest-translator 
3342                             ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
3343         acc-ions)
3344
3345
3346        (let ((cpp-output  (open-output-file (make-output-fname dirname prefix ".cpp")))
3347              (hpp-output  (open-output-file (make-output-fname dirname prefix ".h"))))
3348         
3349          ;; include files and other prelude definitions
3350          (with-output-to-port cpp-output
3351            (lambda ()
3352              (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
3353              (output-prelude sysname indent)
3354              ))
3355         
3356          (with-output-to-port cpp-output
3357            (lambda () (pp indent ,nl "namespace nest {" ,nl)))
3358         
3359          (case method
3360            ((cvode ida) 
3361             (with-output-to-port cpp-output
3362               (lambda ()
3363                 (pp indent ,sundials-prelude)))))
3364
3365
3366          (if (not (null? steady-state-index-map))
3367              (with-output-to-port cpp-output
3368                (lambda ()
3369                  (pp indent ,(fsolve-prelude method)))
3370                ))
3371
3372          ;; user-defined functions
3373          (let ((define-fn  (make-define-fn sysname))
3374                (define-fn-header (make-define-fn-header sysname)))
3375           
3376            (for-each (lambda (fndef) 
3377                        (if (not (member (car fndef) builtin-fns))
3378                            (with-output-to-port cpp-output
3379                              (lambda ()
3380                                (apply define-fn-header (cons indent fndef))
3381                                (pp indent ,nl)))
3382                            ))
3383                      defuns)
3384
3385            (for-each (lambda (fndef) 
3386                        (if (not (member (car fndef) builtin-fns))
3387                            (with-output-to-port cpp-output
3388                              (lambda ()
3389                                (apply define-fn (cons indent fndef))
3390                                (pp indent ,nl)))
3391                            ))
3392                      defuns)
3393            )
3394         
3395         
3396          ;; derivative function
3397          (with-output-to-port cpp-output
3398            (lambda ()
3399              (output-dy sysname method imports-sans-v const-defs state-index-map 
3400                         external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3401                         pool-ions i-eqs v-eq constraints
3402                         indent indent+)
3403              ))
3404
3405          ;; residual function (used by IDA only)
3406          (with-output-to-port cpp-output
3407            (lambda ()
3408              (output-residual sysname method imports-sans-v const-defs state-index-map 
3409                               external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3410                               pool-ions i-eqs v-eq constraints
3411                               indent indent+)
3412              ))
3413
3414          ;; steady state function (used by CVODE and IDA)
3415          (with-output-to-port cpp-output
3416            (lambda ()
3417              (output-steadystate sysname method state-index-map steady-state-index-map 
3418                                  imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3419                                  reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+)
3420              (pp indent ,nl)
3421              ))
3422
3423          ;; Jacobian function
3424          (with-output-to-port cpp-output
3425            (lambda ()
3426              (output-jac  sysname method state-index-map pool-ions i-eqs v-eq indent indent+)
3427              ))
3428
3429          ;;  system recordables
3430          (with-output-to-port cpp-output
3431            (lambda ()
3432              (output-recordables
3433               sysname state-index-map 
3434               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3435               reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
3436              (pp indent ,nl)
3437              ))
3438         
3439          ;; system parameters
3440          (with-output-to-port cpp-output
3441            (lambda ()
3442              (output-parameters sysname imports-sans-v const-defs 
3443                                 external-eq-defs defaults indent indent+)
3444              ))
3445         
3446          ;; initial values function
3447          (with-output-to-port cpp-output
3448            (lambda ()
3449              (output-init sysname method state-index-map steady-state-index-map 
3450                           imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3451                           reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+)
3452              (pp indent ,nl)
3453              ))
3454         
3455         
3456          ;; accessors and modifiers for states and parameters
3457          (with-output-to-port cpp-output
3458            (lambda ()
3459              (output-accessors+modifiers
3460               sysname imports-sans-v defaults state-index-map 
3461               const-defs asgn-eq-defs rate-eq-defs 
3462               reaction-eq-defs i-eqs pool-ions perm-ions 
3463               constraints
3464               indent indent+)
3465              (pp indent ,nl)
3466              ))
3467         
3468         
3469          ;; buffer and node initialization
3470          (with-output-to-port cpp-output
3471            (lambda ()
3472              (output-buffers
3473               sysname method state-index-map 
3474               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3475               reaction-eq-defs i-eqs pool-ions perm-ions 
3476               synapse-info defaults
3477               indent indent+)
3478              (pp indent ,nl)
3479              ))
3480         
3481
3482          (with-output-to-port cpp-output
3483            (lambda ()
3484              ((case method
3485                 ((leapfrog) 
3486                  output-nodes-leapfrog)
3487                 ((gsl)     
3488                  output-nodes-gsl)
3489                 ((cvode)
3490                  output-nodes-cvode)
3491                 ((ida)
3492                  output-nodes-ida)
3493                 )
3494               sysname method state-index-map 
3495               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3496               reaction-eq-defs i-eqs pool-ions perm-ions 
3497               synapse-info defaults abstol reltol maxstep
3498               indent indent+)
3499              (pp indent ,nl)
3500              ))
3501         
3502          ;; spike threshold detect (cvode and ida methods only)
3503          (case method
3504            ((cvode)
3505              (with-output-to-port cpp-output
3506                (lambda ()
3507                  (output-cvode-event
3508                   sysname method imports const-defs state-index-map 
3509                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3510                   defaults i-eqs v-eq indent indent+)
3511                  (pp indent ,nl)
3512                  )))
3513            ((ida)
3514              (with-output-to-port cpp-output
3515                (lambda ()
3516                  (output-ida-event
3517                   sysname method imports const-defs state-index-map 
3518                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
3519                   defaults i-eqs v-eq indent indent+)
3520                  (pp indent ,nl)
3521                  )))
3522            )
3523
3524          ;; synaptic event transients
3525
3526          (with-output-to-port cpp-output
3527            (lambda ()
3528              (output-synapse-transients
3529               sysname method state-index-map
3530               const-defs transient-event-defs event-external-eq-defs 
3531               synapse-info imports constraints
3532               indent indent+)
3533              (pp indent ,nl)
3534              ))
3535
3536          ;; state update
3537          (with-output-to-port cpp-output
3538            (lambda ()
3539              (output-update
3540               sysname method state-index-map 
3541               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3542               reaction-eq-defs transient-event-defs i-eqs
3543               pool-ions perm-ions synapse-info 
3544               event-external-eq-defs
3545               defaults imports abstol reltol
3546               indent indent+)
3547              (pp indent ,nl)
3548              ))
3549         
3550          ;; spike handler
3551          (with-output-to-port cpp-output
3552            (lambda ()
3553              (output-event-handle
3554               sysname state-index-map 
3555               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3556               reaction-eq-defs i-eqs pool-ions perm-ions 
3557               synapse-info
3558               indent indent+)
3559              (pp indent ,nl)
3560              ))
3561         
3562          (with-output-to-port cpp-output
3563            (lambda () (pp indent "}" ,nl)))
3564
3565             (with-output-to-port hpp-output
3566               (lambda ()
3567                 (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
3568                 (output-header
3569                  sysname method state-index-map steady-state-index-map defaults
3570                  external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
3571                  reaction-eq-defs i-eqs pool-ions perm-ions 
3572                  synapse-info
3573                  indent indent+)
3574                 (pp indent ,nl)
3575                 ))
3576             
3577             (close-output-port cpp-output)
3578             (close-output-port hpp-output)
3579             
3580             ))
3581      ))
3582  ))
3583
3584))
Note: See TracBrowser for help on using the repository browser.