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

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

nemo: added several variants of membrane potential equation

File size: 66.3 KB
Line 
1;;       
2;; TODO: implement a new component type to be used for thresholds
3;; TODO: finish gsl-fminimizer
4;;
5;; An extension for translating NEMO models to NEST code.
6;;
7;; Copyright 2011-2013 Ivan Raikov and the Okinawa Institute of Science and Technology
8;;
9;; This program is free software: you can redistribute it and/or
10;; modify it under the terms of the GNU General Public License as
11;; published by the Free Software Foundation, either version 3 of the
12;; License, or (at your option) any later version.
13;;
14;; This program is distributed in the hope that it will be useful, but
15;; WITHOUT ANY WARRANTY; without even the implied warranty of
16;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17;; General Public License for more details.
18;;
19;; A full copy of the GPL license can be found at
20;; <http://www.gnu.org/licenses/>.
21;;
22
23(module nemo-nest
24
25        (nemo:nest-translator)
26
27        (import scheme chicken utils data-structures lolevel ports extras srfi-1 srfi-13 srfi-69)
28       
29        (require-extension lolevel matchable strictly-pretty 
30                           varsubst datatype nemo-core nemo-utils
31                           nemo-gate-complex nemo-synapse nemo-defaults)
32
33(define nest-builtin-consts  `(params))
34
35
36(define C++-ops
37  `(+ - * / > < <= >= =))
38
39
40(define builtin-fns
41  `(+ - * / pow neg abs atan asin acos sin cos exp ln
42      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
43      > < <= >= = and or round ceiling floor max min
44      ))
45
46
47(define (nest-name s)
48  (let ((cs (string->list (->string s))))
49    (let loop ((lst (list)) (cs cs))
50      (if (null? cs) (string->symbol (list->string (reverse lst)))
51          (let* ((c (car cs))
52                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
53                           (else #\_))))
54            (loop (cons c1 lst) (cdr cs)))))))
55
56           
57(define (rhsexpr/C++ expr)
58  (match expr 
59         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/C++ x)) es)))
60         (('pow x y)  (if (and (integer? y)  (positive? y))
61                          (if (> y 1)  (let ((tmp (gensym "x")))
62                                         `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
63                              x)
64                            (if (and (number? y) (zero? y)) 1.0 expr)))
65         ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (nest-name s))
66                                             (map (lambda (x) (rhsexpr/C++ x)) es)) expr))
67         (id          (if (symbol? id) (nest-name id) id))))
68
69
70(define (nest-state-name n s)
71  (nest-name (s+ n s)))
72
73
74(define-syntax pp
75  (syntax-rules ()
76    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
77
78
79(define group/C++   (doc:block 2 (doc:text "(") (doc:text ")")))
80(define block/C++   (doc:block 2 (doc:text "{") (doc:text "}")))
81(define (stmt/C++ x) 
82  (match x
83         (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
84         (else  (doc:cons x (doc:text ";")))))
85
86
87(define (ifthen/C++ c e1 e2)
88  (doc:nest 2
89    (doc:connect (doc:group (doc:connect (doc:text "if") c))
90                 (doc:connect (doc:nest 2 e1)
91                              (doc:nest 2 (doc:connect 
92                                           (doc:text "else") 
93                                           e2))))
94    ))
95
96
97(define (letblk/C++ e1 e2)
98  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
99        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
100        (else (doc:connect (doc:group (doc:nest 2 (stmt/C++ e1)))
101                           (doc:group (doc:nest 2 e2))))))
102       
103
104(define (format-op/C++ indent op args)
105  (let ((op1 (doc:text (->string op))))
106    (if (null? args) op1
107        (match args
108               ((x)      (doc:concat (list op1 x)))
109               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
110               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
111               (lst      (let* ((n   (length lst))
112                                (n/2 (inexact->exact (round (/ n 2)))))
113                           (doc:concat 
114                            (intersperse 
115                             (list (format-op/C++ indent op (take lst n/2 )) op1 
116                                   (format-op/C++ indent op (drop lst n/2 )))
117                             (doc:space)))))))))
118
119(define (format-fncall/C++ indent op args)
120  (let ((op1 (doc:text (->string op))))
121    (doc:cons op1 (group/C++ ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
122
123(define (name-normalize expr)
124  (match expr 
125         (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
126         (('let bs e)
127          `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
128         ((f . es) 
129          (cons (if (member f builtin-fns) f (nest-name f)) (map name-normalize es)))
130         ((? symbol? ) (nest-name expr))
131         ((? atom? ) expr)))
132
133
134(define (canonicalize-expr/C++ expr)
135  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term)))
136    (let* ((expr1 (if-convert expr))
137           (expr2 (subst-convert expr1 subst-empty))
138           (expr3 (let-lift expr2))
139           (expr4 (name-normalize expr3)))
140      expr4)))
141
142
143(define (format-expr/C++ indent expr . rest) 
144  (let-optionals rest ((rv #f))
145   (let ((indent+ (+ 2 indent)))
146    (match expr
147       (('let bindings body)
148        (letblk/C++
149         (fold-right 
150           (lambda (x ax)
151             (letblk/C++
152              (match (second x)
153                     (('if c t e)
154                      (ifthen/C++
155                       (group/C++ (format-expr/C++ indent c))
156                       (block/C++ (format-expr/C++ indent t (first x)))
157                       (block/C++ (format-expr/C++ indent e (first x)))))
158                     (else
159                      (stmt/C++
160                       (format-op/C++ indent+ " = "
161                                         (list (format-expr/C++ indent (first x) )
162                                               (format-expr/C++ indent (second x)))))))
163              ax))
164           (doc:empty) bindings)
165         (match body
166                (('let _ _) (format-expr/C++ indent body rv))
167                (else
168                 (let ((body1 (doc:nest indent (format-expr/C++ indent body))))
169                   (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) body1)))
170                       body1))))))
171       
172       (('if . rest) (error 'format-expr/C++ "invalid if statement " expr))
173
174       ((op . rest) 
175       (let ((op (case op ((ln) 'log) ((abs) 'fabs) (else op))))
176         (let ((fe
177                (if (member op C++-ops)
178                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
179                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
180                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
181                      (case op
182                        ((/) 
183                         (format-op/C++ indent op 
184                                          (map (lambda (x) 
185                                                 (let ((fx (format-expr/C++ indent+ x)))
186                                                   (if (or (symbol? x) (number? x)) fx
187                                                       (if (or mul? plmin?) (group/C++ fx) fx)))) rest)))
188                        ((*) 
189                         (format-op/C++ indent op 
190                                          (map (lambda (x) 
191                                                 (let ((fx (format-expr/C++ indent+ x)))
192                                                   (if (or (symbol? x) (number? x)) fx
193                                                       (if plmin? (group/C++ fx) fx)))) rest)))
194                       
195                        (else
196                         (format-op/C++ indent op 
197                                          (map (lambda (x) 
198                                                 (let ((fx (format-expr/C++ indent+ x))) fx)) rest)))))
199                   
200                    (let ((op (case op ((neg) '-) (else op))))
201                      (format-fncall/C++ indent op (map (lambda (x) (format-expr/C++ indent+ x)) rest))))))
202           (if rv 
203               (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
204               fe))))
205     
206      (else  (let ((fe (doc:text (->string expr))))
207               (if rv 
208                   (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
209                   fe)))))))
210               
211         
212(define (expr->string/C++ x . rest)
213  (let-optionals rest ((rv #f) (width 72))
214    (sdoc->string (doc:format width (format-expr/C++ 2 x rv)))))
215
216
217
218(define (add-params-to-fncall expr #!key (name 'params))
219  (match expr 
220         (('if c t e)  `(if ,(add-params-to-fncall c) ,(add-params-to-fncall t) ,(add-params-to-fncall e)))
221         (('let bs e)
222          `(let ,(map (lambda (b) `(,(car b) ,(add-params-to-fncall (cadr b)))) bs) ,(add-params-to-fncall e)))
223         ((f . es) 
224          (if (member f builtin-fns)
225              (cons  f  (map add-params-to-fncall es))
226              (let ((es1 (map add-params-to-fncall es)))
227                (cons f (if (equal? (last es1) name) es1 (append es1 `(,name))) ))
228              ))
229         ((or (? symbol? ) (? atom? )) expr)))
230
231
232
233(define (state-init n init)
234  (let* ((init  (rhsexpr/C++ init))
235         (init1 (canonicalize-expr/C++ 
236                 (add-params-to-fncall init))))
237    (list (nest-name n) init1)))
238
239
240(define (asgn-eq n rhs)
241  (let* ((fbody   (rhsexpr/C++ rhs))
242         (fbody1  (canonicalize-expr/C++ 
243                   (add-params-to-fncall fbody))))
244    (list (nest-name n) fbody1)))
245
246
247(define (reaction-eq n open transitions conserve)
248  (if (symbol? open)
249      (list (nest-name n) (nest-state-name n open))
250      (list (nest-name n) (sum (map (lambda (x) (nest-state-name n x)) open)))
251      ))
252
253
254(define (external-eq-defs sys)
255  (fold  (lambda (n ax) 
256           (let ((en (hash-table-ref sys n)))
257             (if (nemo:quantity? en)
258                 (cases nemo:quantity en
259                        (EXTERNAL (local-name name namespace) 
260                                  (if (not (equal? local-name name))
261                                      (cons (asgn-eq local-name name) ax) ax))
262                        (else  ax))
263                 ax)))
264         (list) (hash-table-keys sys)))
265
266(define (reaction-transition-eqs n initial open transitions conserve power)
267  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nest-state-name)))
268     (let* ((out-edges  (g 'out-edges))
269            (in-edges   (g 'in-edges))
270            (nodes      ((g 'nodes))))
271       ;; generate differential equations for each state in the transitions system
272
273       (let ((eqs    (fold (lambda (s ax) 
274                            (if (and cnode (= (first cnode) (first s) )) ax
275                                (let* ((out   (out-edges (first s)))
276                                       (in    (in-edges (first s)))
277                                       (open? (eq? (second s) open))
278                                       (name  (nest-name (lookup-def (second s) node-subs))))
279                                  (let* ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
280                                                       `(+ (neg ,(sum (map third out)))
281                                                           ,(sum (map third in))))
282                                                      ((and (not (null? out)) (null? in))
283                                                       `(neg ,(sum (map third out))))
284                                                      ((and (null? out) (not (null? in)))
285                                                       (sum (map third in)))))
286                                         (fbody0 (rhsexpr/C++ rhs1)))
287                       
288                                    (cons (list name (canonicalize-expr/C++ 
289                                                      (add-params-to-fncall fbody0))) ax)
290                                    ))))
291                          (list) nodes)))
292        eqs))))
293           
294       
295
296(define (poset->asgn-eq-defs poset sys)
297  (fold-right
298   (lambda (lst ax)
299     (fold  (lambda (x ax) 
300              (match-let (((i . n)  x))
301                         (let ((en (hash-table-ref sys n)))
302                           (if (nemo:quantity? en)
303                               (cases nemo:quantity en
304                                      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
305                                      (else  ax))
306                               ax))))
307            ax lst))
308   (list) poset))
309
310
311(define (poset->rate-eq-defs poset sys)
312  (fold-right
313   (lambda (lst ax)
314     (fold  (lambda (x ax) 
315              (match-let (((i . n)  x))
316                         (let ((en (hash-table-ref sys n)))
317                           (if (nemo:quantity? en)
318                               (cases nemo:quantity en
319
320                                      (REACTION  (name initial open transitions conserve power) 
321                                                 (append (reaction-transition-eqs name initial open transitions 
322                                                                                  conserve power) ax))
323                                     
324                                      (RATE (name initial rhs power)
325                                            (let ((fbody0  (rhsexpr/C++ rhs))
326                                                  (dy      (nest-name name) ))
327                                              (cons (list dy (canonicalize-expr/C++ 
328                                                              (add-params-to-fncall fbody0))) ax)))
329
330                                      (else  ax))
331                               ax))))
332            ax lst))
333   (list) poset))
334
335
336(define (poset->reaction-eq-defs poset sys)
337  (fold-right
338   (lambda (lst ax)
339     (fold  (lambda (x ax) 
340              (match-let (((i . n)  x))
341                         (let ((en (hash-table-ref sys n)))
342                           (if (nemo:quantity? en)
343                               (cases nemo:quantity en
344                                      (REACTION  (name initial open transitions conserve power) 
345                                                 (cons (reaction-eq name open transitions conserve) ax))
346                                      (else  ax))
347                               ax))))
348            ax lst))
349   (list) poset))
350
351
352(define (poset->init-defs poset sys)
353  (fold-right
354   (lambda (lst ax)
355     (fold-right
356      (lambda (x ax) 
357              (match-let (((i . n)  x))
358                         (let ((en (hash-table-ref sys n)))
359                           (if (nemo:quantity? en)
360                               (cases nemo:quantity en
361                                      (REACTION  (name initial open transitions conserve power)
362                                                 (if (nemo:rhs? initial)
363                                                     (cons* (state-init name initial) 
364                                                            (state-init (nest-state-name name open) name) ax) 
365                                                     ax))
366
367                                      (RATE  (name initial rhs power)
368                                             (if (and initial (nemo:rhs? initial))
369                                                 (cons (state-init name initial) ax) 
370                                                 ax))
371
372                                      (else  ax))
373                               ax))))
374            ax lst))
375   (list) poset))
376
377
378(define (poset->state-conserve-eq-defs poset sys)
379  (fold-right
380   (lambda (lst ax)
381     (fold  (lambda (x ax) 
382              (match-let (((i . n)  x))
383                         (let ((en (hash-table-ref sys n)))
384                           (if (nemo:quantity? en)
385                               (cases nemo:quantity en
386                                      (REACTION (name initial open transitions conserve power)
387                                                (if (and (list? conserve) (every nemo:conseq? conserve))
388                                                    (cons (state-conseqs (nest-name name) transitions conserve
389                                                                        nest-state-name) ax) 
390                                                    ax))
391                                      (else  ax))
392                               ax))))
393            ax lst))
394   (list) poset))
395
396
397(define (find-locals defs)
398  (concatenate
399   (map (lambda (def) 
400          (match def 
401                 (('let bnds body) 
402                  (let ((bexprs (map second bnds)))
403                    (concatenate (list (map first bnds) 
404                                       (find-locals bexprs )
405                                       (find-locals (list body))))))
406                 (('if c t e)      (append (find-locals (list t)) (find-locals (list e))))
407                 ((s . rest)       (find-locals rest))
408                 (else (list)))) 
409        defs)))
410
411
412(define (rate/reaction-power sys n)
413  (let ((en (hash-table-ref sys n)))
414    (if (nemo:quantity? en)
415        (cases nemo:quantity en
416               (REACTION  (name initial open transitions conserve power) 
417                          power)
418               (RATE    (name initial rhs power)
419                        power)
420               (else  #f)) 
421        #f)))
422
423
424(define (bucket-partition p lst)
425  (let loop ((lst lst) (ax (list)))
426    (if (null? lst) ax
427        (let ((x (car lst)))
428          (let bkt-loop ((old-bkts ax) (new-bkts (list)))
429            (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
430                (if (p x (caar old-bkts ))
431                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
432                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
433
434
435(define (make-define-fn sysname )
436  (lambda (indent n proc)
437
438    (let ((lst (procedure-data proc))
439          (indent+ (+ 2 indent)))
440
441      (let ((rt       (or (lookup-def 'rt lst) 'double))
442            (formals  (lookup-def 'formals lst))
443            (vars     (lookup-def 'vars lst))
444            (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
445            (body     (lookup-def 'body lst))
446            (rv       (gensym 'rv)))
447
448        (let ((argument-list 
449               (append
450                (if (null? vars) '() (map (lambda (x) (s+ "double " (nest-name x))) vars))
451                '("const void* params"))))
452          (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) "{" ))
453          (let* ((body0 (rhsexpr/C++ body))
454                 (body1 (canonicalize-expr/C++ (add-params-to-fncall body0)))
455                 (lbs   (enum-bnds body1 (list))))
456            (pp indent+ (double ,rv ";"))
457            (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs) ";")))
458            (if (not (null? consts)) 
459                (begin (pp indent+ (double ,(slp ", " (map (compose nest-name car) consts)) ";")
460                           ("const " ,(s+ sysname "::Parameters_") "& p =  *(reinterpret_cast< const " ,(s+ sysname "::Parameters_") "*>(params));"))
461                       (for-each (lambda (x) (let ((n (car x))) (pp indent+ (,(nest-name n) = ,(s+ "p." (nest-name n) ";"))))) consts)
462                       ))
463            (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
464            (pp indent+ ,(s+ "return " rv ";"))
465            ))
466        (pp indent "}"))
467    )))
468
469
470(define (make-define-fn-header sysname )
471  (lambda (indent n proc)
472
473    (let ((lst (procedure-data proc))
474          (indent+ (+ 2 indent)))
475
476      (let ((rt       (or (lookup-def 'rt lst) 'double))
477            (formals  (lookup-def 'formals lst))
478            (vars     (lookup-def 'vars lst))
479            (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
480            (body     (lookup-def 'body lst))
481            (rv       (gensym 'rv)))
482
483        (let ((argument-list 
484               (append
485                (if (null? vars) '() (map (lambda (x) (s+ "double " (nest-name x))) vars))
486                '("const void* params"))))
487          (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) ";" ))
488          ))
489      )))
490
491
492
493(define (ith v i) (sprintf "Ith(~A,~A)" v i))
494
495
496(define (output-dy sysname method imports const-defs state-index-map 
497                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
498                   pool-ions mcap i-eqs v-eq indent indent+)
499
500
501  (case method
502    ((leapfrog)
503     (pp indent ,nl (
504                     ,(s+ "extern \"C\" int " sysname "_dynamics")
505                     (,(slp ", " `("double t"
506                                   "const double y[]"
507                                   "double f[]"
508                                   "void* pnode"
509                                   )))
510                     #\{
511                     )))
512    ((cvode)
513     (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
514                     (,(slp ", " `("double t"
515                                   "N_Vector y"
516                                   "N_Vector f"
517                                   "void* pnode"
518                                   )))
519                     #\{
520                     )))
521    ((gsl)
522     (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
523                     (,(slp ", " `("double t"
524                                   "const double y[]"
525                                   "double f[]"
526                                   "void* pnode"
527                                   )))
528                     #\{
529                     )))
530     )
531
532  (let* (
533         (asgn-eq-defs 
534          (map
535           (lambda (def) (list (first def) 
536                               (add-params-to-fncall 
537                                (canonicalize-expr/C++ (second def)))))
538           asgn-eq-defs))
539
540         (rate-eq-defs 
541          (map
542           (lambda (def) (list (first def) 
543                               (add-params-to-fncall 
544                                (canonicalize-expr/C++ (second def)))))
545           rate-eq-defs))
546
547         (reaction-eq-defs 
548          (map
549           (lambda (def) (list (first def) 
550                               (add-params-to-fncall 
551                                (canonicalize-expr/C++ (second def)))))
552           reaction-eq-defs))
553
554         (i-eqs 
555          (map
556           (lambda (def) (list (first def) 
557                               (add-params-to-fncall (canonicalize-expr/C++ (second def)))))
558           i-eqs))
559
560         (v-eq
561          (and v-eq
562               (list (first v-eq) 
563                     (add-params-to-fncall (canonicalize-expr/C++ (second v-eq))))))
564
565         (eqs 
566          (append
567           
568           asgn-eq-defs
569           reaction-eq-defs
570                 
571           (map (lambda (pool-ion)
572                  (let ((n (pool-ion-in pool-ion))
573                        (b (pool-ion-inq pool-ion)))
574                    (list n b)))
575                pool-ions)
576           (map (lambda (pool-ion)
577                  (let ((n (pool-ion-out pool-ion))
578                        (b (pool-ion-outq pool-ion)))
579                    (list n b)))
580                pool-ions)
581          ))
582         
583         (eq-dag 
584          (map (lambda (def)
585                 (cons (first def) (enum-freevars (second def) '() '())))
586               eqs))
587
588         (eq-order
589          (reverse
590           (topological-sort eq-dag 
591                             (lambda (x y) (string=? (->string x) (->string y))))))
592
593         (eq-locals  (find-locals 
594                      (map second
595                           (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
596                               (append i-eqs rate-eq-defs eqs)))))
597         )
598
599
600    (pp indent+ (double ,(slp ", " (delete-duplicates 
601                                    (map (compose ->string nest-name)
602                                         (filter (lambda (x) 
603                                                   (not (member x nest-builtin-consts)))
604                                                 (append
605                                                  eq-locals
606                                                  eq-order
607                                                  (map first i-eqs)
608                                                  (map first external-eq-defs)
609                                                  (map first state-index-map)
610                                                  (map first const-defs)
611                                                  )))
612                                    string=?))
613                        ";"))
614
615    (pp indent+ ,nl
616        ("// S is shorthand for the type that describes the model state ")
617        (typedef ,(s+ sysname "::State_") "S;")
618       
619        ,nl
620       
621        ("// cast the node ptr to an object of the proper type")
622        ("assert(pnode);")
623        ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
624        (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
625       
626        ,nl
627
628        ("// params is a reference to the model parameters ")
629        (const struct ,(s+ sysname "::Parameters_") "*params;")
630        ("params = &(node.P_);")
631       
632        ,nl
633
634        ("// state is a reference to the model state ")
635        (struct ,(s+ sysname "::State_") "*state;")
636        ("state = &(vnode.S_);")
637       
638        ,nl
639
640        ("// y[] must be the state vector supplied by the integrator, ")
641        ("// not the state vector in the node, node.S_.y[]. ")
642        ,nl
643        )
644
645   
646    (for-each (lambda (n) (pp indent+ ,(expr->string/C++ (s+ "state->" n) n)))
647              (append
648               (map pool-ion-in pool-ions)
649               (map pool-ion-out pool-ions)
650              ))
651
652    (for-each (lambda (n) (pp indent+ ,(expr->string/C++ (s+ "state->" n) n)))
653              (map first i-eqs))
654
655    (for-each (lambda (def)
656                (let ((n (first def))
657                      (b (second def)))
658                  (pp indent+ ,(expr->string/C++ b (nest-name n)))))
659              external-eq-defs)
660   
661    (for-each (lambda (def)
662                (let ((n (first def)) )
663                  (pp indent+
664                      ,(expr->string/C++ 
665                        (s+ "params->" n) n))))
666              const-defs)
667
668    (let ((vi (lookup-def 'v state-index-map)))
669      (if vi (pp indent+ ,(expr->string/C++ (ith 'y vi) 'v))))
670
671    (for-each (lambda (def)
672                (let* ((n (first def))
673                       (ni (lookup-def n state-index-map)))
674                  (pp indent+ ,(expr->string/C++ (ith 'y ni) (nest-name n)))))
675              rate-eq-defs)
676
677    (for-each (lambda (n)
678                (let ((b  (lookup-def n eqs)))
679                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
680              eq-order)
681
682    (for-each (lambda (def)
683                (let* ((n (first def))
684                       (b (second def))
685                       (fv (ith 'f (lookup-def n state-index-map)))
686                       )
687                  (pp indent+ ,(s+ "// rate equation " n)
688                      ,(expr->string/C++ b fv))))
689              rate-eq-defs)
690   
691    (for-each (lambda (def) 
692                (pp indent+ ,(expr->string/C++ (second def) (first def)))) 
693            i-eqs)
694   
695    (if v-eq
696        (let ((vi (lookup-def 'v state-index-map)))
697          (if vi
698              (pp indent+ ,(expr->string/C++ (second v-eq) (ith 'f vi)))
699              )))
700
701   
702    (for-each (lambda (n) (pp indent+ ,(expr->string/C++ n (s+ "state->" n))))
703              (append (map pool-ion-in pool-ions)
704                      (map pool-ion-out pool-ions)))
705
706   
707    (for-each (lambda (n) (pp indent+ ,(expr->string/C++ n (s+ "state->" n))))
708              (map first i-eqs))
709
710    (for-each (lambda (n) (pp indent+ ,(expr->string/C++ (s+ "state->" n) n)))
711              (map first i-eqs))
712
713    (case method
714      ((leapfrog)
715        (pp indent+ ,nl ("return 0;")))
716      ((cvode)
717        (pp indent+ ,nl ("return 0;")))
718      ((gsl)
719        (pp indent+ ,nl ("return GSL_SUCCESS;")))
720      )
721    (pp indent  #\})
722
723    ))
724
725
726
727(define (output-parameters sysname imports const-defs external-eq-defs defaults indent indent+)
728
729  (let* ((parameter-defs
730          (map
731           (lambda (def) (list (first def) 
732                               (add-params-to-fncall (canonicalize-expr/C++ (second def)))))
733           const-defs))
734         (parameter-locals  (find-locals (map second parameter-defs))))
735
736    (pp indent ,nl (,(s+ sysname "::Parameters_::Parameters_") () ":"))
737
738    (let ((const-exprs
739            (map (lambda (def)
740                   (let* ((n  (first def)) (b (second def)))
741                     (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
742                 const-defs) )
743             
744          (default-exprs 
745            (map (lambda (def)
746                   (let ((n (first def))
747                         (b (second def)))
748                     (expr->string/C++ (nest-name b) (nest-name n))))
749                 defaults))
750           )
751     
752      (if (not (null? parameter-locals)) 
753          (pp indent+ (double ,(slp ", " parameter-locals) ";")))
754
755      (pp indent+ ,(slp ",\n" const-exprs))
756     
757      (pp indent  #\{ ,(slp "\n" default-exprs) #\})
758     
759      ))
760  )
761
762
763(define (output-init sysname state-index-map steady-state-index-map 
764                     imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
765                     reaction-eq-defs i-eqs pool-ions perm-ions defaults indent indent+)
766
767 
768  (let* ((N (length state-index-map))
769
770         (asgn-eq-defs
771          (map
772           (lambda (def) (list (first def) 
773                               (add-params-to-fncall 
774                                (canonicalize-expr/C++ (second def)))))
775           asgn-eq-defs))
776         
777         (init-eq-defs
778          (map
779           (lambda (def) (list (first def) 
780                               (add-params-to-fncall 
781                                (canonicalize-expr/C++ (second def)))))
782           init-eq-defs))
783
784         (reaction-eq-defs 
785          (map
786           (lambda (def) (list (first def) 
787                               (add-params-to-fncall 
788                                (canonicalize-expr/C++ (second def)))))
789           reaction-eq-defs))
790
791         (i-eqs 
792          (map
793           (lambda (def) (list (first def) 
794                               (add-params-to-fncall 
795                                (canonicalize-expr/C++ (second def)))))
796           i-eqs))
797
798         (init-eqs 
799          (append
800             
801           (map (lambda (def)
802                  (let ((n (first def))
803                        (b (second def)))
804                    (list (nest-name n) (nest-name b))))
805                external-eq-defs)
806           
807           asgn-eq-defs
808           init-eq-defs
809
810           (map (lambda (pool-ion)
811                  (let ((n (pool-ion-in pool-ion))
812                        (b (pool-ion-inq pool-ion)))
813                    (list n b)))
814                pool-ions)
815
816           (map (lambda (pool-ion)
817                  (let ((n (pool-ion-out pool-ion))
818                        (b (pool-ion-outq pool-ion)))
819                    (list n b)))
820                pool-ions)
821           ))
822
823         (init-dag 
824          (map (lambda (def)
825                 (cons (first def) (enum-freevars (second def) '() '())))
826               init-eqs))
827
828         (init-order
829          (reverse
830           (topological-sort init-dag 
831                             (lambda (x y) (string=? (->string x) (->string y))))))
832
833         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
834
835         )
836   
837    (pp indent ,nl (,(s+ sysname "::State_::State_") (const Parameters_& p) ": r_(0)"))
838   
839    (pp indent  #\{)
840
841    (pp indent+
842        (double ,(slp ", " (delete-duplicates
843                            (map ->string 
844                                 (filter (lambda (x) (and (not (member x nest-builtin-consts))
845                                                          (not (member x (map second imports)))))
846                                         (append
847                                          init-locals
848                                          init-order
849                                          (map first external-eq-defs)
850                                          (map first i-eqs)
851                                          (map first state-index-map) 
852                                          (map first const-defs)
853                                          (map first reaction-eq-defs)
854                                          )))
855                                 string=?))
856                      ";")
857        ("const Parameters_ *params = &p;")
858        )
859
860    (pp indent+ ,(sprintf "memset(y_,0,~A*sizeof(double));" N))
861
862    (for-each
863     (lambda (n) (pp indent+ ,(expr->string/C++ (s+ "p." n) n)))
864     (map (compose nest-name first) const-defs) )
865
866    (let ((vi (lookup-def 'v state-index-map))
867          (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
868      (if (and vi vrest) 
869          (pp indent+ ,(expr->string/C++  vrest 'v ))))
870
871    (for-each (lambda (n)
872                (let ((b  (lookup-def n init-eqs)))
873                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
874              init-order)
875   
876    (for-each (lambda (def)
877                (let* ((n  (first def)) 
878                       (ni (lookup-def n state-index-map)))
879                  (if ni (pp indent+ ,(expr->string/C++ n  (s+ "y_[" ni "]"))))))
880              init-eq-defs)
881
882   
883#|
884    (if (not (null? steady-state-index-map))
885        (begin
886          (gsl-fminimizer sysname N indent+)
887          (for-each
888           (lambda (def)
889             (let* ((n (first def))
890                    (si (lookup-def n steady-state-index-map))
891                    (ni (lookup-def n state-index-map)))
892               (if si (begin
893                        (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n))
894                        (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")")
895                                                       (s+ "y_[" ni "]")))))))
896           rate-eq-defs)
897      (pp indent+ "gsl_vector_free (ys);")
898      ))
899|#
900
901    (for-each
902     (lambda (def)
903       (let ((n (first def)) (b (second def)))
904         (if (not (lookup-def n init-eq-defs))
905             (pp indent+ ,(expr->string/C++ b n)))))
906     reaction-eq-defs)
907   
908    (for-each
909     (lambda (def) 
910       (pp indent+ ,(expr->string/C++ (second def) (first def))))
911     i-eqs)
912
913    (let ((vi (lookup-def 'v state-index-map)))
914      (if vi (pp indent+ ,(expr->string/C++ 'v (s+ "y_[" vi "]")))))
915     
916    (pp indent  #\})
917
918    (pp indent ,nl
919        (,(s+ sysname "::State_::State_") (const State_& s) ": " )
920        (,(slp ",\n" (map (lambda (n) (sprintf "~A(s.~A)" n n))
921                          (cons "r_"
922                                (delete-duplicates
923                                 (append
924                                  (map pool-ion-in pool-ions)
925                                  (map pool-ion-out pool-ions)
926                                  )
927                                 )
928                                ))
929               ))
930        )
931
932    (pp indent  #\{)
933    (pp indent+ (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N)))
934
935    (pp indent  #\})
936
937    (pp indent ,nl (,(s+ sysname "::State_& " sysname "::State_::operator=(const State_& s)")))
938
939    (pp indent  #\{)
940
941
942    (pp indent+ (,(sprintf #<<EOF
943   assert(this != &s)
944     for ( size_t i = 0 ; i < ~A ; ++i )
945       y_[i] = s.y_[i];
946EOF
947N)))
948
949    (pp indent+
950        (,(slp "\n" (map (lambda (n) (expr->string/C++ (sprintf "s.~A" n) n))
951                          (cons "r_"
952                                (delete-duplicates
953                                 (append
954                                  (map pool-ion-in pool-ions)
955                                  (map pool-ion-out pool-ions)
956                                  ))
957                                ))
958               ))
959        )
960
961     (pp indent+ "return *this;")
962
963     (pp indent  #\})
964))
965
966
967
968(define (gsl-fminimizer sysname N indent)
969
970  (pp indent (,(sprintf #<<EOF
971
972  int minimizer_status, minimizer_iterc;
973
974  const gsl_multimin_fdfminimizer_type *T;
975  gsl_multimin_fdfminimizer *S;
976
977  gsl_vector *ys;
978  gsl_multimin_function_fdf F;
979
980  F.f   = &~A_dynamics;
981  F.df  = &~A_dynamics;
982  F.fdf = &~A_dynamics;
983  F.n   = ~A;
984  F.params = (void *)&p;
985
986  ys = gsl_vector_alloc (~A);
987
988  for (int i = 0; i < ~A; i++)
989  {
990     gsl_vector_set (ys, i, 0.0);
991  }
992
993  T = gsl_multimin_fdfminimizer_conjugate_fr;
994  S = gsl_multimin_fdfminimizer_alloc (T, ~A);
995
996  gsl_multimin_fdfminimizer_set (S, &F, ys, 0.01, 1e-4);
997
998  do
999    {
1000      minimizer_iterc++;
1001      minimizer_status = gsl_multimin_fdfminimizer_iterate (S);
1002
1003      if (minimizer_status)
1004        break;
1005
1006      minimizer_status = gsl_multimin_test_gradient (S->gradient, 1e-3);
1007
1008    }
1009  while (minimizer_status == GSL_CONTINUE && minimizer_iterc < 100);
1010
1011  gsl_multimin_fdfminimizer_free (S);
1012
1013EOF
1014sysname sysname sysname N N N N)))
1015)
1016
1017
1018
1019(define (output-accessors+modifiers
1020         sysname imports state-index-map const-defs asgn-eq-defs rate-eq-defs 
1021         reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
1022
1023  (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
1024  (pp indent  #\{)
1025
1026  (for-each
1027   (lambda (def)
1028     (let ((n (first def)))
1029       (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
1030   const-defs)
1031
1032  (pp indent  #\})
1033
1034  (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
1035  (pp indent  #\{)
1036
1037  (for-each
1038   (lambda (def)
1039     (let ((n (first def)))
1040       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))))
1041   const-defs)
1042
1043  (pp indent  #\})
1044
1045  (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
1046  (pp indent  #\{)
1047
1048  (for-each
1049   (lambda (n)
1050     (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n))))
1051   (delete-duplicates
1052    (append (map first i-eqs)
1053            (map pool-ion-in pool-ions)
1054            (map pool-ion-out pool-ions))
1055    ))
1056
1057  (for-each
1058   (lambda (def)
1059     (let ((n (first def)) (i (second def)))
1060       (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
1061   state-index-map)
1062
1063  (pp indent  #\})
1064
1065  (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
1066  (pp indent  #\{)
1067
1068  (for-each
1069   (lambda (n)
1070     (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n))))
1071   (delete-duplicates
1072    (append (map first i-eqs)
1073            (map pool-ion-in pool-ions)
1074            (map pool-ion-out pool-ions))))
1075
1076  (for-each
1077   (lambda (def)
1078     (let ((n (first def)) (i (second def)))
1079       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
1080   state-index-map)
1081
1082
1083  (pp indent  #\})
1084
1085)
1086
1087
1088
1089(define (output-handle
1090         sysname state-index-map steady-state-index-map 
1091         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1092         reaction-eq-defs i-eqs pool-ions perm-ions 
1093         synapse-info
1094         indent indent+)
1095 
1096  (pp indent  ,nl (,(sprintf #<<EOF
1097void ~A::handle(SpikeEvent & e)
1098  {
1099    int flag;
1100    assert(e.get_delay() > 0);
1101    flag = 0;
1102EOF
1103   sysname)))
1104
1105   (let ((isyns (lookup-def 'i-synapses synapse-info))
1106          (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
1107      (for-each
1108       (lambda (isyn psc)
1109         (pp indent+ 
1110             (,(sprintf #<<EOF
1111    if ((flag==0) && (e.get_weight() > P_.~A))
1112      {
1113        B_.spike_~A.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
1114                              abs(e.get_weight()) * e.get_multiplicity());
1115        flag = 1;
1116      }
1117EOF
1118             (nest-name (second isyn))
1119             (second psc)))))
1120       isyns pscs))
1121
1122  (pp indent ,nl #\} ,nl)
1123
1124  (pp indent  ,nl (,(sprintf #<<EOF
1125void ~A::handle(CurrentEvent& e)
1126  {
1127    assert(e.get_delay() > 0);
1128
1129    const double_t c=e.get_current();
1130    const double_t w=e.get_weight();
1131
1132    B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), 
1133                        w *c);
1134  }
1135
1136void ~A::handle(DataLoggingRequest& e)
1137  {
1138    B_.logger_.handle(e);
1139  }
1140EOF
1141  sysname sysname)))
1142
1143)
1144
1145
1146(define (output-recordables
1147         sysname state-index-map steady-state-index-map 
1148         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1149         reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
1150 
1151  (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
1152  (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
1153  (pp indent #\{)
1154  (for-each
1155   (lambda (def)
1156     (let ((n (first def)) (i (second def)))
1157       (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);" 
1158                              (->string n) sysname sysname (string-upcase (->string n)))))
1159       ))
1160   state-index-map)
1161
1162  (if (lookup-def 'v state-index-map)
1163      (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);" 
1164                             sysname sysname (string-upcase (->string 'v))))))
1165
1166  (pp indent #\})
1167  )
1168
1169
1170(define (output-prelude sysname steady-state-index-map indent)
1171  (pp indent (,(sprintf "#include \"~A.h\"
1172#include \"exceptions.h\"
1173#include \"network.h\"
1174#include \"dict.h\"
1175#include \"integerdatum.h\"
1176#include \"doubledatum.h\"
1177#include \"dictutils.h\"
1178#include \"numerics.h\"
1179#include <limits>
1180
1181#include \"universal_data_logger_impl.h\"
1182
1183#include <iomanip>
1184#include <iostream>
1185#include <cstdio>
1186#include <cstring>
1187#include <cmath>
1188"  sysname)))
1189  )
1190
1191(define (output-header 
1192         sysname method state-index-map steady-state-index-map 
1193         defaults external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1194         reaction-eq-defs i-eqs pool-ions perm-ions 
1195         synapse-info
1196         indent indent+)
1197
1198  (pp indent ("#include \"nest.h\"
1199#include \"event.h\"
1200#include \"archiving_node.h\"
1201#include \"ring_buffer.h\"
1202#include \"connection.h\"
1203#include \"universal_data_logger.h\"
1204#include \"recordables_map.h\"
1205"))
1206
1207  (case method
1208    ((leapfrog) 
1209     (pp indent
1210         ("#define Ith(v,i)    (v[i])   /* Ith component in a vector */")
1211         ))
1212    ((cvode) 
1213     (pp indent
1214         ("#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */")
1215         ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
1216         ("#include <cvode/cvode_dense.h>       /* prototype for CVDense */")
1217         ("#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */")
1218         ("#include <sundials/sundials_types.h> /* definition of type realtype */")
1219         ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
1220         ))
1221    ((gsl) 
1222     (pp indent
1223         ("#include <gsl/gsl_errno.h>")
1224         ("#include <gsl/gsl_matrix.h>")
1225         ("#include <gsl/gsl_sf_exp.h>")
1226         ("#include <gsl/gsl_odeiv.h>")
1227         ("#define Ith(v,i)    (v[i])")
1228         ))
1229    )
1230
1231  (pp indent ,nl)
1232  (pp indent "namespace nest {" ,nl)
1233
1234 
1235  (case method
1236    ((cvode) 
1237     (pp indent #<<EOF
1238
1239/* 
1240 * Check function return value.
1241 *    opt == 0 means SUNDIALS function allocates memory so check if
1242 *             returned NULL pointer
1243 *    opt == 1 means SUNDIALS function returns a flag so check if
1244 *             flag >= 0
1245 *    opt == 2 means function allocates memory so check if returned
1246 *             NULL pointer 
1247 */
1248
1249static int check_flag(void *flagvalue, const char *funcname, int opt)
1250{
1251  int *errflag;
1252
1253  /* Check if CVode function returned NULL pointer - no memory allocated */
1254  if (opt == 0 && flagvalue == NULL) {
1255    fprintf(stderr, 
1256            "\nCVode error: %s() failed - returned NULL pointer\n\n",
1257            funcname);
1258    return(1); }
1259
1260  /* Check if flag < 0 */
1261  else if (opt == 1) {
1262    errflag = (int *) flagvalue;
1263    if (*errflag < 0) {
1264      fprintf(stderr, 
1265              "\nCVode error: %s() failed with flag = %d\n\n",
1266              funcname, *errflag);
1267      return(1); }}
1268
1269  /* Check if function returned NULL pointer - no memory allocated */
1270  else if (opt == 2 && flagvalue == NULL) {
1271    fprintf(stderr, 
1272            "\nMemory error: %s() failed - returned NULL pointer\n\n",
1273            funcname);
1274    return(1); }
1275
1276  return(0);
1277}
1278EOF
1279)))
1280
1281 
1282  (case method
1283    ((leapfrog) 
1284     (pp indent 
1285         ("typedef int (*rhsfn_t)(double t, const double *y, double *ydot, void *user_data);") ,nl
1286         (extern "\"C\""
1287                     int ,(s+ sysname "_dynamics")\
1288                     "(double, const double *, double *, void*)" #\;) ,nl))
1289    ((cvode) 
1290     (pp indent (extern "\"C\""
1291                     int ,(s+ sysname "_dynamics")
1292                     "(double, const N_Vector, N_Vector, void*)" #\;) ,nl))
1293    (else
1294     (pp indent (extern "\"C\""
1295                     int ,(s+ sysname "_dynamics")
1296                     "(double, const double*, double*, void*)" #\;) ,nl))
1297    )
1298
1299  (pp indent 
1300      (,(sprintf "class ~A:" sysname)
1301       "public Archiving_Node { "))
1302
1303  (pp indent+
1304       ("public: ")
1305       (,(s+ "~" sysname "();"))
1306       (,(s+ sysname "(const " sysname "&);"))
1307       (,(s+ sysname "();")))
1308
1309  (pp indent (#<<EOF
1310    using Node::connect_sender;
1311    using Node::handle;
1312
1313    port check_connection(Connection&, port);
1314   
1315    void handle(SpikeEvent &);
1316    void handle(CurrentEvent &);
1317    void handle(DataLoggingRequest &);
1318   
1319    port connect_sender(SpikeEvent &, port);
1320    port connect_sender(CurrentEvent &, port);
1321    port connect_sender(DataLoggingRequest &, port);
1322   
1323    void get_status(DictionaryDatum &) const;
1324    void set_status(const DictionaryDatum &);
1325   
1326    void init_node_(const Node& proto);
1327    void init_state_(const Node& proto);
1328    void init_buffers_();
1329    void calibrate();
1330   
1331    void update(Time const &, const long_t, const long_t);
1332EOF
1333))
1334
1335 
1336  (case method
1337    ((cvode) 
1338     (pp indent (friend int ,(s+ sysname "_dynamics")
1339                     "(double, const N_Vector, N_Vector, void*)" #\;) ,nl))
1340    (else
1341     (pp indent (friend int ,(s+ sysname "_dynamics")
1342                     "(double, const double*, double*, void*)" #\;) ,nl))
1343    )
1344
1345  (pp indent (,(sprintf #<<EOF
1346    // The next two classes need to be friends to access the State_ class/member
1347    friend class RecordablesMap<~A>;
1348    friend class UniversalDataLogger<~A>;
1349EOF
1350sysname sysname)))
1351
1352 
1353  (pp indent+ ,nl "struct Parameters_ { ")
1354
1355  (for-each
1356   (lambda (x) (pp indent+ (,(sprintf "double ~A;" x))))
1357   (append
1358    (map (compose nest-name first) const-defs)
1359    (map (compose nest-name first) defaults)))
1360
1361  (pp indent+ 
1362      ("Parameters_();")
1363      ("void get(DictionaryDatum&) const;")
1364      ("void set(const DictionaryDatum&);"))
1365
1366  (pp indent+ "};")
1367
1368  (pp indent+ ,nl "struct State_ { ")
1369
1370  (pp indent+ ,nl 
1371      "enum StateVecElems {"
1372      (,(slp ", " (map (lambda (x) 
1373                         (let ((n (string-upcase (->string (first x))))
1374                               (ni (second x)))
1375                           (s+ n " = " ni)))
1376                       state-index-map)))
1377      "};"
1378      "int_t     r_;"
1379
1380      (,(sprintf "double y_[~A];" (length state-index-map)))
1381
1382      ,(if (not (and (null? external-eq-defs) (null? pool-ions)))
1383          `(,(sprintf "double ~A;" 
1384                      (slp ", " 
1385                           (delete-duplicates
1386                            (append (map pool-ion-in pool-ions)
1387                                    (map pool-ion-out pool-ions)
1388                                    (map first i-eqs)
1389                                    ))
1390                           ))
1391            )
1392          `())
1393
1394      "State_(const Parameters_& p);" 
1395      "State_(const State_& s);"
1396      "State_& operator=(const State_& s);"
1397      "void get(DictionaryDatum&) const;"
1398      "void set(const DictionaryDatum&, const Parameters_&);")
1399
1400  (pp indent+ "};")
1401
1402
1403  (pp indent+ ,nl #<<EOF
1404struct Variables_ {
1405      int_t    RefractoryCounts_;
1406      double   U_old_; // for spike-detection
1407};
1408EOF
1409)
1410
1411  (pp indent+ ,nl "struct Buffers_ { ")
1412
1413  (pp indent+ ,nl 
1414      (,(sprintf "Buffers_(~A&);" sysname))
1415      (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
1416      (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))
1417
1418  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
1419    (for-each
1420     (lambda (psc)
1421       (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
1422     pscs))
1423
1424  (pp indent+ ,nl
1425  (,(case method
1426
1427    ((leapfrog)
1428#<<EOF
1429  rhsfn_t sys_; // pointer to function that computes dv/dt
1430  double *u, *v, *dvdt;  // intermediate state vectors used by leapfrog method
1431  unsigned int  N;  // size of state vector used by leapfrog method
1432EOF
1433)
1434
1435    ((cvode)
1436#<<EOF
1437  N_Vector y;    //!< current state vector used by CVode
1438  void *   sys_;  //!< CVode control structure
1439EOF
1440)
1441
1442    ((gsl)
1443#<<EOF
1444  gsl_odeiv_step*    s_;    //!< stepping function
1445  gsl_odeiv_control* c_;    //!< adaptive stepsize control function
1446  gsl_odeiv_evolve*  e_;    //!< evolution function
1447  gsl_odeiv_system   sys_;  //!< struct describing system
1448EOF
1449)
1450
1451   )))
1452
1453
1454  (pp indent+ ,nl
1455#<<EOF
1456  RingBuffer currents_;
1457
1458  double_t step_;           //!< step size in ms
1459  double   IntegrationStep_;//!< current integration time step, updated by solver
1460
1461  /** 
1462   * Input current injected by CurrentEvent.
1463   * This variable is used to transport the current applied into the
1464   * _dynamics function computing the derivative of the state vector.
1465   * It must be a part of Buffers_, since it is initialized once before
1466   * the first simulation, but not modified before later Simulate calls.
1467   */
1468  double_t I_stim_;
1469EOF
1470)
1471  (pp indent+ "};")
1472
1473  (pp indent+
1474      "template <State_::StateVecElems elem>"
1475      "double_t get_y_elem_() const { return S_.y_[elem]; }"
1476      "Parameters_ P_;"
1477      "State_      S_;"
1478      "Variables_  V_;"
1479      "Buffers_    B_;"
1480      (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
1481      "}; ")
1482
1483  (pp indent+ (,(sprintf #<<EOF
1484  inline
1485  port ~A::check_connection(Connection& c, port receptor_type)
1486  {
1487    SpikeEvent e;
1488    e.set_sender(*this);
1489    c.check_event(e);
1490    return c.get_target()->connect_sender(e, receptor_type);
1491  }
1492
1493  inline
1494  port ~A::connect_sender(SpikeEvent&, port receptor_type)
1495  {
1496    if (receptor_type != 0)
1497      throw UnknownReceptorType(receptor_type, get_name());
1498    return 0;
1499  }
1500 
1501  inline
1502  port ~A::connect_sender(CurrentEvent&, port receptor_type)
1503  {
1504    if (receptor_type != 0)
1505      throw UnknownReceptorType(receptor_type, get_name());
1506    return 0;
1507  }
1508
1509  inline
1510  port ~A::connect_sender(DataLoggingRequest& dlr, 
1511                                      port receptor_type)
1512  {
1513    if (receptor_type != 0)
1514      throw UnknownReceptorType(receptor_type, get_name());
1515    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
1516  }
1517
1518  inline
1519    void ~A::get_status(DictionaryDatum &d) const
1520  {
1521    P_.get(d);
1522    S_.get(d);
1523    Archiving_Node::get_status(d);
1524
1525    (*d)[names::recordables] = recordablesMap_.get_list();
1526
1527    def<double_t>(d, names::t_spike, get_spiketime_ms());
1528  }
1529
1530  inline
1531    void ~A::set_status(const DictionaryDatum &d)
1532  {
1533    Parameters_ ptmp = P_;  // temporary copy in case of errors
1534    ptmp.set(d);                       // throws if BadProperty
1535    State_      stmp = S_;  // temporary copy in case of errors
1536    stmp.set(d, ptmp);                 // throws if BadProperty
1537
1538    // We now know that (ptmp, stmp) are consistent. We do not 
1539    // write them back to (P_, S_) before we are also sure that 
1540    // the properties to be set in the parent class are internally 
1541    // consistent.
1542    Archiving_Node::set_status(d);
1543
1544    // if we get here, temporaries contain consistent set of properties
1545    P_ = ptmp;
1546    S_ = stmp;
1547
1548    calibrate();
1549  }
1550EOF
1551  sysname sysname sysname 
1552  sysname sysname sysname)))
1553
1554  (pp indent "}" ,nl)
1555  )
1556
1557(define (output-update
1558         sysname method state-index-map steady-state-index-map 
1559         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1560         reaction-eq-defs i-eqs pool-ions perm-ions 
1561         synapse-info
1562         indent indent+)
1563 
1564  (let* ((vi (lookup-def 'v state-index-map))
1565         (threshold (sprintf #<<EOF
1566            // (threshold && maximum)
1567            if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
1568              {
1569                S_.r_ = V_.RefractoryCounts_;
1570               
1571                set_spiketime(Time::step(origin.get_steps()+lag+1));
1572               
1573                SpikeEvent se;
1574                network()->send(*this, se, lag);
1575              }
1576EOF
1577        vi vi)))
1578   
1579    (pp indent  (,(sprintf #<<EOF
1580void ~A::update(Time const & origin, const long_t from, const long_t to)
1581EOF
1582  sysname)))
1583
1584    (pp indent  #\{)
1585
1586    (pp indent+ (,(sprintf #<<EOF
1587assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
1588    assert(from < to);
1589
1590    double tout;
1591    long_t current_steps = origin.get_steps();
1592
1593    for ( long_t lag = from ; lag < to ; ++lag )
1594      {
1595        double h = B_.step_;   
1596        double tt = 0.0 ;
1597        V_.U_old_ = S_.y_[~A];
1598EOF
1599vi)))
1600
1601    (case method
1602      ((leapfrog)
1603       (pp indent+ (,(sprintf #<<EOF
1604
1605        tt = (current_steps+lag)*h;
1606        tout = (current_steps+lag+1)*h;
1607
1608        while (tt < tout)
1609        {
1610           B_.sys_ (tt, S_.y_, B_.dvdt, reinterpret_cast<void*>(this));
1611
1612           for (int i = 0; i < B_.N; i++)
1613           {
1614              B_.u[i] = S_.y_[i] + 0.5 * h * B_.dvdt[i];
1615           }
1616
1617           B_.sys_ (tt+0.5*h, B_.u, B_.dvdt, reinterpret_cast<void*>(this));
1618
1619           for (int i = 0; i < B_.N; i++)
1620           { 
1621              B_.v[i] = S_.y_[i] + h * B_.dvdt[i];
1622           }
1623
1624           tt += h;
1625        }
1626
1627        for (int i = 0; i < B_.N; i++)
1628        {
1629           S_.y_[i] = B_.v[i];
1630        }
1631
1632EOF
1633 ))
1634))
1635      ((cvode)
1636       (pp indent+ (,(sprintf #<<EOF
1637
1638        int N = NV_LENGTH_S (B_.y);
1639        tout = (current_steps+lag+1)*h;
1640
1641        // adaptive step integration
1642        while (tt < tout)
1643        {
1644          const int status = CVode (B_.sys_, tout, B_.y, &tt, CV_NORMAL);
1645
1646          if ( status != CV_SUCCESS )
1647            throw GSLSolverFailure(get_name(), status);
1648        }
1649
1650        for (int i = 0; i < N; i++)
1651        {
1652           S_.y_[i] = Ith(B_.y,i);
1653        }
1654EOF
1655 ))
1656))
1657      ((gsl)
1658       (pp indent+ (,(sprintf #<<EOF
1659                             
1660        // adaptive step integration
1661        while (tt < h)
1662        {
1663          const int status = gsl_odeiv_evolve_apply(B_.e_, B_.c_, B_.s_, 
1664                                 &B_.sys_,              // system of ODE
1665                                 &tt,                   // from t...
1666                                  h,             // ...to t=t+h
1667                                 &B_.IntegrationStep_ , // integration window (written on!)
1668                                  S_.y_);               // neuron state
1669
1670          if ( status != GSL_SUCCESS )
1671            throw GSLSolverFailure(get_name(), status);
1672        }
1673EOF
1674))))
1675      )
1676
1677    (let ((isyns (lookup-def 'i-synapses synapse-info))
1678          (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
1679
1680      (for-each (lambda (isyn psc)
1681                  (let ((gi (lookup-def (nest-name (fourth isyn)) state-index-map)))
1682                    (if gi
1683                        (pp indent+ (,(sprintf "      S_.y_[~A] += B_.spike_~A.get_value(lag);"
1684                                               gi (second psc)))
1685                        ))))
1686                isyns pscs))
1687
1688    (pp indent+ (,(sprintf #<<EOF
1689        // sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
1690        // refractory?
1691        if (S_.r_)
1692          {
1693            --S_.r_;
1694          }
1695        else
1696          {
1697           ~A
1698          }
1699   
1700        // set new input current
1701        B_.I_stim_ = B_.currents_.get_value(lag);
1702
1703        // log state data
1704        B_.logger_.record_data(current_steps + lag);
1705
1706      }
1707EOF
1708  (if (lookup-def 'V_t const-defs) threshold "")
1709  )))
1710       
1711  (pp indent  #\})
1712
1713))
1714
1715(define (output-buffers+nodes
1716         sysname method state-index-map steady-state-index-map 
1717         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1718         reaction-eq-defs i-eqs pool-ions perm-ions 
1719         synapse-info
1720         indent indent+)
1721
1722  (let ((N (length state-index-map)))
1723
1724    (case method
1725
1726
1727      ((leapfrog)
1728
1729       (pp indent ,nl
1730           (,(sprintf #<<EOF
1731~A::Buffers_::Buffers_(~A& n)
1732    : logger_(n),
1733      sys_(0),
1734      N(0),
1735      u(0),
1736      v(0),
1737      dvdt(0)
1738{
1739    // Initialization of the remaining members is deferred to
1740    // init_buffers_().
1741}
1742EOF
1743  sysname sysname))
1744
1745  (,(sprintf #<<EOF
1746~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
1747    : logger_(n),
1748      sys_(0),
1749      N(0),
1750      u(0),
1751      v(0),
1752      dvdt(0)
1753{
1754    // Initialization of the remaining members is deferred to
1755    // init_buffers_().
1756}
1757EOF
1758  sysname sysname))
1759    ))
1760
1761      ((cvode)
1762
1763       (pp indent ,nl
1764           (,(sprintf #<<EOF
1765~A::Buffers_::Buffers_(~A& n)
1766    : logger_(n),
1767      sys_(0)
1768{
1769    // Initialization of the remaining members is deferred to
1770    // init_buffers_().
1771}
1772EOF
1773  sysname sysname))
1774
1775  (,(sprintf #<<EOF
1776~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
1777    : logger_(n),
1778      sys_(0)
1779{
1780    // Initialization of the remaining members is deferred to
1781    // init_buffers_().
1782}
1783EOF
1784  sysname sysname))
1785    ))
1786
1787    ((gsl)
1788
1789     (pp indent ,nl
1790         (,(sprintf #<<EOF
1791~A::Buffers_::Buffers_(~A& n)
1792    : logger_(n),
1793      s_(0),
1794      c_(0),
1795      e_(0)
1796{
1797    // Initialization of the remaining members is deferred to
1798    // init_buffers_().
1799}
1800EOF
1801  sysname sysname))
1802
1803  (,(sprintf #<<EOF
1804~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
1805    : logger_(n),
1806      s_(0),
1807      c_(0),
1808      e_(0)
1809{
1810    // Initialization of the remaining members is deferred to
1811    // init_buffers_().
1812}
1813EOF
1814  sysname sysname))
1815    ))
1816
1817
1818   )
1819
1820
1821
1822  (pp indent  ,nl  (,(sprintf #<<EOF
1823~A::~A()
1824    : Archiving_Node(), 
1825      P_(), 
1826      S_(P_),
1827      B_(*this)
1828{
1829    recordablesMap_.create();
1830}
1831EOF
1832  sysname sysname)))
1833
1834  (pp indent  ,nl  (,(sprintf #<<EOF
1835~A::~A(const ~A& n)
1836    : Archiving_Node(n), 
1837      P_(n.P_), 
1838      S_(n.S_),
1839      B_(n.B_, *this)
1840{
1841}
1842EOF
1843  sysname sysname sysname)))
1844
1845
1846
1847  (pp indent (,(s+ sysname "::~" sysname "()")))
1848
1849  (case method
1850
1851   ((leapfrog)
1852    (pp indent #<<EOF
1853{
1854    if ( B_.u ) free(B_.u);
1855    if ( B_.v ) free(B_.v);
1856    if ( B_.dvdt ) free(B_.dvdt);
1857}
1858EOF
1859  ))
1860
1861
1862   ((gsl)
1863    (pp indent #<<EOF
1864{
1865    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
1866    if ( B_.s_ ) gsl_odeiv_step_free(B_.s_);
1867    if ( B_.c_ ) gsl_odeiv_control_free(B_.c_);
1868    if ( B_.e_ ) gsl_odeiv_evolve_free(B_.e_);
1869}
1870EOF
1871  ))
1872
1873   ((cvode)
1874    (pp indent #<<EOF
1875{
1876
1877  if ( B_.sys_ )
1878  {
1879    /* Free y vector */
1880    N_VDestroy_Serial(B_.y);
1881
1882    /* Free integrator memory */
1883//    CVodeFree(&B_.sys_);
1884  }
1885}
1886EOF
1887    ))
1888   )
1889
1890
1891  (pp indent ,nl (,(sprintf #<<EOF
1892  void ~A::init_node_(const Node& proto)
1893{
1894    const ~A& pr = downcast<~A>(proto);
1895    P_ = pr.P_;
1896    S_ = State_(P_);
1897}
1898EOF
1899  sysname sysname sysname)))
1900
1901  (pp indent  ,nl (,(sprintf #<<EOF
1902void ~A::init_state_(const Node& proto)
1903{
1904    const ~A& pr = downcast<~A>(proto);
1905    S_ = State_(pr.P_);
1906}
1907EOF
1908  sysname sysname sysname)))
1909
1910
1911  (pp indent ,nl (,(sprintf #<<EOF
1912void ~A::init_buffers_()
1913{
1914EOF
1915 sysname)))
1916
1917 (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
1918   (for-each (lambda (psc)
1919               (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
1920       pscs))
1921
1922       (pp indent+ (#<<EOF
1923    B_.currents_.clear();           
1924    Archiving_Node::clear_history();
1925
1926    B_.logger_.reset();
1927
1928    B_.step_ = Time::get_resolution().get_ms();
1929    B_.IntegrationStep_ = B_.step_;
1930
1931    B_.I_stim_ = 0.0;
1932EOF
1933))
1934
1935
1936    (case method
1937     ((leapfrog)
1938
1939      (pp indent+ (,(sprintf #<<EOF
1940
1941    B_.N = ~A;
1942    B_.sys_ = ~A_dynamics;
1943    B_.u = (double *)malloc(sizeof(double) * B_.N);
1944    assert (B_.u);
1945    B_.v = (double *)malloc(sizeof(double) * B_.N);
1946    assert (B_.v);
1947    B_.dvdt = (double *)malloc(sizeof(double) * B_.N);
1948    assert (B_.dvdt);
1949
1950EOF
1951N sysname))))
1952
1953     ((cvode)
1954
1955      (pp indent+ (,(sprintf #<<EOF
1956
1957    int status, N;
1958
1959    N = ~A;
1960
1961    /* Creates serial vector of length N */
1962    B_.y = N_VNew_Serial(N);
1963    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) abort();
1964
1965    for (int i = 0; i < N; i++)
1966    {
1967       Ith(B_.y,i) = S_.y_[i];
1968    }
1969 
1970    /* Calls CVodeCreate to create the solver memory and specify the 
1971     * Backward Differentiation Formula and the use of a Newton iteration */
1972    B_.sys_ = CVodeCreate(CV_BDF, CV_NEWTON);
1973    if (check_flag((void *)B_.sys_, "CVodeCreate", 0)) abort();
1974
1975   /* Calls CVodeInit to initialize the integrator memory and specify the
1976    * right hand side function in y'=f(t,y), the initial time, and
1977    * the initial values. */
1978    status = CVodeInit (B_.sys_, ~A_dynamics, 0.0, B_.y);
1979    if (check_flag(&status, "CVodeInit", 1)) abort();
1980
1981    /* Sets the relative and absolute error tolerances of CVode  */
1982    status = CVodeSStolerances (B_.sys_, 1e-4, 1e-4);
1983    if (check_flag(&status, "CVodeSStolerances", 1)) abort();
1984
1985    /* Turn on CVode stability limit detection (only applicable for order 3 and greater) */
1986    status = CVodeSetStabLimDet (B_.sys_,true);
1987    if (check_flag(&status, "CVodeSetStabLimDet", 1)) abort();
1988
1989    /* Sets the maximum order of CVode  */
1990    status = CVodeSetMaxOrd (B_.sys_,1);
1991    if (check_flag(&status, "CVodeSetMaxOrd", 1)) abort();
1992
1993    /* Sets minimum and maximum step size. */
1994    //    status = CVodeSetMinStep (B_.sys_,1e-6);
1995    //    if (check_flag(&status, "CVodeSetMinStep", 1)) abort();
1996    status = CVodeSetMaxStep (B_.sys_,B_.step_);
1997    if (check_flag(&status, "CVodeSetMaxStep", 1)) abort();
1998
1999   /* Calls CVodeSetUserData to configure the integrator to pass the 
2000    * params structure to the right-hand function */
2001   status = CVodeSetUserData(B_.sys_, reinterpret_cast<void*>(this));
2002   if (check_flag(&status, "CVodeSetUserData", 1)) abort();
2003
2004    /* Initializes dense linear solver. */
2005    status = CVDense (B_.sys_, N);
2006    if (check_flag(&status, "CVDense", 1)) abort();
2007EOF
2008    N sysname))))
2009
2010    ((gsl)
2011
2012    (pp indent+ (,(sprintf #<<EOF
2013
2014    int N = ~A;
2015    static const gsl_odeiv_step_type* T1 = gsl_odeiv_step_rkf45;
2016 
2017    if ( B_.s_ == 0 )
2018      B_.s_ = gsl_odeiv_step_alloc (T1, N);
2019    else
2020      gsl_odeiv_step_reset(B_.s_);
2021   
2022    if ( B_.c_ == 0 ) 
2023      B_.c_ = gsl_odeiv_control_y_new (1e-3, 0.0);
2024    else
2025      gsl_odeiv_control_init(B_.c_, 1e-3, 0.0, 1.0, 0.0);
2026   
2027    if ( B_.e_ == 0 ) 
2028      B_.e_ = gsl_odeiv_evolve_alloc(N);
2029    else
2030      gsl_odeiv_evolve_reset(B_.e_);
2031 
2032    B_.sys_.function  = ~A_dynamics;
2033    B_.sys_.jacobian  = 0;
2034    B_.sys_.dimension = N;
2035    B_.sys_.params    = reinterpret_cast<void*>(this);
2036EOF
2037  N sysname))))
2038  )
2039  (pp indent ,nl #\})
2040
2041  (let ((iv (lookup-def 'v state-index-map)))
2042    (if iv
2043        (pp indent ,nl (,(sprintf #<<EOF
2044void ~A::calibrate()
2045{
2046    B_.logger_.init()
2047    V_.RefractoryCounts_ = 20;
2048    V_.U_old_ = S_.y_[~A];
2049}
2050EOF
2051  sysname iv)))
2052  ))
2053 ))
2054
2055
2056
2057(define (nemo:nest-translator sys . rest)
2058
2059  (define (cid x)  (second x))
2060  (define (cn x)   (first x))
2061
2062  (let-optionals rest ((dirname ".") (method #f))
2063
2064    (let ((method (or method 'gsl)))
2065
2066      (if (not (member method '(gsl cvode leapfrog)))
2067          (nemo:error 'nemo:nest-translator ": unknown method " method))
2068
2069  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
2070    (let ((imports  ((dis 'imports)  sys))
2071          (exports  ((dis 'exports)  sys)))
2072      (let* ((indent      0)
2073             (indent+     (+ 2 indent ))
2074
2075             (sysname     (nest-name ((dis 'sysname) sys)))
2076             (prefix      (->string sysname))
2077             (deps*       ((dis 'depgraph*)  sys))
2078             (consts      ((dis 'consts)     sys))
2079             (asgns       ((dis 'asgns)      sys))
2080             (states      ((dis 'states)     sys))
2081             (reactions   ((dis 'reactions)  sys))
2082             (defuns      ((dis 'defuns)     sys))
2083             (components  ((dis 'components) sys))
2084             
2085             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
2086             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
2087
2088             (const-defs       (filter-map
2089                                (lambda (nv)
2090                                  (and (not (member (first nv) nest-builtin-consts))
2091                                       (let ((v1 (canonicalize-expr/C++ (second nv))))
2092                                         (list (nest-name (first nv)) v1))))
2093                                consts))
2094             
2095             (gate-complex-info    (nemo:gate-complex-query sys))
2096             (synapse-info         (nemo:post-synaptic-conductance-query sys))
2097             (i-synapses           (lookup-def 'i-synapses synapse-info))
2098             (defaults             (nemo:defaults-query sys))
2099
2100             (gate-complexes  (lookup-def 'gate-complexes gate-complex-info))
2101             (perm-ions       (map (match-lambda ((comp i e erev val) `(,comp ,(nest-name i) ,(nest-name e) ,erev)))
2102                                   (lookup-def 'perm-ions gate-complex-info)))
2103             (acc-ions        (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out)))))
2104                                   (lookup-def 'acc-ions gate-complex-info)))
2105             (epools          (lookup-def 'pool-ions gate-complex-info))
2106             (pool-ions       (pool-ion-name-map nest-name  epools))
2107
2108             (i-gates         (lookup-def 'i-gates gate-complex-info))
2109
2110             (comprc         (any (match-lambda ((name 'membrane-tau id) (list name id)) (else #f)) components))
2111             (mrc            (and comprc (car ((dis 'component-exports) sys (cid comprc)))))
2112               
2113             (areacomp        (any (match-lambda ((name 'membrane-area id) (list name id)) (else #f)) components))
2114             (marea           (and areacomp (car ((dis 'component-exports) sys (cid areacomp)))))
2115               
2116             (i-eqs (filter-map
2117                     (lambda (gate-complex) 
2118                       
2119                       (let* ((label             (first gate-complex))
2120                              (n                 (second gate-complex))
2121                              (subcomps          ((dis 'component-subcomps) sys n))
2122                              (acc               (lookup-def 'accumulating-substance subcomps))
2123                              (perm              (lookup-def 'permeating-ion subcomps))
2124                              (permqs            (and perm ((dis 'component-exports) sys (cid perm))))
2125                              (pore              (lookup-def 'pore subcomps))
2126                              (permeability      (lookup-def 'permeability subcomps))
2127                              (gates             (filter (lambda (x) (equal? (car x) 'gate)) subcomps))
2128                              (sts               (map (lambda (gate) ((dis 'component-exports) sys (cid gate))) gates))
2129                              )
2130
2131
2132                         (if (and pore (null? permqs))
2133                             (nemo:error 'nemo:nest-translator ": ion channel definition " label
2134                                         "permeating-ion component lacks exported quantities"))
2135
2136                         (for-each
2137                          (lambda (st)
2138                            (if (null? st)
2139                                (nemo:error 'nemo:nest-translator: "ion channel definition " label
2140                                            "gate component lacks exported quantities")))
2141                          sts)
2142                         
2143                         (if (not (or pore permeability))
2144                             (nemo:error 'nemo:nest-translator ": ion channel definition " label
2145                                         "lacks any pore or permeability components"))
2146
2147
2148                         (cond ((and perm permeability (pair? gates))
2149                                     (let* ((i     (nest-name (s+ 'i (cn perm))))
2150                                            (pmax  (car ((dis 'component-exports) sys (cid permeability))))
2151                                            (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
2152                                            (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
2153                                            (gion  `(* ,pmax ,(sum (map (lambda (gpwr) 
2154                                                                          (match gpwr ((x)  x) (else `(* ,@gpwr))))
2155                                                                        gpwrs))))
2156                                            )
2157                                         (list i #f gion (nest-name (s+ 'i_ label) ))))
2158
2159                               ((and perm pore (pair? gates))
2160                                (case (cn perm)
2161                                  ((non-specific)
2162                                   (let* ((i     (nest-name 'i))
2163                                          (e     (car permqs))
2164                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
2165                                          (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
2166                                          (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
2167                                          (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
2168                                                                        (match gpwr ((x)  x) (else `(* ,@gpwr))))
2169                                                                      gpwrs))))
2170                                          )
2171                                     (list i e gion  (nest-name (s+ 'i_ label) ))))
2172
2173                                  (else
2174                                   (let* ((i     (nest-name (s+ 'i (cn perm))))
2175                                          (e     (car permqs))
2176                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
2177                                          (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
2178                                          (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
2179                                          (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
2180                                                                        (match gpwr ((x)  x) (else `(* ,@gpwr))))
2181                                                                      gpwrs))))
2182                                          )
2183                                     (list i e gion (nest-name (s+ 'i_ label) ))))))
2184                               
2185                               ((and perm pore)
2186                                (case (cn perm)
2187                                  ((non-specific)
2188                                   (let* ((i     (nest-name 'i))
2189                                          (e     (car permqs))
2190                                          (gmax  (car ((dis 'component-exports) sys (cid pore)))))
2191                                     (list i e gmax (nest-name (s+ 'i_ label) ))))
2192                                  (else
2193                                   (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))))
2194                               
2195                               ((and acc pore (pair? gates))
2196                                (let* ((i     (nest-name (s+ 'i (cn acc))))
2197                                       (gmax  (car ((dis 'component-exports) sys (cid pore))))
2198                                       (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
2199                                       (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
2200                                       (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
2201                                                                     (match gpwr ((x)  x) (else `(* ,@gpwr))))
2202                                                                   gpwrs))))
2203                                       )
2204                                  (list i #f gion  (nest-name (s+ 'i_ label) ))))
2205                               (else (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))
2206                               )))
2207                     gate-complexes))
2208
2209             (i-names (delete-duplicates (map first i-eqs)))
2210               
2211             (i-eqs  (fold  (lambda (i-gate ax) 
2212                              (let ((i-gate-var (first i-gate)))
2213                                (cons (list (nest-name 'i) #f i-gate-var (s+ 'i_ (second i-gate))) ax)))
2214                            i-eqs i-gates))
2215
2216             (i-eqs  (fold  (lambda (i-synapse ax) 
2217                              (cons (list (first i-synapse) (third i-synapse) (fourth i-synapse) (gensym 'i_syn )) ax))
2218                            i-eqs 
2219                            i-synapses))
2220
2221             (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
2222             
2223             (i-eqs  (fold (lambda (b ax) 
2224                             (match b 
2225                                    ((and ps ((i e gion ii) . rst)) 
2226                                     (let loop ((ps ps) (summands (list)) (eqs (list)))
2227                                       (if (null? ps)
2228                                           
2229                                           (let* ((sum0  (sum summands))
2230                                                  (sum1  (rhsexpr/C++ sum0))
2231                                                  (sum2  (add-params-to-fncall 
2232                                                          (canonicalize-expr/C++ sum1))))
2233                                             (append eqs (list (list i sum2)) ax))
2234                                           
2235                                           (match-let (((i e gion ii) (car ps)))
2236                                                      (loop (cdr ps) 
2237                                                            (cons ii summands) 
2238                                                            (let* ((expr0 (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
2239                                                                   (expr1 (canonicalize-expr/C++ expr0)))
2240                                                              (cons (list ii expr1) eqs)))))))
2241                                   
2242                                    ((i e gion ii)
2243                                     (let* ((expr0  (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
2244                                            (expr1  (canonicalize-expr/C++ expr0)))
2245                                       (cons (list i expr1) ax)))
2246                                   
2247                                    (else ax)))
2248                           (list) i-bkts))
2249
2250             (external-eq-defs   (external-eq-defs sys))
2251
2252             (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
2253             
2254             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys)))
2255             
2256             (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
2257             
2258             (init-eq-defs       (poset->init-defs poset sys))
2259             
2260             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
2261                                      (poset->state-conserve-eq-defs poset sys)))
2262             
2263             (imports-sans-v (filter (lambda (x) (not (equal? 'v (first x)))) imports))
2264
2265             (v-eq    (cond ((and mcap mres marea (lookup-def 'v imports))
2266                             (list 'v (rhsexpr/C++ `(+ "(node.B_.I_stim_)" 
2267                                                       (/ (* (neg ,(sum i-names)) ,marea) (* 2. ,mres ,mcap))))))
2268                            ((and marea (lookup-def 'v imports))
2269                             (list 'v (rhsexpr/C++ `(+ "(node.B_.I_stim_)" 
2270                                                       (* (neg ,(sum i-names)) ,marea)))))
2271                            (else (list 'v 0.0))))
2272
2273             (v-eq    (and (member 'v globals) (not (null? i-names))
2274                           (let ((istim "(node.B_.I_stim_)" )) 
2275                             (cond ((and mrc marea)
2276                                    (list 'v (rhsexpr/C++ `(/ (+ ,istim (* (neg ,(sum i-names)) ,marea)) ,mrc))))
2277                                   (marea
2278                                    (list 'v (rhsexpr/C++ `(+ istim (* (neg ,(sum i-names)) ,marea)))))
2279                                   (mrc
2280                                    (list 'v (rhsexpr/C++ `(/ (+ istim (neg ,(sum i-names))) ,mrc))))
2281                                   (else
2282                                    (list 'v (rhsexpr/C++ `(+istim (neg ,(sum i-names))))))
2283                                   ))))
2284
2285             
2286             (state-index-map  (let ((acc (fold (lambda (def ax)
2287                                                  (let ((st-name (first def)))
2288                                                    (list (+ 1 (first ax)) 
2289                                                          (cons `(,st-name ,(first ax)) (second ax)))))
2290                                                (list 0 (list)) 
2291                                                (if (lookup-def 'v imports)
2292                                                    (cons (list 'v) rate-eq-defs)
2293                                                    rate-eq-defs)
2294                                                )))
2295                                 (second acc)))
2296             
2297             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
2298                                                         (let ((st-name (first def)))
2299                                                           (if (not (alist-ref st-name init-eq-defs))
2300                                                               (list (+ 1 (first ax)) 
2301                                                                     (cons `(,st-name ,(first ax)) (second ax)))
2302                                                               ax)))
2303                                                       (list 0 (list)) 
2304                                                       rate-eq-defs)))
2305                                        (second acc)))
2306             
2307             (dfenv 
2308              (map (lambda (x) (let ((n (first x))) (list n (nest-name (s+ "d_" n )))))
2309                   defuns))
2310
2311             )
2312       
2313       
2314       
2315        (for-each
2316         (lambda (a)
2317           (let ((acc-ion   (car a)))
2318             (if (assoc acc-ion perm-ions)
2319                 (nemo:error 'nemo:nest-translator 
2320                             ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
2321         acc-ions)
2322
2323
2324        (let ((cpp-output  (open-output-file (make-output-fname dirname prefix ".cpp")))
2325              (hpp-output  (open-output-file (make-output-fname dirname prefix ".h"))))
2326         
2327          ;; include files and other prelude definitions
2328          (with-output-to-port cpp-output
2329            (lambda ()
2330              (output-prelude sysname steady-state-index-map indent)
2331              ))
2332         
2333          (with-output-to-port cpp-output
2334            (lambda () (pp indent "namespace nest {" ,nl)))
2335         
2336          ;; user-defined functions
2337          (let ((define-fn  (make-define-fn sysname))
2338                (define-fn-header (make-define-fn-header sysname)))
2339           
2340            (for-each (lambda (fndef) 
2341                        (if (not (member (car fndef) builtin-fns))
2342                            (with-output-to-port cpp-output
2343                              (lambda ()
2344                                (apply define-fn-header (cons indent fndef))
2345                                (pp indent ,nl)))
2346                            ))
2347                      defuns)
2348
2349            (for-each (lambda (fndef) 
2350                        (if (not (member (car fndef) builtin-fns))
2351                            (with-output-to-port cpp-output
2352                              (lambda ()
2353                                (apply define-fn (cons indent fndef))
2354                                (pp indent ,nl)))
2355                            ))
2356                      defuns)
2357            )
2358         
2359         
2360          ;; derivative function
2361          (with-output-to-port cpp-output
2362            (lambda ()
2363              (output-dy sysname method imports-sans-v const-defs state-index-map 
2364                         external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
2365                         pool-ions mcap i-eqs v-eq
2366                         indent indent+)
2367              ))
2368
2369          ;;  system recordables
2370          (with-output-to-port cpp-output
2371            (lambda ()
2372              (output-recordables
2373               sysname state-index-map steady-state-index-map 
2374               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2375               reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
2376              (pp indent ,nl)
2377              ))
2378         
2379          ;; system parameters
2380          (with-output-to-port cpp-output
2381            (lambda ()
2382              (output-parameters sysname imports-sans-v const-defs 
2383                                 external-eq-defs defaults indent indent+)
2384              ))
2385         
2386          ;; initial values function
2387          (with-output-to-port cpp-output
2388            (lambda ()
2389              (output-init sysname state-index-map steady-state-index-map 
2390                              imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2391                              reaction-eq-defs i-eqs pool-ions perm-ions defaults indent indent+)
2392              (pp indent ,nl)
2393              ))
2394         
2395          ;; accessors and modifiers for states and parameters
2396          (with-output-to-port cpp-output
2397            (lambda ()
2398              (output-accessors+modifiers
2399               sysname imports-sans-v state-index-map 
2400               const-defs asgn-eq-defs rate-eq-defs 
2401               reaction-eq-defs i-eqs pool-ions perm-ions 
2402               indent indent+)
2403              (pp indent ,nl)
2404              ))
2405         
2406         
2407          ;; buffer and node initialization
2408          (with-output-to-port cpp-output
2409            (lambda ()
2410              (output-buffers+nodes
2411               sysname method state-index-map steady-state-index-map 
2412               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2413               reaction-eq-defs i-eqs pool-ions perm-ions 
2414               synapse-info
2415               indent indent+)
2416              (pp indent ,nl)
2417              ))
2418         
2419          ;; state update
2420          (with-output-to-port cpp-output
2421            (lambda ()
2422              (output-update
2423               sysname method state-index-map steady-state-index-map 
2424               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2425               reaction-eq-defs i-eqs pool-ions perm-ions 
2426               synapse-info indent indent+)
2427              (pp indent ,nl)
2428              ))
2429         
2430          ;; spike handler
2431          (with-output-to-port cpp-output
2432            (lambda ()
2433              (output-handle
2434               sysname state-index-map steady-state-index-map 
2435               const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2436               reaction-eq-defs i-eqs pool-ions perm-ions 
2437               synapse-info
2438               indent indent+)
2439              (pp indent ,nl)
2440              ))
2441         
2442          (with-output-to-port cpp-output
2443            (lambda () (pp indent "}" ,nl)))
2444
2445
2446             (with-output-to-port hpp-output
2447               (lambda ()
2448                 (output-header
2449                  sysname method state-index-map steady-state-index-map defaults
2450                  external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
2451                  reaction-eq-defs i-eqs pool-ions perm-ions 
2452                  synapse-info
2453                  indent indent+)
2454                 (pp indent ,nl)
2455                 ))
2456             
2457             (close-output-port cpp-output)
2458             (close-output-port hpp-output)
2459           
2460           ))
2461      ))
2462  ))
2463)
2464)
Note: See TracBrowser for help on using the repository browser.