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

Last change on this file since 27093 was 27093, checked in by Ivan Raikov, 9 years ago

nemo: eliminated dependency on environments

File size: 50.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-2012 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)
32
33(define nest-builtin-consts
34  `())
35
36(define C++-ops
37  `(+ - * / > < <= >= =))
38
39(define builtin-fns
40  `(+ - * / pow neg abs atan asin acos sin cos exp ln
41      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
42      > < <= >= = and or round ceiling floor max min
43      ))
44
45(define (nest-name s)
46  (let ((cs (string->list (->string s))))
47    (let loop ((lst (list)) (cs cs))
48      (if (null? cs) (string->symbol (list->string (reverse lst)))
49          (let* ((c (car cs))
50                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
51                           (else #\_))))
52            (loop (cons c1 lst) (cdr cs)))))))
53           
54(define (rhsexpr/C++ expr)
55  (match expr 
56         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/C++ x)) es)))
57         (('pow x y)  (if (and (integer? y)  (positive? y))
58                          (if (> y 1)  (let ((tmp (gensym "x")))
59                                         `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
60                              x)
61                            (if (and (number? y) (zero? y)) 1.0 expr)))
62         ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (nest-name s)) (map (lambda (x) (rhsexpr/C++ x)) es)) expr))
63         (id          (if (symbol? id) (nest-name id) id))))
64
65(define (nest-state-name n s)
66  (nest-name (s+ n s)))
67
68(define-syntax pp
69  (syntax-rules ()
70    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
71
72
73(define group/C++   (doc:block 2 (doc:text "(") (doc:text ")")))
74(define block/C++   (doc:block 2 (doc:text "{") (doc:text "}")))
75(define (stmt/C++ x) 
76  (match x
77         (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
78         (else  (doc:cons x (doc:text ";")))))
79
80
81(define (ifthen/C++ c e1 e2)
82  (doc:nest 2
83    (doc:connect (doc:group (doc:connect (doc:text "if") c))
84                 (doc:connect (doc:nest 2 e1)
85                              (doc:nest 2 (doc:connect 
86                                           (doc:text "else") 
87                                           e2))))
88    ))
89
90(define (letblk/C++ e1 e2)
91  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
92        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
93        (else (doc:connect (doc:group (doc:nest 2 (stmt/C++ e1)))
94                           (doc:group (doc:nest 2 e2))))))
95       
96
97(define (format-op/C++ indent op args)
98  (let ((op1 (doc:text (->string op))))
99    (if (null? args) op1
100        (match args
101               ((x)      (doc:concat (list op1 x)))
102               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
103               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
104               (lst      (let* ((n   (length lst))
105                                (n/2 (inexact->exact (round (/ n 2)))))
106                           (doc:concat 
107                            (intersperse 
108                             (list (format-op/C++ indent op (take lst n/2 )) op1 
109                                   (format-op/C++ indent op (drop lst n/2 )))
110                             (doc:space)))))))))
111
112(define (format-fncall/C++ indent op args)
113  (let ((op1 (doc:text (->string op))))
114    (doc:cons op1 (group/C++ ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
115
116(define (name-normalize expr)
117  (match expr 
118         (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
119         (('let bs e)
120          `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
121         ((f . es) 
122          (cons f (map name-normalize es)))
123         ((? symbol? ) (nest-name expr))
124         ((? atom? ) expr)))
125
126
127(define (canonicalize-expr/C++ expr)
128  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term)))
129    (let* ((expr1 (if-convert expr))
130           (expr2 (subst-convert expr1 subst-empty))
131           (expr3 (let-lift expr2))
132           (expr4 (name-normalize expr3)))
133      expr4)))
134
135
136(define (format-expr/C++ indent expr . rest) 
137  (let-optionals rest ((rv #f))
138   (let ((indent+ (+ 2 indent)))
139    (match expr
140       (('let bindings body)
141        (letblk/C++
142         (fold-right 
143           (lambda (x ax)
144             (letblk/C++
145              (match (second x)
146                     (('if c t e)
147                      (ifthen/C++
148                       (group/C++ (format-expr/C++ indent c))
149                       (block/C++ (format-expr/C++ indent t (first x)))
150                       (block/C++ (format-expr/C++ indent e (first x)))))
151                     (else
152                      (stmt/C++
153                       (format-op/C++ indent+ " = "
154                                         (list (format-expr/C++ indent (first x) )
155                                               (format-expr/C++ indent (second x)))))))
156              ax))
157           (doc:empty) bindings)
158         (match body
159                (('let _ _) (format-expr/C++ indent body rv))
160                (else
161                 (let ((body1 (doc:nest indent (format-expr/C++ indent body))))
162                   (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) body1)))
163                       body1))))))
164       
165       (('if . rest) (error 'format-expr/C++ "invalid if statement " expr))
166
167       ((op . rest) 
168       (let ((op (case op ((ln) 'log) (else op))))
169         (let ((fe
170                (if (member op C++-ops)
171                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
172                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
173                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
174                      (case op
175                        ((/) 
176                         (format-op/C++ indent op 
177                                          (map (lambda (x) 
178                                                 (let ((fx (format-expr/C++ indent+ x)))
179                                                   (if (or (symbol? x) (number? x)) fx
180                                                       (if (or mul? plmin?) (group/C++ fx) fx)))) rest)))
181                        ((*) 
182                         (format-op/C++ indent op 
183                                          (map (lambda (x) 
184                                                 (let ((fx (format-expr/C++ indent+ x)))
185                                                   (if (or (symbol? x) (number? x)) fx
186                                                       (if plmin? (group/C++ fx) fx)))) rest)))
187                       
188                        (else
189                         (format-op/C++ indent op 
190                                          (map (lambda (x) 
191                                                 (let ((fx (format-expr/C++ indent+ x))) fx)) rest)))))
192                   
193                    (let ((op (case op ((neg) '-) (else op))))
194                      (format-fncall/C++ indent op (map (lambda (x) (format-expr/C++ indent+ x)) rest))))))
195           (if rv 
196               (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
197               fe))))
198     
199      (else  (let ((fe (doc:text (->string expr))))
200               (if rv 
201                   (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
202                   fe)))))))
203               
204         
205(define (expr->string/C++ x . rest)
206  (let-optionals rest ((rv #f) (width 72))
207    (sdoc->string (doc:format width (format-expr/C++ 2 x rv)))))
208
209
210
211(define (state-init n init)
212  (let* ((init  (rhsexpr/C++ init))
213         (init1 (canonicalize-expr/C++ init)))
214    (list (nest-name n) init1)))
215
216
217(define (asgn-eq n rhs)
218  (let* ((fbody   (rhsexpr/C++ rhs))
219         (fbody1  (canonicalize-expr/C++ fbody)))
220    (list (nest-name n) fbody1)))
221
222
223(define (reaction-eq n open transitions conserve)
224  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nest-state-name)))
225    (let ((nodes ((g 'nodes))))
226      (list (nest-name n) (nest-state-name n open)))))
227
228
229(define (reaction-transition-eqs n initial open transitions conserve power)
230  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nest-state-name)))
231     (let* ((out-edges  (g 'out-edges))
232            (in-edges   (g 'in-edges))
233            (nodes      ((g 'nodes))))
234       ;; generate differential equations for each state in the transitions system
235      (let ((eqs    (fold (lambda (s ax) 
236                            (if (= (first cnode) (first s) ) ax
237                                (let* ((out   (out-edges (first s)))
238                                       (in    (in-edges (first s)))
239                                       (open? (eq? (second s) open))
240                                       (name  (nest-name (lookup-def (second s) node-subs))))
241                                  (let* ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
242                                                       `(+ (neg ,(sum (map third out)))
243                                                           ,(sum (map third in))))
244                                                      ((and (not (null? out)) (null? in))
245                                                       `(neg ,(sum (map third out))))
246                                                      ((and (null? out) (not (null? in)))
247                                                       (sum (map third in)))))
248                                         (fbody0 (rhsexpr/C++ rhs1)))
249                                   
250                                    (cons (list name (canonicalize-expr/C++ fbody0)) ax)
251                                    ))))
252                          (list) nodes)))
253        eqs))))
254           
255       
256
257(define (poset->asgn-eq-defs poset sys)
258  (fold-right
259   (lambda (lst ax)
260     (fold  (lambda (x ax) 
261              (match-let (((i . n)  x))
262                         (let ((en (hash-table-ref sys n)))
263                           (if (nemo:quantity? en)
264                               (cases nemo:quantity en
265                                      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
266                                      (else  ax))
267                               ax))))
268            ax lst))
269   (list) poset))
270
271(define (poset->rate-eq-defs poset sys)
272  (fold-right
273   (lambda (lst ax)
274     (fold  (lambda (x ax) 
275              (match-let (((i . n)  x))
276                         (let ((en (hash-table-ref sys n)))
277                           (if (nemo:quantity? en)
278                               (cases nemo:quantity en
279
280                                      (REACTION  (name initial open transitions conserve power) 
281                                                 (append (reaction-transition-eqs name initial open transitions 
282                                                                                  conserve power) ax))
283                                     
284                                      (RATE (name initial rhs power)
285                                            (let ((fbody0  (rhsexpr/C++ rhs))
286                                                  (dy      (nest-name name) ))
287                                              (cons (list dy (canonicalize-expr/C++ fbody0)) ax)))
288
289                                      (else  ax))
290                               ax))))
291            ax lst))
292   (list) poset))
293
294
295(define (poset->reaction-eq-defs poset sys)
296  (fold-right
297   (lambda (lst ax)
298     (fold  (lambda (x ax) 
299              (match-let (((i . n)  x))
300                         (let ((en (hash-table-ref sys n)))
301                           (if (nemo:quantity? en)
302                               (cases nemo:quantity en
303                                      (REACTION  (name initial open transitions conserve power) 
304                                                 (cons (reaction-eq name open transitions conserve) ax))
305                                      (else  ax))
306                               ax))))
307            ax lst))
308   (list) poset))
309
310
311(define (poset->init-defs poset sys)
312  (fold-right
313   (lambda (lst ax)
314     (fold-right
315      (lambda (x ax) 
316              (match-let (((i . n)  x))
317                         (let ((en (hash-table-ref sys n)))
318                           (if (nemo:quantity? en)
319                               (cases nemo:quantity en
320                                      (REACTION  (name initial open transitions conserve power)
321                                                 (if (nemo:rhs? initial)
322                                                     (cons* (state-init name initial) 
323                                                            (state-init (nest-state-name name open) name) ax) 
324                                                     ax))
325
326                                      (RATE  (name initial rhs power)
327                                             (if (and initial (nemo:rhs? initial))
328                                                 (cons (state-init name initial) ax) 
329                                                 ax))
330
331                                      (else  ax))
332                               ax))))
333            ax lst))
334   (list) poset))
335
336
337(define (poset->state-conserve-eq-defs poset sys)
338  (fold-right
339   (lambda (lst ax)
340     (fold  (lambda (x ax) 
341              (match-let (((i . n)  x))
342                         (let ((en (hash-table-ref sys n)))
343                           (if (nemo:quantity? en)
344                               (cases nemo:quantity en
345                                      (REACTION (name initial open transitions conserve power)
346                                                (if (and (list? conserve) (every nemo:conseq? conserve))
347                                                    (cons (state-conseqs (nest-name name) transitions conserve
348                                                                        nest-state-name) ax) 
349                                                    ax))
350                                      (else  ax))
351                               ax))))
352            ax lst))
353   (list) poset))
354
355
356(define (find-locals defs)
357  (concatenate
358   (map (lambda (def) 
359          (match def 
360                 (('let bnds body) 
361                  (let ((bexprs (map second bnds)))
362                    (concatenate (list (map first bnds) 
363                                       (find-locals bexprs )
364                                       (find-locals (list body))))))
365                 (('if c t e)      (append (find-locals (list t)) (find-locals (list e))))
366                 ((s . rest)       (find-locals rest))
367                 (else (list)))) 
368        defs)))
369
370
371(define (rate/reaction-power sys n)
372  (let ((en (hash-table-ref sys n)))
373    (if (nemo:quantity? en)
374        (cases nemo:quantity en
375               (REACTION  (name initial open transitions conserve power) 
376                          power)
377               (RATE    (name initial rhs power)
378                        power)
379               (else  #f)) 
380        #f)))
381
382
383(define (bucket-partition p lst)
384  (let loop ((lst lst) (ax (list)))
385    (if (null? lst) ax
386        (let ((x (car lst)))
387          (let bkt-loop ((old-bkts ax) (new-bkts (list)))
388            (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
389                (if (p x (caar old-bkts ))
390                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
391                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
392
393
394(define (make-define-fn sysname )
395  (lambda (indent n proc)
396
397    (let ((lst (procedure-data proc))
398          (indent+ (+ 2 indent)))
399
400      (let ((rt       (or (lookup-def 'rt lst) 'double))
401            (formals  (lookup-def 'formals lst))
402            (vars     (lookup-def 'vars lst))
403            (consts   (lookup-def 'consts lst))
404            (body     (lookup-def 'body lst))
405            (rv       (gensym 'rv)))
406
407        (let ((argument-list 
408               (append
409                (if (null? vars) '() (map (lambda (x) (s+ "double " (nest-name x))) vars))
410                '("void* pnode"))))
411          (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) "{" ))
412          (let* ((body0 (rhsexpr/C++ body))
413                 (body1 (canonicalize-expr/C++ body0))
414                 (lbs   (enum-bnds body1 (list))))
415            (pp indent+ (double ,rv ";"))
416            (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs) ";")))
417            (if (not (null? consts)) 
418                (begin (pp indent+ (double ,(slp ", " (map nest-name consts)) ";")
419                           ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));"))
420                       (for-each (lambda (x) (pp indent+ (,(nest-name x) = ,(s+ "node.P_." (nest-name x) ";")))) consts)
421                       ))
422            (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
423            (pp indent+ ,(s+ "return " rv ";"))
424            ))
425        (pp indent "}"))
426    )))
427
428
429(define (output-dy sysname const-defs state-index-map 
430                   rate-eq-defs reaction-eq-defs asgn-eq-defs
431                   pool-ions mcap i-eqs v-eq indent indent+)
432
433
434  (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
435                      (,(slp ", " `("double t"
436                                    "const double y[]"
437                                    "double f[]"
438                                    "void* pnode"
439                                    )))
440                      #\{
441                      ))
442
443
444  (let* (
445         (asgn-eq-defs 
446          (map
447           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
448           asgn-eq-defs))
449
450         (rate-eq-defs 
451          (map
452           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
453           rate-eq-defs))
454
455         (reaction-eq-defs 
456          (map
457           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
458           reaction-eq-defs))
459
460         (i-eqs 
461          (map
462           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
463           i-eqs))
464
465         (v-eq
466          (and v-eq
467               (list (first v-eq) (canonicalize-expr/C++ (second v-eq)))))
468
469         (eqs 
470          (append
471           
472           asgn-eq-defs
473           reaction-eq-defs
474                 
475           (map (lambda (pool-ion)
476                  (let ((n (third pool-ion))
477                        (b (first pool-ion)))
478                    (list n b)))
479                pool-ions)))
480         
481         (eq-dag 
482          (map (lambda (def)
483                 (cons (first def) (enum-freevars (second def) '() '())))
484               eqs))
485
486         (eq-order
487          (reverse
488           (topological-sort eq-dag 
489                             (lambda (x y) (string=? (->string x) (->string y))))))
490
491         (eq-locals  (find-locals 
492                      (map second
493                           (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
494                               (append i-eqs rate-eq-defs eqs)))))
495         )
496
497
498
499    (pp indent+ (double ,(slp ", " (delete-duplicates 
500                                    (map ->string
501                                         (append
502                                          eq-locals
503                                          eq-order
504                                          (map first state-index-map)
505                                          (map first i-eqs)
506                                          (map first const-defs)))
507                                    string=?))
508                        ";"))
509
510    (pp indent+ ,nl
511        ("// S is shorthand for the type that describes the model state ")
512        (typedef ,(s+ sysname "::State_") "S;")
513       
514        ,nl
515       
516        ("// cast the node ptr to an object of the proper type")
517        ("assert(pnode);")
518        ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
519       
520        ,nl
521       
522        ("// y[] must be the state vector supplied by the integrator, ")
523        ("// not the state vector in the node, node.S_.y[]. ")
524        ,nl
525        )
526   
527    (for-each (lambda (def)
528                (let ((n (first def)) )
529                  (pp indent+
530                      ,(expr->string/C++ 
531                        (s+ "node.P_." n) n))))
532              const-defs)
533
534    (let ((vi (lookup-def 'v state-index-map)))
535      (if vi (pp indent+ ,(expr->string/C++ (s+ "y[" vi "]") 'v))))
536
537    (for-each (lambda (def)
538                (let ((n (first def)) )
539                  (pp indent+
540                      ,(expr->string/C++ 
541                        (s+ "y[" (lookup-def n state-index-map) "]") n))))
542              rate-eq-defs)
543
544    (for-each (lambda (n)
545                (let ((b  (lookup-def n eqs)))
546                  (if b (pp indent+ 
547                            ,(expr->string/C++ b (nest-name n))))))
548              eq-order)
549
550    (for-each (lambda (def)
551                (let ((n (s+ "f[" (lookup-def (first def) state-index-map) "]")) 
552                      (b (second def)))
553                  (pp indent+ ,(s+ "// rate equation " (first def))
554                      ,(expr->string/C++ b n))))
555              rate-eq-defs)
556   
557    (for-each (lambda (def) 
558                (pp indent+ ,(expr->string/C++ (second def) (first def)))) 
559            i-eqs)
560   
561    (if v-eq
562        (let ((vi (lookup-def 'v state-index-map)))
563          (if vi
564              (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "f[" vi "]")))
565              )))
566   
567    (pp indent+ ,nl ("return GSL_SUCCESS;"))
568    (pp indent  #\})
569
570    ))
571
572
573
574(define (output-parameters sysname const-defs indent indent+)
575
576  (let* ((parameter-defs
577          (map
578           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
579           const-defs))
580         (parameter-locals  (find-locals (map second parameter-defs))))
581
582    (pp indent ,nl (,(s+ sysname "::Parameters_::Parameters_") () ":"))
583   
584    (let ((const-exprs
585           (slp ",\n"
586                (map (lambda (def)
587                       (let* ((n  (first def)) (b (second def)))
588                         (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
589                     const-defs) )))
590     
591      (if (not (null? parameter-locals)) 
592          (pp indent+ (double ,(slp ", " parameter-locals) ";")))
593
594      (pp indent+ ,const-exprs)
595     
596      (pp indent ("{}"))
597     
598      )))
599
600
601(define (output-init sysname state-index-map steady-state-index-map 
602                     const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
603                     reaction-eq-defs i-eqs v-eq pool-ions perm-ions indent indent+)
604
605 
606  (let* ((N (length state-index-map))
607
608         (asgn-eq-defs
609          (map
610           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
611           asgn-eq-defs))
612         
613         (init-eq-defs
614          (map
615           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
616           init-eq-defs))
617
618         (reaction-eq-defs 
619          (map
620           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
621           reaction-eq-defs))
622
623         (i-eqs 
624          (map
625           (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
626           i-eqs))
627
628         (v-eq
629          (and v-eq
630               (list (first v-eq) (canonicalize-expr/C++ (second v-eq)))))
631         
632         (init-eqs 
633          (append
634           
635           asgn-eq-defs
636           init-eq-defs
637
638           (map (lambda (pool-ion)
639                  (let ((n (third pool-ion))
640                        (b (first pool-ion)))
641                    (list n b)))
642                pool-ions)))
643
644         (init-dag 
645          (map (lambda (def)
646                 (cons (first def) (enum-freevars (second def) '() '())))
647               init-eqs))
648
649         (init-order
650          (reverse
651           (topological-sort init-dag 
652                             (lambda (x y) (string=? (->string x) (->string y))))))
653
654         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
655         )
656   
657    (pp indent ,nl (,(s+ sysname "::State_::State_") (const Parameters_& p) ": r_(0)"))
658   
659    (pp indent  #\{)
660
661    (pp indent+ (double ,(slp ", " (delete-duplicates
662                                    (map ->string 
663                                         (append
664                                          init-locals
665                                          init-order
666                                          (map first state-index-map) 
667                                          (map first i-eqs)
668                                          (map first const-defs)
669                                          (map first reaction-eq-defs)
670                                          ))
671                                    string=?))
672                              ";"))
673
674    (for-each (lambda (def)
675                (let ((n (first def)) )
676                  (pp indent+
677                      ,(expr->string/C++ 
678                        (s+ "p." n) n))))
679              const-defs)
680
681    (let ((vi (lookup-def 'v state-index-map))
682          (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
683      (if (and vi vrest) 
684          (pp indent+ ,(expr->string/C++  vrest 'v ))))
685
686    (for-each (lambda (n)
687                (let ((b  (lookup-def n init-eqs)))
688                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
689              init-order)
690   
691    (for-each (lambda (def)
692                (let* ((n  (first def)) 
693                       (ni (lookup-def n state-index-map)))
694                  (if ni (pp indent+ ,(expr->string/C++ n  (s+ "y_[" ni "]"))))))
695              init-eq-defs)
696#|
697    (if (not (null? steady-state-index-map))
698        (begin
699          (gsl-fminimizer sysname N indent+)
700          (for-each
701           (lambda (def)
702             (let* ((n (first def))
703                    (si (lookup-def n steady-state-index-map))
704                    (ni (lookup-def n state-index-map)))
705               (if si (begin
706                        (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n))
707                        (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")")
708                                                       (s+ "y_[" ni "]")))))))
709           rate-eq-defs)
710      (pp indent+ "gsl_vector_free (ys);")
711      ))
712|#
713
714    (for-each
715     (lambda (def)
716       (let ((n (first def)) (b (second def)))
717         (if (not (lookup-def n init-eq-defs))
718             (pp indent+ ,(expr->string/C++ b n)))))
719     reaction-eq-defs)
720   
721    (for-each
722     (lambda (def) 
723       (pp indent+ ,(expr->string/C++ (second def) (first def))))
724     i-eqs)
725
726    (if v-eq
727        (let ((vi (lookup-def 'v state-index-map)))
728          (if vi
729              (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "y_[" vi "]"))))))
730     
731    (pp indent  #\})
732   
733    (pp indent ,nl (,(s+ sysname "::State_::State_") (const State_& s) ": r_(s.r_)"))
734    (pp indent  #\{)
735    (pp indent+ (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N)))
736    (pp indent  #\})
737
738    (pp indent ,nl (,(s+ sysname "::State_& " sysname "::State_::operator=(const State_& s)")))
739
740    (pp indent  #\{)
741
742
743    (pp indent+ (,(sprintf #<<EOF
744   assert(this != &s)
745     for ( size_t i = 0 ; i < ~A ; ++i )
746       y_[i] = s.y_[i];
747     r_ = s.r_;
748     return *this;
749EOF
750N)))
751
752     (pp indent  #\})
753
754))
755
756
757
758(define (gsl-fminimizer sysname N indent)
759
760  (pp indent (,(sprintf #<<EOF
761
762  int minimizer_status, minimizer_iterc;
763
764  const gsl_multimin_fdfminimizer_type *T;
765  gsl_multimin_fdfminimizer *S;
766
767  gsl_vector *ys;
768  gsl_multimin_function_fdf F;
769
770  F.f   = &~A_dynamics;
771  F.df  = &~A_dynamics;
772  F.fdf = &~A_dynamics;
773  F.n   = ~A;
774  F.params = (void *)&p;
775
776  ys = gsl_vector_alloc (~A);
777
778  for (int i = 0; i < ~A; i++)
779  {
780     gsl_vector_set (ys, i, 0.0);
781  }
782
783  T = gsl_multimin_fdfminimizer_conjugate_fr;
784  S = gsl_multimin_fdfminimizer_alloc (T, ~A);
785
786  gsl_multimin_fdfminimizer_set (S, &F, ys, 0.01, 1e-4);
787
788  do
789    {
790      minimizer_iterc++;
791      minimizer_status = gsl_multimin_fdfminimizer_iterate (S);
792
793      if (minimizer_status)
794        break;
795
796      minimizer_status = gsl_multimin_test_gradient (S->gradient, 1e-3);
797
798    }
799  while (minimizer_status == GSL_CONTINUE && minimizer_iterc < 100);
800
801  gsl_multimin_fdfminimizer_free (S);
802
803EOF
804sysname sysname sysname N N N N)))
805)
806
807
808
809(define (output-accessors+modifiers
810         sysname state-index-map const-defs asgn-eq-defs rate-eq-defs 
811         reaction-eq-defs i-eqs v-eq pool-ions perm-ions indent indent+)
812
813  (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
814  (pp indent  #\{)
815
816  (for-each
817   (lambda (def)
818     (let ((n (first def)))
819       (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
820   const-defs)
821
822  (pp indent  #\})
823
824  (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
825  (pp indent  #\{)
826
827  (for-each
828   (lambda (def)
829     (let ((n (first def)))
830       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))))
831   const-defs)
832
833  (pp indent  #\})
834
835  (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
836  (pp indent  #\{)
837
838  (for-each
839   (lambda (def)
840     (let ((n (first def)) (i (second def)))
841       (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
842   state-index-map)
843
844  (pp indent  #\})
845
846  (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
847  (pp indent  #\{)
848
849  (for-each
850   (lambda (def)
851     (let ((n (first def)) (i (second def)))
852       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
853   state-index-map)
854
855
856  (pp indent  #\})
857
858)
859
860(define (output-buffers+nodes
861         sysname state-index-map steady-state-index-map 
862         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
863         reaction-eq-defs i-eqs pool-ions perm-ions 
864         synapse-info
865         indent indent+)
866
867  (let ((N (length state-index-map)))
868
869  (pp indent ,nl  (,(sprintf #<<EOF
870~A::Buffers_::Buffers_(~A& n)
871    : logger_(n),
872      s_(0),
873      c_(0),
874      e_(0)
875{
876    // Initialization of the remaining members is deferred to
877    // init_buffers_().
878}
879EOF
880  sysname sysname)))
881
882  (pp indent  ,nl  (,(sprintf #<<EOF
883~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
884    : logger_(n),
885      s_(0),
886      c_(0),
887      e_(0)
888{
889    // Initialization of the remaining members is deferred to
890    // init_buffers_().
891}
892EOF
893  sysname sysname)))
894
895
896  (pp indent  ,nl  (,(sprintf #<<EOF
897~A::~A()
898    : Archiving_Node(), 
899      P_(), 
900      S_(P_),
901      B_(*this)
902{
903    recordablesMap_.create();
904}
905EOF
906  sysname sysname)))
907
908
909  (pp indent  ,nl  (,(sprintf #<<EOF
910~A::~A(const ~A& n)
911    : Archiving_Node(n), 
912      P_(n.P_), 
913      S_(n.S_),
914      B_(n.B_, *this)
915{
916}
917EOF
918  sysname sysname sysname)))
919
920  (pp indent (,(s+ sysname "::~" sysname "()")))
921  (pp indent #<<EOF
922{
923    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
924    if ( B_.s_ ) gsl_odeiv_step_free(B_.s_);
925    if ( B_.c_ ) gsl_odeiv_control_free(B_.c_);
926    if ( B_.e_ ) gsl_odeiv_evolve_free(B_.e_);
927}
928EOF
929  )
930
931  (pp indent ,nl (,(sprintf #<<EOF
932  void ~A::init_node_(const Node& proto)
933{
934    const ~A& pr = downcast<~A>(proto);
935    P_ = pr.P_;
936    S_ = pr.S_;
937}
938EOF
939  sysname sysname sysname)))
940
941  (pp indent  ,nl (,(sprintf #<<EOF
942void ~A::init_state_(const Node& proto)
943{
944    const ~A& pr = downcast<~A>(proto);
945    S_ = pr.S_;
946}
947EOF
948  sysname sysname sysname)))
949
950  (pp indent ,nl #<<EOF
951void ~A::init_buffers_()
952{
953EOF
954)
955     (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
956      (for-each
957       (lambda (psc)
958         (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
959       pscs))
960
961    (pp indent+ (,(sprintf #<<EOF
962    B_.currents_.clear();           
963    Archiving_Node::clear_history();
964
965    B_.logger_.reset();
966
967    B_.step_ = Time::get_resolution().get_ms();
968    B_.IntegrationStep_ = B_.step_;
969
970    B_.I_stim_ = 0.0;
971
972    static const gsl_odeiv_step_type* T1 = gsl_odeiv_step_rkf45;
973 
974    if ( B_.s_ == 0 )
975      B_.s_ = gsl_odeiv_step_alloc (T1, ~A);
976    else
977      gsl_odeiv_step_reset(B_.s_);
978   
979    if ( B_.c_ == 0 ) 
980      B_.c_ = gsl_odeiv_control_y_new (1e-3, 0.0);
981    else
982      gsl_odeiv_control_init(B_.c_, 1e-3, 0.0, 1.0, 0.0);
983   
984    if ( B_.e_ == 0 ) 
985      B_.e_ = gsl_odeiv_evolve_alloc(~A);
986    else
987      gsl_odeiv_evolve_reset(B_.e_);
988 
989    B_.sys_.function  = ~A_dynamics;
990    B_.sys_.jacobian  = 0;
991    B_.sys_.dimension = ~A;
992    B_.sys_.params    = reinterpret_cast<void*>(this);
993}
994EOF
995  N N sysname N)))
996
997  (let ((iv (lookup-def 'v state-index-map)))
998    (if iv
999        (pp indent ,nl (,(sprintf #<<EOF
1000void ~A::calibrate()
1001{
1002    B_.logger_.init()
1003    V_.RefractoryCounts_ = 20;
1004    V_.U_old_ = S_.y_[~A];
1005}
1006EOF
1007  sysname iv)))
1008  ))
1009
1010))
1011
1012
1013(define (output-update
1014         sysname state-index-map steady-state-index-map 
1015         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1016         reaction-eq-defs i-eqs pool-ions perm-ions 
1017         synapse-info
1018         indent indent+)
1019 
1020  (let* ((vi (lookup-def 'v state-index-map))
1021         (threshold (sprintf #<<EOF
1022            // (threshold && maximum)
1023            if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
1024              {
1025                S_.r_ = V_.RefractoryCounts_;
1026               
1027                set_spiketime(Time::step(origin.get_steps()+lag+1));
1028               
1029                SpikeEvent se;
1030                network()->send(*this, se, lag);
1031              }
1032EOF
1033        vi vi)))
1034   
1035    (pp indent  (,(sprintf #<<EOF
1036void ~A::update(Time const & origin, const long_t from, const long_t to)
1037EOF
1038  sysname)))
1039
1040    (pp indent  #\{)
1041
1042    (pp indent+ (,(sprintf #<<EOF
1043assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
1044    assert(from < to);
1045
1046    for ( long_t lag = from ; lag < to ; ++lag )
1047      {
1048   
1049        double tt = 0.0 ; //it's all relative!
1050        V_.U_old_ = S_.y_[~A];
1051
1052   
1053        // adaptive step integration
1054        while (tt < B_.step_)
1055        {
1056          const int status = gsl_odeiv_evolve_apply(B_.e_, B_.c_, B_.s_, 
1057                                 &B_.sys_,              // system of ODE
1058                                 &tt,                   // from t...
1059                                  B_.step_,             // ...to t=t+h
1060                                 &B_.IntegrationStep_ , // integration window (written on!)
1061                                  S_.y_);               // neuron state
1062
1063          if ( status != GSL_SUCCESS )
1064            throw GSLSolverFailure(get_name(), status);
1065        }
1066EOF
1067vi)))
1068
1069    (let ((isyns (lookup-def 'i-synapses synapse-info))
1070          (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
1071
1072      (for-each (lambda (isyn psc)
1073                  (let ((gi (lookup-def (nest-name (fourth isyn)) state-index-map)))
1074                    (if gi
1075                        (pp indent+ (,(sprintf "      S_.y_[~A] += B_.spike_~A.get_value(lag);"
1076                                               gi (second psc)))
1077                        ))))
1078                isyns pscs))
1079
1080    (pp indent+ (,(sprintf #<<EOF
1081        // sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
1082        // refractory?
1083        if (S_.r_)
1084          {
1085            --S_.r_;
1086          }
1087        else
1088          {
1089           ~A
1090          }
1091   
1092        // set new input current
1093        B_.I_stim_ = B_.currents_.get_value(lag);
1094
1095        // log state data
1096        B_.logger_.record_data(origin.get_steps() + lag);
1097
1098      }
1099EOF
1100  (if (lookup-def 'V_t const-defs) threshold "")
1101  )))
1102 
1103  (pp indent  #\})
1104
1105))
1106
1107
1108(define (output-handle
1109         sysname state-index-map steady-state-index-map 
1110         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1111         reaction-eq-defs i-eqs pool-ions perm-ions 
1112         synapse-info
1113         indent indent+)
1114 
1115  (pp indent  ,nl (,(sprintf #<<EOF
1116void ~A::handle(SpikeEvent & e)
1117  {
1118    int flag;
1119    assert(e.get_delay() > 0);
1120    flag = 0;
1121EOF
1122   sysname)))
1123
1124    (let ((isyns (lookup-def 'i-synapses synapse-info))
1125          (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
1126      (for-each
1127       (lambda (isyn psc)
1128         (pp indent+ 
1129             (,(sprintf #<<EOF
1130    if ((flag==0) && (e.get_weight() > P_.~A))
1131      {
1132        B_.spike_~A.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
1133                              abs(e.get_weight()) * e.get_multiplicity());
1134        flag = 1;
1135      }
1136EOF
1137             (nest-name (second isyn))
1138             (second psc)))))
1139       isyns pscs))
1140
1141  (pp indent ,nl #\} ,nl)
1142
1143  (pp indent  ,nl (,(sprintf #<<EOF
1144void ~A::handle(CurrentEvent& e)
1145  {
1146    assert(e.get_delay() > 0);
1147
1148    const double_t c=e.get_current();
1149    const double_t w=e.get_weight();
1150
1151    B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), 
1152                        w *c);
1153  }
1154
1155void ~A::handle(DataLoggingRequest& e)
1156  {
1157    B_.logger_.handle(e);
1158  }
1159EOF
1160  sysname sysname)))
1161
1162)
1163
1164
1165(define (output-recordables
1166         sysname state-index-map steady-state-index-map 
1167         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1168         reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
1169 
1170  (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
1171  (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
1172  (pp indent #\{)
1173  (for-each
1174   (lambda (def)
1175     (let ((n (first def)) (i (second def)))
1176       (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);" 
1177                              (->string n) sysname sysname (string-upcase (->string n)))))
1178       ))
1179   state-index-map)
1180
1181  (if (lookup-def 'v state-index-map)
1182      (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);" 
1183                             sysname sysname (string-upcase (->string 'v))))))
1184
1185  (pp indent #\})
1186  )
1187
1188
1189(define (output-prelude sysname steady-state-index-map indent)
1190  (pp indent (,(sprintf "#include \"~A.h\"
1191#include \"exceptions.h\"
1192#include \"network.h\"
1193#include \"dict.h\"
1194#include \"integerdatum.h\"
1195#include \"doubledatum.h\"
1196#include \"dictutils.h\"
1197#include \"numerics.h\"
1198#include <limits>
1199
1200#include \"universal_data_logger_impl.h\"
1201
1202#include <iomanip>
1203#include <iostream>
1204#include <cstdio>
1205#include <cmath>
1206
1207#include <gsl/gsl_errno.h>
1208#include <gsl/gsl_matrix.h>
1209#include <gsl/gsl_sf_exp.h>
1210"  sysname)))
1211  (if (not (null? steady-state-index-map))
1212      (pp indent ("#include <gsl/gsl_multimin.h>")))
1213  )
1214
1215
1216(define (output-header 
1217         sysname state-index-map steady-state-index-map 
1218         const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1219         reaction-eq-defs i-eqs pool-ions perm-ions 
1220         synapse-info
1221         indent indent+)
1222
1223  (pp indent ("#include \"nest.h\"
1224#include \"event.h\"
1225#include \"archiving_node.h\"
1226#include \"ring_buffer.h\"
1227#include \"connection.h\"
1228#include \"universal_data_logger.h\"
1229#include \"recordables_map.h\"
1230#include <gsl/gsl_odeiv.h>
1231"))
1232 
1233  (pp indent "namespace nest {" ,nl)
1234
1235  (pp indent (extern "\"C\""
1236                     int ,(s+ sysname "_dynamics")
1237                     "(double, const double*, double*, void*)" #\;) ,nl)
1238
1239  (pp indent 
1240      (,(sprintf "class ~A:" sysname)
1241       "public Archiving_Node { "))
1242  (pp indent+
1243       ("public: ")
1244       (,(s+ "~" sysname "();"))
1245       (,(s+ sysname "(const " sysname "&);"))
1246       (,(s+ sysname "();")))
1247
1248  (pp indent (#<<EOF
1249    using Node::connect_sender;
1250    using Node::handle;
1251
1252    port check_connection(Connection&, port);
1253   
1254    void handle(SpikeEvent &);
1255    void handle(CurrentEvent &);
1256    void handle(DataLoggingRequest &);
1257   
1258    port connect_sender(SpikeEvent &, port);
1259    port connect_sender(CurrentEvent &, port);
1260    port connect_sender(DataLoggingRequest &, port);
1261   
1262    void get_status(DictionaryDatum &) const;
1263    void set_status(const DictionaryDatum &);
1264   
1265    void init_node_(const Node& proto);
1266    void init_state_(const Node& proto);
1267    void init_buffers_();
1268    void calibrate();
1269   
1270    void update(Time const &, const long_t, const long_t);
1271EOF
1272))
1273
1274  (pp indent (,(sprintf #<<EOF
1275    // make dynamics function quasi-member
1276    friend int ~A_dynamics(double, const double*, double*, void*);
1277
1278    // The next two classes need to be friends to access the State_ class/member
1279    friend class RecordablesMap<~A>;
1280    friend class UniversalDataLogger<~A>;
1281EOF
1282sysname sysname sysname)))
1283
1284 
1285  (pp indent+ ,nl "struct Parameters_ { ")
1286
1287  (for-each
1288   (lambda (x) (pp indent+ (,(sprintf "double ~A;" (first x)))))
1289   const-defs)
1290
1291  (pp indent+ 
1292      ("Parameters_();")
1293      ("void get(DictionaryDatum&) const;")
1294      ("void set(const DictionaryDatum&);"))
1295
1296  (pp indent+ "};")
1297
1298  (pp indent+ ,nl "struct State_ { ")
1299
1300  (pp indent+ ,nl 
1301      "enum StateVecElems {"
1302      (,(slp ", " (map (lambda (x) 
1303                         (let ((n (string-upcase (->string (first x))))
1304                               (ni (second x)))
1305                           (s+ n " = " ni)))
1306                       state-index-map)))
1307      "};"
1308      (,(sprintf "double y_[~A];" (length state-index-map)))
1309      "int_t     r_;"
1310      "State_(const Parameters_& p);"
1311      "State_(const State_& s);"
1312     
1313      "State_& operator=(const State_& s);"
1314
1315      "void get(DictionaryDatum&) const;"
1316      "void set(const DictionaryDatum&, const Parameters_&);")
1317  (pp indent+ "};")
1318
1319
1320  (pp indent+ ,nl #<<EOF
1321    struct Variables_ {
1322      int_t    RefractoryCounts_;
1323      double   U_old_; // for spike-detection
1324    };
1325EOF
1326)
1327
1328  (pp indent+ ,nl "struct Buffers_ { ")
1329
1330  (pp indent+ ,nl 
1331      (,(sprintf "Buffers_(~A&);" sysname))
1332      (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
1333      (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))
1334
1335  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
1336    (for-each
1337     (lambda (psc)
1338       (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
1339     pscs))
1340
1341  (pp indent+ ,nl
1342#<<EOF
1343
1344  RingBuffer currents_;
1345
1346  gsl_odeiv_step*    s_;    //!< stepping function
1347  gsl_odeiv_control* c_;    //!< adaptive stepsize control function
1348  gsl_odeiv_evolve*  e_;    //!< evolution function
1349  gsl_odeiv_system   sys_;  //!< struct describing system
1350
1351  double_t step_;           //!< step size in ms
1352  double   IntegrationStep_;//!< current integration time step, updated by GSL
1353
1354  /** 
1355   * Input current injected by CurrentEvent.
1356   * This variable is used to transport the current applied into the
1357   * _dynamics function computing the derivative of the state vector.
1358   * It must be a part of Buffers_, since it is initialized once before
1359   * the first simulation, but not modified before later Simulate calls.
1360   */
1361  double_t I_stim_;
1362EOF
1363)
1364  (pp indent+ "};")
1365
1366  (pp indent+
1367      "template <State_::StateVecElems elem>"
1368      "double_t get_y_elem_() const { return S_.y_[elem]; }"
1369      "Parameters_ P_;"
1370      "State_      S_;"
1371      "Variables_  V_;"
1372      "Buffers_    B_;"
1373      (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
1374      "}; ")
1375
1376  (pp indent+ (,(sprintf #<<EOF
1377  inline
1378  port ~A::check_connection(Connection& c, port receptor_type)
1379  {
1380    SpikeEvent e;
1381    e.set_sender(*this);
1382    c.check_event(e);
1383    return c.get_target()->connect_sender(e, receptor_type);
1384  }
1385
1386  inline
1387  port ~A::connect_sender(SpikeEvent&, port receptor_type)
1388  {
1389    if (receptor_type != 0)
1390      throw UnknownReceptorType(receptor_type, get_name());
1391    return 0;
1392  }
1393 
1394  inline
1395  port ~A::connect_sender(CurrentEvent&, port receptor_type)
1396  {
1397    if (receptor_type != 0)
1398      throw UnknownReceptorType(receptor_type, get_name());
1399    return 0;
1400  }
1401
1402  inline
1403  port ~A::connect_sender(DataLoggingRequest& dlr, 
1404                                      port receptor_type)
1405  {
1406    if (receptor_type != 0)
1407      throw UnknownReceptorType(receptor_type, get_name());
1408    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
1409  }
1410
1411  inline
1412    void ~A::get_status(DictionaryDatum &d) const
1413  {
1414    P_.get(d);
1415    S_.get(d);
1416    Archiving_Node::get_status(d);
1417
1418    (*d)[names::recordables] = recordablesMap_.get_list();
1419
1420    def<double_t>(d, names::t_spike, get_spiketime_ms());
1421  }
1422
1423  inline
1424    void ~A::set_status(const DictionaryDatum &d)
1425  {
1426    Parameters_ ptmp = P_;  // temporary copy in case of errors
1427    ptmp.set(d);                       // throws if BadProperty
1428    State_      stmp = S_;  // temporary copy in case of errors
1429    stmp.set(d, ptmp);                 // throws if BadProperty
1430
1431    // We now know that (ptmp, stmp) are consistent. We do not 
1432    // write them back to (P_, S_) before we are also sure that 
1433    // the properties to be set in the parent class are internally 
1434    // consistent.
1435    Archiving_Node::set_status(d);
1436
1437    // if we get here, temporaries contain consistent set of properties
1438    P_ = ptmp;
1439    S_ = stmp;
1440
1441    calibrate();
1442  }
1443EOF
1444  sysname sysname sysname 
1445  sysname sysname sysname)))
1446
1447  (pp indent "}" ,nl)
1448  )
1449
1450(define (nemo:nest-translator sys . rest)
1451
1452  (define (cid x)  (second x))
1453  (define (cn x)   (first x))
1454
1455  (let-optionals rest ((dirname "."))
1456  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
1457    (let ((imports  ((dis 'imports)  sys))
1458          (exports  ((dis 'exports)  sys)))
1459      (let* ((indent      0)
1460             (indent+     (+ 2 indent ))
1461
1462             (eval-const  (dis 'eval-const))
1463             (sysname     (nest-name ((dis 'sysname) sys)))
1464             (prefix      (->string sysname))
1465             (deps*       ((dis 'depgraph*)  sys))
1466             (consts      ((dis 'consts)     sys))
1467             (asgns       ((dis 'asgns)      sys))
1468             (states      ((dis 'states)     sys))
1469             (reactions   ((dis 'reactions)  sys))
1470             (defuns      ((dis 'defuns)     sys))
1471             (components  ((dis 'components) sys))
1472             
1473             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
1474             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
1475
1476             (const-defs       (filter-map
1477                                (lambda (nv)
1478                                  (and (not (member (first nv) nest-builtin-consts))
1479                                       (let ((v1 (canonicalize-expr/C++ (second nv))))
1480                                         (list (nest-name (first nv)) v1))))
1481                                consts))
1482             
1483             (gate-complex-info    (nemo:gate-complex-query sys))
1484             (synapse-info         (nemo:post-synaptic-conductance-query sys))
1485             (i-synapses           (lookup-def 'i-synapses synapse-info))
1486
1487             (gate-complexes  (lookup-def 'gate-complexes gate-complex-info))
1488             (perm-ions       (map (match-lambda ((comp i e erev) `(,comp ,(nest-name i) ,(nest-name e) ,erev)))
1489                                   (lookup-def 'perm-ions gate-complex-info)))
1490             (acc-ions        (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out)))))
1491                                   (lookup-def 'acc-ions gate-complex-info)))
1492             (epools          (lookup-def 'pool-ions gate-complex-info))
1493             (pool-ions       (map (lambda (lst) (map nest-name lst)) epools))
1494
1495             (i-gates         (lookup-def 'i-gates gate-complex-info))
1496
1497             (capcomp         (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
1498             (mcap            (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
1499               
1500             (i-eqs (filter-map
1501                     (lambda (gate-complex) 
1502                       
1503                       (let* ((label             (first gate-complex))
1504                              (n                 (second gate-complex))
1505                              (subcomps          ((dis 'component-subcomps) sys n))
1506                              (acc               (lookup-def 'accumulating-substance subcomps))
1507                              (perm              (lookup-def 'permeating-ion subcomps))
1508                              (permqs            (and perm ((dis 'component-exports) sys (cid perm))))
1509                              (pore              (lookup-def 'pore subcomps))
1510                              (permeability      (lookup-def 'permeability subcomps))
1511                              (gate              (lookup-def 'gate subcomps))
1512                              (sts               (and gate ((dis 'component-exports) sys (cid gate)))))
1513                         
1514                         (if (not (or pore permeability))
1515                             (nemo:error 'nemo:nest-translator ": ion channel definition " label
1516                                         "lacks any pore or permeability components"))
1517
1518
1519                         (cond ((and perm permeability gate)
1520                                     (let* ((i     (nest-name (s+ 'i (cn perm))))
1521                                            (pmax  (car ((dis 'component-exports) sys (cid permeability))))
1522                                           
1523                                            (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
1524                                            (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
1525
1526                                            (gion  `(* ,pmax ,@sptms)))
1527                                         (list i #f gion (nest-name (s+ 'i_ label) ))))
1528
1529                               ((and perm pore gate)
1530                                (case (cn perm)
1531                                  ((non-specific)
1532                                   (let* ((i     (nest-name 'i))
1533                                          (e     (car permqs))
1534                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
1535
1536                                          (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
1537                                          (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
1538
1539                                          (gion  `(* ,gmax ,@sptms)))
1540                                     (list i e gion  (nest-name (s+ 'i_ label) ))))
1541
1542                                  (else
1543                                   (let* ((i     (nest-name (s+ 'i (cn perm))))
1544                                          (e     (car permqs))
1545                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
1546
1547                                          (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
1548                                          (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
1549
1550                                          (gion  `(* ,gmax ,@sptms)))
1551                                     (list i e gion (nest-name (s+ 'i_ label) ))))))
1552                               
1553                               ((and perm pore)
1554                                (case (cn perm)
1555                                  ((non-specific)
1556                                   (let* ((i     (nest-name 'i))
1557                                          (e     (car permqs))
1558                                          (gmax  (car ((dis 'component-exports) sys (cid pore)))))
1559                                     (list i e gmax (nest-name (s+ 'i_ label) ))))
1560                                  (else
1561                                   (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))))
1562                               
1563                               ((and acc pore gate)
1564                                (let* ((i     (nest-name (s+ 'i (cn acc))))
1565                                       (gmax  (car ((dis 'component-exports) sys (cid pore))))
1566
1567                                       (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
1568                                       (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
1569                                       
1570                                       (gion  `(* ,gmax ,@sptms)))
1571                                  (list i #f gion  (nest-name (s+ 'i_ label) ))))
1572                               (else (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))
1573                               )))
1574                     gate-complexes))
1575
1576             (i-names (delete-duplicates (map first i-eqs)))
1577               
1578             (i-eqs  (fold  (lambda (i-gate ax) 
1579                              (let ((i-gate-var (first i-gate)))
1580                                (cons (list (nest-name 'i) #f i-gate-var (s+ 'i_ (second i-gate))) ax)))
1581                            i-eqs i-gates))
1582
1583             (i-eqs  (fold  (lambda (i-synapse ax) 
1584                              (cons (list (first i-synapse) (third i-synapse) (fourth i-synapse) (gensym 'i_syn )) ax))
1585                            i-eqs 
1586                            i-synapses))
1587
1588             (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
1589             
1590             (i-eqs  (fold (lambda (b ax) 
1591                             (match b 
1592                                    ((and ps ((i e gion ii) . rst)) 
1593                                     (let loop ((ps ps) (summands (list)) (eqs (list)))
1594                                       (if (null? ps)
1595                                           
1596                                           (let* ((sum0  (sum summands))
1597                                                  (sum1  (rhsexpr/C++ sum0))
1598                                                  (sum2  (canonicalize-expr/C++ sum1)))
1599                                             (append eqs (list (list i sum2)) ax))
1600                                           
1601                                           (match-let (((i e gion ii) (car ps)))
1602                                                      (loop (cdr ps) 
1603                                                            (cons ii summands) 
1604                                                            (let* ((expr0 (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
1605                                                                   (expr1 (canonicalize-expr/C++ expr0)))
1606                                                              (cons (list ii expr1) eqs)))))))
1607                                   
1608                                    ((i e gion ii)
1609                                     (let* ((expr0  (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
1610                                            (expr1  (canonicalize-expr/C++ expr0)))
1611                                       (cons (list i expr1) ax)))
1612                                   
1613                                    (else ax)))
1614                           (list) i-bkts))
1615
1616             (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
1617             
1618             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys)))
1619             
1620             (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
1621             
1622             (init-eq-defs       (poset->init-defs poset sys))
1623             
1624             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
1625                                      (poset->state-conserve-eq-defs poset sys)))
1626             
1627             (globals          (map nest-name 
1628                                    (delete-duplicates (append
1629                                                        exports
1630                                                        (map second  perm-ions)
1631                                                        (map third   perm-ions)
1632                                                        (map second  acc-ions)
1633                                                        (map third   acc-ions)
1634                                                        (map fourth  acc-ions)
1635                                                        (map second  pool-ions)
1636                                                        (map third   pool-ions)
1637                                                        (map first   imports) 
1638                                                        (map first   const-defs)
1639                                                        (map first   asgn-eq-defs)
1640                                                        (map (lambda (gate-complex) (nest-name (s+ 'i_ (first gate-complex)))) gate-complexes )
1641                                                        (map (lambda (i-gate) (nest-name (s+ 'i_ (second i-gate)))) i-gates )
1642                                                        ))))
1643
1644             (parameters (map nest-name (map first const-defs)))
1645
1646             (v-eq    (if (and mcap (member 'v globals))
1647                          (list 'v (rhsexpr/C++ `(/ (neg ,(sum (cons "(-node.B_.I_stim_)" i-names))) ,mcap)))
1648                          (list 'v 0.0)))
1649             
1650             (v-init-eq  (if (and mcap (member 'v globals))
1651                             (list 'v (rhsexpr/C++ `(/ (neg ,(sum i-names)) ,mcap)))
1652                             (list 'v 0.0)))
1653             
1654             (state-index-map  (let ((acc (fold (lambda (def ax)
1655                                                  (let ((st-name (first def)))
1656                                                    (list (+ 1 (first ax)) 
1657                                                          (cons `(,st-name ,(first ax)) (second ax)))))
1658                                                (list 0 (list)) 
1659                                                (if (member 'v globals)
1660                                                    (cons (list 'v) rate-eq-defs)
1661                                                    rate-eq-defs)
1662                                                )))
1663                                 (second acc)))
1664             
1665             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
1666                                                         (let ((st-name (first def)))
1667                                                           (if (not (alist-ref st-name init-eq-defs))
1668                                                               (list (+ 1 (first ax)) 
1669                                                                     (cons `(,st-name ,(first ax)) (second ax)))
1670                                                               ax)))
1671                                                       (list 0 (list)) 
1672                                                       rate-eq-defs)))
1673                                        (second acc)))
1674             
1675             (dfenv 
1676              (map (lambda (x) (let ((n (first x)))
1677                                 (list n (nest-name (s+ "d_" n )))))
1678                   defuns))
1679
1680             )
1681       
1682       
1683       
1684        (for-each
1685         (lambda (a)
1686           (let ((acc-ion   (car a)))
1687             (if (assoc acc-ion perm-ions)
1688                 (nemo:error 'nemo:nest-translator 
1689                             ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
1690         acc-ions)
1691
1692           (for-each
1693            (lambda (p)
1694              (let ((pool-ion  (car p)))
1695                (if (assoc pool-ion perm-ions)
1696                    (nemo:error 'nemo:nest-translator 
1697                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
1698            pool-ions)
1699
1700           (let ((cpp-output  (open-output-file (make-output-fname dirname prefix ".cpp")))
1701                 (hpp-output  (open-output-file (make-output-fname dirname prefix ".h"))))
1702
1703             ;; include files and other prelude definitions
1704             (with-output-to-port cpp-output
1705               (lambda ()
1706                 (output-prelude sysname steady-state-index-map indent)
1707                 ))
1708
1709             (with-output-to-port cpp-output
1710               (lambda () (pp indent "namespace nest {" ,nl)))
1711
1712             ;; user-defined functions
1713             (let ((define-fn  (make-define-fn sysname)))
1714               (for-each (lambda (fndef) 
1715                           (if (not (member (car fndef) builtin-fns))
1716                               (with-output-to-port cpp-output
1717                                 (lambda ()
1718                                   (apply define-fn (cons indent fndef))
1719                                   (pp indent ,nl)))
1720                               ))
1721                         defuns))
1722
1723
1724             ;; derivative function
1725             (with-output-to-port cpp-output
1726               (lambda ()
1727                 (output-dy sysname const-defs state-index-map 
1728                            rate-eq-defs reaction-eq-defs asgn-eq-defs
1729                            pool-ions mcap i-eqs v-eq
1730                            indent indent+)
1731                 ))
1732
1733             ;;  system recordables
1734             (with-output-to-port cpp-output
1735               (lambda ()
1736                 (output-recordables
1737                  sysname state-index-map steady-state-index-map 
1738                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1739                  reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
1740                 (pp indent ,nl)
1741                 ))
1742
1743             ;; system parameters
1744             (with-output-to-port cpp-output
1745               (lambda ()
1746                 (output-parameters sysname const-defs
1747                                    indent indent+)
1748                 ))
1749
1750             ;; initial values function
1751             (with-output-to-port cpp-output
1752               (lambda ()
1753                 (output-init sysname state-index-map steady-state-index-map 
1754                              const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1755                              reaction-eq-defs i-eqs v-init-eq
1756                              pool-ions perm-ions indent indent+)
1757                 (pp indent ,nl)
1758                 ))
1759
1760             ;; accessors and modifiers for states and parameters
1761             (with-output-to-port cpp-output
1762               (lambda ()
1763                 (output-accessors+modifiers
1764                  sysname state-index-map 
1765                  const-defs asgn-eq-defs rate-eq-defs 
1766                  reaction-eq-defs i-eqs v-eq pool-ions perm-ions 
1767                  indent indent+)
1768                 (pp indent ,nl)
1769                 ))
1770
1771             
1772             ;; buffer and node initialization
1773             (with-output-to-port cpp-output
1774               (lambda ()
1775                 (output-buffers+nodes
1776                  sysname state-index-map steady-state-index-map 
1777                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1778                  reaction-eq-defs i-eqs pool-ions perm-ions 
1779                  synapse-info
1780                  indent indent+)
1781                 (pp indent ,nl)
1782                 ))
1783
1784             ;; state update
1785             (with-output-to-port cpp-output
1786               (lambda ()
1787                 (output-update
1788                  sysname state-index-map steady-state-index-map 
1789                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1790                  reaction-eq-defs i-eqs pool-ions perm-ions 
1791                  synapse-info indent indent+)
1792                 (pp indent ,nl)
1793                 ))
1794
1795             ;; spike handler
1796             (with-output-to-port cpp-output
1797               (lambda ()
1798                 (output-handle
1799                  sysname state-index-map steady-state-index-map 
1800                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1801                  reaction-eq-defs i-eqs pool-ions perm-ions 
1802                  synapse-info
1803                  indent indent+)
1804                 (pp indent ,nl)
1805                 ))
1806
1807             (with-output-to-port cpp-output
1808               (lambda () (pp indent "}" ,nl)))
1809
1810
1811             (with-output-to-port hpp-output
1812               (lambda ()
1813                 (output-header
1814                  sysname state-index-map steady-state-index-map 
1815                  const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
1816                  reaction-eq-defs i-eqs pool-ions perm-ions 
1817                  synapse-info
1818                  indent indent+)
1819                 (pp indent ,nl)
1820                 ))
1821             
1822             (close-output-port cpp-output)
1823             (close-output-port hpp-output)
1824           
1825           ))
1826      ))
1827  ))
1828
1829
1830)
Note: See TracBrowser for help on using the repository browser.