;;       
;; TODO: finish gsl-fminimizer
;; 
;; An extension for translating NEMO models to NEST code.
;;
;; Copyright 2011-2013 Ivan Raikov and the Okinawa Institute of Science and Technology
;;
;; This program is free software: you can redistribute it and/or
;; modify it under the terms of the GNU General Public License as
;; published by the Free Software Foundation, either version 3 of the
;; License, or (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful, but
;; WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;; General Public License for more details.
;;
;; A full copy of the GPL license can be found at
;; <http://www.gnu.org/licenses/>.
;;

(module nemo-nest

	(nemo:nest-translator)

	(import scheme chicken utils data-structures ports extras srfi-1 srfi-13 srfi-69)
	
	(require-extension lolevel posix matchable strictly-pretty 
			   varsubst datatype nemo-core nemo-utils
			   nemo-gate-complex nemo-geometry nemo-synapse nemo-defaults)

(define nest-builtin-consts  `(params))


(define C++-ops
  `(+ - * / > < <= >= =))


(define builtin-fns
  `(+ - * / pow neg abs atan asin acos sin cos exp ln
      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
      > < <= >= = and or round ceiling floor max min
      ))


(define (nest-name s)
  (let ((cs (string->list (->string s))))
    (let loop ((lst (list)) (cs cs))
      (if (null? cs) (string->symbol (list->string (reverse lst)))
	  (let* ((c (car cs))
		 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
			   (else #\_))))
	    (loop (cons c1 lst) (cdr cs)))))))

	    
(define (rhsexpr/C++ expr)
  (match expr 
	 (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/C++ x)) es)))
	 (('pow x y)  (if (and (integer? y)  (positive? y))
			  (if (> y 1)  (let ((tmp (gensym "x")))
					 `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
			      x)
			    (if (and (number? y) (zero? y)) 1.0 expr)))
	 ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (nest-name s))
					     (map (lambda (x) (rhsexpr/C++ x)) es)) expr))
	 (id          (if (symbol? id) (nest-name id) id))))


(define (nest-state-name n s)
  (nest-name (s+ n s)))


(define-syntax pp
  (syntax-rules ()
    ((pp indent val ...) (ppf indent (quasiquote val) ...))))


(define group/C++   (doc:block 2 (doc:text "(") (doc:text ")")))
(define block/C++   (doc:block 2 (doc:text "{") (doc:text "}")))
(define (stmt/C++ x) 
  (match x
	 (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
	 (else  (doc:cons x (doc:text ";")))))


(define (ifthen/C++ c e1 e2)
  (doc:nest 2
    (doc:connect (doc:group (doc:connect (doc:text "if") c))
		 (doc:connect (doc:nest 2 e1)
			      (doc:nest 2 (doc:connect 
					   (doc:text "else") 
					   e2))))
    ))


(define (letblk/C++ e1 e2)
  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
	((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
	(else (doc:connect (doc:group (doc:nest 2 (stmt/C++ e1)))
			   (doc:group (doc:nest 2 e2))))))
	

(define (format-op/C++ indent op args)
  (let ((op1 (doc:text (->string op))))
    (if (null? args) op1
	(match args
	       ((x)      (doc:concat (list op1 x)))
	       ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
	       ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
	       (lst      (let* ((n   (length lst))
				(n/2 (inexact->exact (round (/ n 2)))))
			   (doc:concat 
			    (intersperse 
			     (list (format-op/C++ indent op (take lst n/2 )) op1 
				   (format-op/C++ indent op (drop lst n/2 )))
			     (doc:space)))))))))

(define (format-fncall/C++ indent op args)
  (let ((op1 (doc:text (->string op))))
    (doc:cons op1 (group/C++ ((doc:list indent identity (lambda () (doc:text ", "))) args)))))

(define (name-normalize expr)
  (match expr 
	 (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
	 (('let bs e)
	  `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
	 ((f . es) 
	  (cons (if (member f builtin-fns) f (nest-name f)) (map name-normalize es)))
	 ((? symbol? ) (nest-name expr))
	 ((? atom? ) expr)))


(define (canonicalize-expr/C++ expr)
  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term)))
    (let* ((expr1 (if-convert expr))
	   (expr2 (subst-convert expr1 subst-empty))
	   (expr3 (let-lift expr2))
	   (expr4 (name-normalize expr3)))
      expr4)))


(define (format-expr/C++ indent expr . rest)  
  (let-optionals rest ((rv #f))
   (let ((indent+ (+ 2 indent)))
    (match expr
       (('let bindings body)
	(letblk/C++
	 (fold-right 
	   (lambda (x ax)
	     (letblk/C++
	      (match (second x)
		     (('if c t e)
		      (ifthen/C++
		       (group/C++ (format-expr/C++ indent c))
		       (block/C++ (format-expr/C++ indent t (first x)))
		       (block/C++ (format-expr/C++ indent e (first x)))))
		     (else 
		      (stmt/C++
		       (format-op/C++ indent+ " = "
					 (list (format-expr/C++ indent (first x) )
					       (format-expr/C++ indent (second x)))))))
	      ax))
	   (doc:empty) bindings)
	 (match body
		(('let _ _) (format-expr/C++ indent body rv))
		(else
		 (let ((body1 (doc:nest indent (format-expr/C++ indent body))))
		   (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) body1)))
		       body1))))))
       
       (('if . rest) (error 'format-expr/C++ "invalid if statement " expr))

       ((op . rest)  
       (let ((op (case op ((ln) 'log) ((abs) 'fabs) (else op))))
	 (let ((fe
		(if (member op C++-ops)
		    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
			  (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
			  (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
		      (case op
			((/)  
			 (format-op/C++ indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/C++ indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if (or mul? plmin?) (group/C++ fx) fx)))) rest)))
			((*)  
			 (format-op/C++ indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/C++ indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if plmin? (group/C++ fx) fx)))) rest)))
			
			(else
			 (format-op/C++ indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/C++ indent+ x))) fx)) rest)))))
		    
		    (let ((op (case op ((neg) '-) (else op))))
		      (format-fncall/C++ indent op (map (lambda (x) (format-expr/C++ indent+ x)) rest))))))
	   (if rv 
	       (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
	       fe))))
      
      (else  (let ((fe (doc:text (->string expr))))
	       (if rv 
		   (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
		   fe)))))))
	       
	  
(define (expr->string/C++ x . rest)
  (let-optionals rest ((rv #f) (width 72))
    (sdoc->string (doc:format width (format-expr/C++ 2 x rv)))))


(define (make-define-fn sysname )
  (lambda (indent n proc)

    (let ((lst (procedure-data proc))
	  (indent+ (+ 2 indent)))

      (let ((rt       (or (lookup-def 'rt lst) 'double))
	    (formals  (lookup-def 'formals lst))
	    (vars     (lookup-def 'vars lst))
	    (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
	    (body     (lookup-def 'body lst))
	    (rv       (gensym 'rv)))

	(let ((argument-list 
	       (append
		(if (null? vars) '() (map (lambda (x) (s+ "double " (nest-name x))) vars))
		'("const void* params"))))
	  (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) "{" ))
	  (let* ((body0 (rhsexpr/C++ body))
		 (body1 (canonicalize-expr/C++ (add-params-to-fncall body0 builtin-fns)))
		 (lbs   (enum-bnds body1 (list))))
	    (pp indent+ (double ,rv ";"))
	    (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs) ";")))
	    (if (not (null? consts)) 
		(begin (pp indent+ (double ,(slp ", " (map (compose nest-name car) consts)) ";")
			   ("const " ,(s+ sysname "::Parameters_") "& p =  " 
                            " *(reinterpret_cast< const " ,(s+ sysname "::Parameters_") "*>(params));"))
		       (for-each (lambda (x) (let ((n (car x))) 
                                               (pp indent+ (,(nest-name n) = ,(s+ "p." (nest-name n) ";"))))) 
                                 consts)
		       ))
	    (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
	    (pp indent+ ,(s+ "return " rv ";"))
	    ))
	(pp indent "}"))
    )))


(define (make-define-fn-header sysname )
  (lambda (indent n proc)

    (let ((lst (procedure-data proc))
	  (indent+ (+ 2 indent)))

      (let ((rt       (or (lookup-def 'rt lst) 'double))
	    (formals  (lookup-def 'formals lst))
	    (vars     (lookup-def 'vars lst))
	    (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
	    (body     (lookup-def 'body lst))
	    (rv       (gensym 'rv)))

	(let ((argument-list 
	       (append
		(if (null? vars) '() (map (lambda (x) (s+ "double " (nest-name x))) vars))
		'("const void* params"))))
	  (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) ";" ))
	  ))
      )))



(define (ith v i) (sprintf "Ith(~A,~A)" v i))


(define (output-dy sysname method imports const-defs state-index-map 
		   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
		   pool-ions i-eqs v-eq indent indent+)


  (case method
    ((leapfrog)
     (pp indent ,nl (
		     ,(s+ "extern \"C\" int " sysname "_dynamics")
		     (,(slp ", " `("double t"
				   "const double y[]"
				   "double f[]"
				   "void* pnode"
				   )))
		     #\{
		     )))
    ((cvode)
     (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
		     (,(slp ", " `("double t"
				   "N_Vector y"
				   "N_Vector f"
				   "void* pnode"
				   )))
		     #\{
		     )))
    ((gsl)
     (pp indent ,nl (,(s+ "extern \"C\" int " sysname "_dynamics")
		     (,(slp ", " `("double t"
				   "const double y[]"
				   "double f[]"
				   "void* pnode"
				   )))
		     #\{
		     )))
     )

  (let* (
	 (i-eqs 
	  (map 
	   (lambda (def) (list (first def) 
			       (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
	   i-eqs))

	 (v-eq
	  (and v-eq
	       (list (first v-eq) 
		     (add-params-to-fncall (canonicalize-expr/C++ (second v-eq)) builtin-fns))))

	 (eqs 
	  (append 

	   external-eq-defs
	   asgn-eq-defs
	   reaction-eq-defs
		 
	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-in pool-ion))
			(b (pool-ion-inq pool-ion)))
		    (list n b)))
		pool-ions)

	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-out pool-ion))
			(b (pool-ion-outq pool-ion)))
		    (list n b)))
		pool-ions)

           i-eqs
	  ))
	 
	 (eq-dag 
	  (map (lambda (def)
		 (cons (first def) (enum-freevars (second def) '() '())))
	       eqs))

	 (eq-order
	  (reverse
	   (topological-sort eq-dag 
			     (lambda (x y) (string=? (->string x) (->string y))))))

	 (eq-locals  (find-locals 
		      (map second
			   (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
			       (append i-eqs rate-eq-defs eqs)))))
	 )


    (pp indent+ (double ,(slp ", " (delete-duplicates 
				    (map (compose ->string nest-name)
					 (filter (lambda (x) 
						   (not (member x nest-builtin-consts)))
						 (append 
						  eq-locals
						  eq-order
						  (map first i-eqs)
						  (map first external-eq-defs)
						  (map first state-index-map)
						  (map first const-defs)
						  )))
				    string=?))
			";"))

    (pp indent+ ,nl
	("// S is shorthand for the type that describes the model state ")
	(typedef ,(s+ sysname "::State_") "S;")
	
	,nl
	
	("// cast the node ptr to an object of the proper type")
	("assert(pnode);")
	("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
	(,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
	
	,nl

	("// params is a reference to the model parameters ")
	(const struct ,(s+ sysname "::Parameters_") "*params;")
	("params = &(node.P_);")
	
	,nl

	("// state is a reference to the model state ")
	(struct ,(s+ sysname "::State_") "*state;")
	("state = &(vnode.S_);")
	
	,nl

	("// y[] must be the state vector supplied by the integrator, ")
	("// not the state vector in the node, node.S_.y[]. ")
	,nl
	)

    (for-each (lambda (def)
		(let ((n (first def)) )
		  (pp indent+
		      ,(expr->string/C++ 
			(s+ "params->" n) n))))
	      const-defs)

    (for-each (lambda (def) 
		(pp indent+ ,(expr->string/C++ (second def) (first def)))) 
	    i-eqs)

    (let ((vi (lookup-def 'v state-index-map)))
      (if vi (pp indent+ ,(expr->string/C++ (ith 'y vi) 'v))))

    (for-each (lambda (def)
		(let* ((n (first def))
		       (ni (lookup-def n state-index-map)))
		  (pp indent+ ,(expr->string/C++ (ith 'y ni) (nest-name n)))))
	      rate-eq-defs)

    (for-each (lambda (n)
		(let ((b  (lookup-def n eqs)))
		  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
	      eq-order)

    (for-each (lambda (def)
		(let* ((n (first def))
		       (b (second def))
		       (fv (ith 'f (lookup-def n state-index-map)))
		       )
		  (pp indent+ ,(s+ "// rate equation " n)
		      ,(expr->string/C++ b fv))))
	      rate-eq-defs)
    
    
    (if v-eq
	(let ((vi (lookup-def 'v state-index-map)))
	  (if vi
	      (pp indent+ ,(expr->string/C++ (second v-eq) (ith 'f vi)))
	      )))

    

    (case method
      ((leapfrog)
	(pp indent+ ,nl ("return 0;")))
      ((cvode)
	(pp indent+ ,nl ("return 0;")))
      ((gsl)
	(pp indent+ ,nl ("return GSL_SUCCESS;")))
      )
    (pp indent  #\})

    ))


(define (output-jac sysname method state-index-map pool-ions i-eqs v-eq indent indent+)

;; Diagonal Jacobian approximation: (f(s+.01) - f(s))/.001 

;; CVODE:
;;  typedef int (*CVDlsDenseJacFn)(long int N, realtype t,
;;                               N_Vector y, N_Vector fy, 
;;                               DlsMat Jac, void *user_data,
;;                               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);

;; GSL:
;; int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);



  (case method

    ((gsl)
     (pp indent ,nl
         (,(s+ "extern \"C\" int " sysname "_jacobian")
          (,(slp ", " `("double t"
                        "const double y[]"
                        "double *dfdy"
                        "double dfdt[]"
                        "void* pnode"
                        ))))
          #\{

          ("// cast the node ptr to an object of the proper type")
          ("assert(pnode);")
          ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
          (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")

          ("// state is a reference to the model state ")
          (struct ,(s+ sysname "::Buffers_") "*b;")
          ("b = &(vnode.B_);")
          (,(sprintf
#<<EOF
   for (int i = 0; i < b->N; i++) 
   {
       b->u[i] = y[i] + 0.01;
   }
   ~A_dynamics(t, b->u, b->jac, pnode);
   for (int i = 0; i < b->N; i++) 
   {
       dfdt[i*b->N + i] = (b->jac[i] - dfdy[i]) / .001;
   }
   return 0;
EOF
   sysname))
		     #\}
		     ))
     )
     
  )



(define (output-parameters sysname imports const-defs external-eq-defs defaults indent indent+)

  (let* ((parameter-defs
	  (map 
	   (lambda (def) (list (first def) 
			       (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
	   const-defs))
	 (parameter-locals  (find-locals (map second parameter-defs))))

    (pp indent ,nl (,(s+ sysname "::Parameters_::Parameters_") () ":"))

    (let ((const-exprs
	    (map (lambda (def)
		   (let* ((n  (first def)) (b (second def)))
		     (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
		 const-defs) )
	     
	  (default-exprs 
	    (map (lambda (def)
		   (let ((n (first def))
			 (b (second def)))
		     (expr->string/C++ (nest-name b) (nest-name n))))
		 defaults))
	   )
      
      (if (not (null? parameter-locals)) 
	  (pp indent+ (double ,(slp ", " parameter-locals) ";")))

      (pp indent+ ,(slp ",\n" const-exprs))
      
      (pp indent  #\{ ,(slp "\n" default-exprs) #\})
      
      ))
  )


(define (output-init sysname state-index-map steady-state-index-map 
		     imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
		     reaction-eq-defs i-eqs pool-ions perm-ions defaults indent indent+)

  
  (let* ((N (length state-index-map))

	 (i-eqs 
	  (map 
	   (lambda (def) (list (first def) 
			       (add-params-to-fncall 
				(canonicalize-expr/C++ (second def)) builtin-fns)))
	   i-eqs))

	 (init-eqs 
	  (append 
	     
	   (map (lambda (def)
		  (let ((n (first def))
			(b (second def)))
		    (list (nest-name n) (nest-name b))))
		external-eq-defs)
	   
	   asgn-eq-defs
	   init-eq-defs

	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-in pool-ion))
			(b (pool-ion-inq pool-ion)))
		    (list n b)))
		pool-ions)

	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-out pool-ion))
			(b (pool-ion-outq pool-ion)))
		    (list n b)))
		pool-ions)
	   ))

	 (init-dag 
	  (map (lambda (def)
		 (cons (first def) (enum-freevars (second def) '() '())))
	       init-eqs))

	 (init-order
	  (reverse
	   (topological-sort init-dag 
			     (lambda (x y) (string=? (->string x) (->string y))))))

	 (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))

	 )
    
    (pp indent ,nl (,(s+ sysname "::State_::State_") (const Parameters_& p) ": r_(0)"))
    
    (pp indent  #\{)

    (pp indent+
	(double ,(slp ", " (delete-duplicates
			    (map ->string 
				 (filter (lambda (x) (not (member x nest-builtin-consts)))
					 (append 
					  init-locals
					  init-order
					  (map first external-eq-defs)
					  (map pool-ion-in pool-ions)
					  (map pool-ion-out pool-ions)
					  (map first i-eqs)
					  (map first state-index-map) 
					  (map first const-defs)
					  (map first reaction-eq-defs)
					  )))
				 string=?))
		      ";")
	("const Parameters_ *params = &p;")
	)

    (pp indent+ ,(sprintf "memset(y_,0,~A*sizeof(double));" N))

    (for-each
     (lambda (n) (pp indent+ ,(expr->string/C++ (s+ "p." n) n)))
     (map (compose nest-name first) const-defs) )

    (let ((vi (lookup-def 'v state-index-map))
	  (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
      (if (and vi vrest) 
	  (pp indent+ ,(expr->string/C++  vrest 'v ))))

    (for-each (lambda (n)
		(let ((b  (lookup-def n init-eqs)))
		  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
	      init-order)
    
    (for-each (lambda (def)
		(let* ((n  (first def)) 
		       (ni (lookup-def n state-index-map)))
		  (if ni (pp indent+ ,(expr->string/C++ n  (s+ "y_[" ni "]"))))))
	      init-eq-defs)

    
#|
    (if (not (null? steady-state-index-map))
	(begin
	  (gsl-fminimizer sysname N indent+)
	  (for-each
	   (lambda (def)
	     (let* ((n (first def)) 
		    (si (lookup-def n steady-state-index-map))
		    (ni (lookup-def n state-index-map)))
	       (if si (begin 
			(pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n))
			(pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")") 
						       (s+ "y_[" ni "]")))))))
	   rate-eq-defs)
      (pp indent+ "gsl_vector_free (ys);")
      ))
|#

    (for-each
     (lambda (def)
       (let ((n (first def)) (b (second def)))
	 (if (not (lookup-def n init-eq-defs))
	     (pp indent+ ,(expr->string/C++ b n)))))
     reaction-eq-defs)
    
    (for-each
     (lambda (def) 
       (pp indent+ ,(expr->string/C++ (second def) (first def))))
     i-eqs)

    (let ((vi (lookup-def 'v state-index-map)))
      (if vi (pp indent+ ,(expr->string/C++ 'v (s+ "y_[" vi "]")))))
      
    (pp indent  #\})

    (pp indent ,nl
	(,(s+ sysname "::State_::State_") (const State_& s) ": r_(0)" )
	)

    (pp indent  #\{)
    (pp indent+ (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N)))

    (pp indent  #\})

    (pp indent ,nl (,(s+ sysname "::State_& " sysname "::State_::operator=(const State_& s)")))

    (pp indent  #\{)


    (pp indent+ (,(sprintf #<<EOF
   assert(this != &s);  
     for ( size_t i = 0 ; i < ~A ; ++i )
       y_[i] = s.y_[i];
EOF
N)))
     (pp indent+ "return *this;")

     (pp indent  #\})
))



(define (gsl-fminimizer sysname N indent)

  (pp indent (,(sprintf #<<EOF

  int minimizer_status, minimizer_iterc;

  const gsl_multimin_fdfminimizer_type *T;
  gsl_multimin_fdfminimizer *S;

  gsl_vector *ys;
  gsl_multimin_function_fdf F;

  F.f   = &~A_dynamics;
  F.df  = &~A_dynamics;
  F.fdf = &~A_dynamics;
  F.n   = ~A;
  F.params = (void *)&p;

  ys = gsl_vector_alloc (~A);

  for (int i = 0; i < ~A; i++)
  {
     gsl_vector_set (ys, i, 0.0);
  }

  T = gsl_multimin_fdfminimizer_conjugate_fr;
  S = gsl_multimin_fdfminimizer_alloc (T, ~A);

  gsl_multimin_fdfminimizer_set (S, &F, ys, 0.01, 1e-4);

  do
    {
      minimizer_iterc++;
      minimizer_status = gsl_multimin_fdfminimizer_iterate (S);

      if (minimizer_status)
        break;

      minimizer_status = gsl_multimin_test_gradient (S->gradient, 1e-3);

    }
  while (minimizer_status == GSL_CONTINUE && minimizer_iterc < 100);

  gsl_multimin_fdfminimizer_free (S);

EOF
sysname sysname sysname N N N N)))
)



(define (output-accessors+modifiers
	 sysname imports state-index-map const-defs asgn-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
         indent indent+)

  (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
  (pp indent  #\{)

  (for-each
   (lambda (def)
     (let ((n (first def)))
       (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
   const-defs)

  (pp indent  #\})

  (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
  (pp indent  #\{)

  (for-each
   (lambda (def)
     (let ((n (first def)))
       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))))
   const-defs)

  (pp indent  #\})

  (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
  (pp indent  #\{)

  (for-each
   (lambda (def)
     (let ((n (first def)) (i (second def)))
       (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
   state-index-map)

  (let ((vi (lookup-def 'v state-index-map)))
    (if vi
        (pp indent+ (,(sprintf "def<double_t>(d, names::V_m, y_[~A]);"  vi) ))
        ))

  (pp indent  #\})

  (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
  (pp indent  #\{)

  (for-each
   (lambda (def)
     (let ((n (first def)) (i (second def)))
       (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))))
   state-index-map)

  (let ((vi (lookup-def 'v state-index-map)))
    (if vi
        (pp indent+ (,(sprintf "updateValue<double_t>(d, names::V_m, y_[~A]);"  vi) ))
        ))

  (pp indent  #\})

)



(define (output-handle
	 sysname state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info
	 indent indent+)
  
  (pp indent  ,nl (,(sprintf #<<EOF
void ~A::handle(SpikeEvent & e)
  {
    int flag;
    assert(e.get_delay() > 0);
    flag = 0;
EOF
   sysname)))

   (let ((isyns (lookup-def 'i-synapses synapse-info))
	  (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
      (for-each
       (lambda (isyn psc)
	 (pp indent+ 
	     (,(sprintf #<<EOF
    if ((flag==0) && (e.get_weight() > P_.~A))
      {
	B_.spike_~A.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
			      abs(e.get_weight()) * e.get_multiplicity());
        flag = 1;
      }
EOF
             (nest-name (second isyn))
	     (second psc)))))
       isyns pscs))

  (pp indent ,nl #\} ,nl)

  (pp indent  ,nl (,(sprintf #<<EOF
void ~A::handle(CurrentEvent& e)
  {
    assert(e.get_delay() > 0);

    const double_t c=e.get_current();
    const double_t w=e.get_weight();

    B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), 
			w *c);
  }

void ~A::handle(DataLoggingRequest& e)
  {
    B_.logger_.handle(e);
  }
EOF
  sysname sysname)))

)


(define (output-recordables
	 sysname state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
  
  (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
  (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
  (pp indent #\{)
  (for-each 
   (lambda (def)
     (let ((n (first def)) (i (second def)))
       (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);" 
			      (->string n) sysname sysname (string-upcase (->string n)))))
       ))
   state-index-map)

  (if (lookup-def 'v state-index-map)
      (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);" 
			     sysname sysname (string-upcase (->string 'v))))))

  (pp indent #\})
  )


(define (output-prelude sysname steady-state-index-map indent)
  (pp indent (,(sprintf "#include \"~A.h\"
#include \"exceptions.h\"
#include \"network.h\"
#include \"dict.h\"
#include \"integerdatum.h\"
#include \"doubledatum.h\"
#include \"dictutils.h\"
#include \"numerics.h\"
#include <limits>

#include \"universal_data_logger_impl.h\"

#include <iomanip>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
"  sysname)))
  )

(define (output-header 
	 sysname method state-index-map steady-state-index-map 
	 defaults external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info
	 indent indent+)

  (pp indent ("#include \"nest.h\"
#include \"event.h\"
#include \"archiving_node.h\"
#include \"ring_buffer.h\"
#include \"connection.h\"
#include \"universal_data_logger.h\"
#include \"recordables_map.h\"
"))

  (case method
    ((leapfrog) 
     (pp indent
	 ("#define Ith(v,i)    (v[i])   /* Ith component in a vector */")
	 ))
    ((cvode) 
     (pp indent
	 ("#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */")
	 ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
	 ("#include <cvode/cvode_diag.h>        /* prototype for CVDiag */")
	 ("#include <sundials/sundials_types.h> /* definition of type realtype */")
	 ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
	 ))
    ((gsl) 
     (pp indent
	 ("#include <gsl/gsl_errno.h>")
	 ("#include <gsl/gsl_matrix.h>")
	 ("#include <gsl/gsl_sf_exp.h>")
	 ("#include <gsl/gsl_odeiv2.h>")
	 ("#define Ith(v,i)    (v[i])")
	 ))
    )

  (pp indent ("#define PI M_PI"))

  (pp indent ,nl)

  (pp indent "namespace nest {" ,nl)

  
  (case method
    ((cvode) 
     (pp indent #<<EOF

/* 
 * Check function return value.
 *    opt == 0 means SUNDIALS function allocates memory so check if
 *             returned NULL pointer
 *    opt == 1 means SUNDIALS function returns a flag so check if
 *             flag >= 0
 *    opt == 2 means function allocates memory so check if returned
 *             NULL pointer 
 */

static int check_flag(void *flagvalue, const char *funcname, int opt)
{
  int *errflag;

  /* Check if CVode function returned NULL pointer - no memory allocated */
  if (opt == 0 && flagvalue == NULL) {
    fprintf(stderr, 
            "\nCVode error: %s() failed - returned NULL pointer\n\n",
            funcname);
    return(1); }

  /* Check if flag < 0 */
  else if (opt == 1) {
    errflag = (int *) flagvalue;
    if (*errflag < 0) {
      fprintf(stderr, 
              "\nCVode error: %s() failed with flag = %d\n\n",
              funcname, *errflag);
      return(1); }}

  /* Check if function returned NULL pointer - no memory allocated */
  else if (opt == 2 && flagvalue == NULL) {
    fprintf(stderr, 
            "\nMemory error: %s() failed - returned NULL pointer\n\n",
            funcname);
    return(1); }

  return(0);
}
EOF
)))

  
  (case method
    ((leapfrog) 
     (pp indent 
	 ("typedef int (*rhsfn_t)(double t, const double *y, double *ydot, void *user_data);") ,nl
	 (extern "\"C\""
		     int ,(s+ sysname "_dynamics")\
		     "(double, const double *, double *, void*)" #\;) ,nl))
    ((cvode) 
     (pp indent (extern "\"C\""
		     int ,(s+ sysname "_dynamics")
		     "(double, const N_Vector, N_Vector, void*)" #\;) ,nl))
    (else
     (pp indent (extern "\"C\""
		     int ,(s+ sysname "_dynamics")
		     "(double, const double*, double*, void*)" #\;) ,nl))
    )

  (pp indent 
      (,(sprintf "class ~A:" sysname)
       "public Archiving_Node { "))

  (pp indent+
       ("public: ")
       (,(s+ "~" sysname "();"))
       (,(s+ sysname "(const " sysname "&);"))
       (,(s+ sysname "();")))

  (pp indent (#<<EOF
    using Node::connect_sender;
    using Node::handle;

    port check_connection(Connection&, port);
    
    void handle(SpikeEvent &);
    void handle(CurrentEvent &);
    void handle(DataLoggingRequest &); 
    
    port connect_sender(SpikeEvent &, port);
    port connect_sender(CurrentEvent &, port);
    port connect_sender(DataLoggingRequest &, port);
    
    void get_status(DictionaryDatum &) const;
    void set_status(const DictionaryDatum &);
    
    void init_node_(const Node& proto);
    void init_state_(const Node& proto);
    void init_buffers_();
    void calibrate();
    
    void update(Time const &, const long_t, const long_t);
EOF
))

  
  (case method
    ((cvode) 
     (pp indent (friend int ,(s+ sysname "_dynamics")
		     "(double, const N_Vector, N_Vector, void*)" #\;) ,nl))
    (else
     (pp indent (friend int ,(s+ sysname "_dynamics")
		     "(double, const double*, double*, void*)" #\;) ,nl))
    )

  (pp indent (,(sprintf #<<EOF
    // The next two classes need to be friends to access the State_ class/member
    friend class RecordablesMap<~A>;
    friend class UniversalDataLogger<~A>;
EOF
sysname sysname)))

  
  (pp indent+ ,nl "struct Parameters_ { ")

  (for-each 
   (lambda (x) (pp indent+ (,(sprintf "double ~A;" x))))
   (append 
    (map (compose nest-name first) const-defs)
    (map (compose nest-name first) defaults)))

  (pp indent+ 
      ("Parameters_();")
      ("void get(DictionaryDatum&) const;")
      ("void set(const DictionaryDatum&);"))

  (pp indent+ "};")

  (pp indent+ ,nl "struct State_ { ")

  (pp indent+ ,nl 
      "enum StateVecElems {"
      (,(slp ", " (map (lambda (x) 
			 (let ((n (string-upcase (->string (first x))))
			       (ni (second x)))
			   (s+ n " = " ni)))
		       state-index-map)))
      "};"
      "int_t     r_;"

      (,(sprintf "double y_[~A];" (length state-index-map)))

      "State_(const Parameters_& p);" 
      "State_(const State_& s);"
      "State_& operator=(const State_& s);"
      "void get(DictionaryDatum&) const;"
      "void set(const DictionaryDatum&, const Parameters_&);")

  (pp indent+ "};")


  (pp indent+ ,nl #<<EOF
struct Variables_ {
      int_t    RefractoryCounts_;
      double   U_old_; // for spike-detection
};
EOF
)

  (pp indent+ ,nl "struct Buffers_ { ")

  (pp indent+ ,nl 
      (,(sprintf "Buffers_(~A&);" sysname))
      (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
      (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))

  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
    (for-each
     (lambda (psc)
       (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
     pscs))

  (pp indent+ ,nl
  (,(case method

    ((leapfrog)
#<<EOF
  rhsfn_t sys_; // pointer to function that computes dv/dt
  double *u, *v, *dvdt;  // intermediate state vectors used by leapfrog method
  unsigned int  N;  // size of state vector used by leapfrog method
EOF
)

    ((cvode)
#<<EOF
  N_Vector y;    //!< current state vector used by CVode
  void *   sys_;  //!< CVode control structure
EOF
)

    ((gsl)
#<<EOF
  gsl_odeiv2_step*    s_;    //!< stepping function
  gsl_odeiv2_control* c_;    //!< adaptive stepsize control function
  gsl_odeiv2_evolve*  e_;    //!< evolution function
  gsl_odeiv2_system   sys_;  //!< struct describing system
  unsigned int N;  // size of state vector used by Jacobian
  double *u, *jac;  // intermediate state vectors used for Jacobian approximation
EOF
)

   )))


  (pp indent+ ,nl
#<<EOF
  RingBuffer currents_;

  double_t step_;           //!< step size in ms
  double   IntegrationStep_;//!< current integration time step, updated by solver

  /** 
   * Input current injected by CurrentEvent.
   * This variable is used to transport the current applied into the
   * _dynamics function computing the derivative of the state vector.
   * It must be a part of Buffers_, since it is initialized once before
   * the first simulation, but not modified before later Simulate calls.
   */
  double_t I_stim_;
EOF
)
  (pp indent+ "};")

  (pp indent+
      "template <State_::StateVecElems elem>"
      "double_t get_y_elem_() const { return S_.y_[elem]; }"
      "Parameters_ P_;"
      "State_      S_;"
      "Variables_  V_;"
      "Buffers_    B_;"
      (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
      "}; ")

  (pp indent+ (,(sprintf #<<EOF
  inline
  port ~A::check_connection(Connection& c, port receptor_type)
  {
    SpikeEvent e;
    e.set_sender(*this);
    c.check_event(e);
    return c.get_target()->connect_sender(e, receptor_type);
  }

  inline
  port ~A::connect_sender(SpikeEvent&, port receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return 0;
  }
 
  inline
  port ~A::connect_sender(CurrentEvent&, port receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return 0;
  }

  inline
  port ~A::connect_sender(DataLoggingRequest& dlr, 
				      port receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
  }

  inline
    void ~A::get_status(DictionaryDatum &d) const
  {
    P_.get(d);
    S_.get(d);
    Archiving_Node::get_status(d);

    (*d)[names::recordables] = recordablesMap_.get_list();

    def<double_t>(d, names::t_spike, get_spiketime_ms());
  }

  inline
    void ~A::set_status(const DictionaryDatum &d)
  {
    Parameters_ ptmp = P_;  // temporary copy in case of errors
    ptmp.set(d);                       // throws if BadProperty
    State_      stmp = S_;  // temporary copy in case of errors
    stmp.set(d, ptmp);                 // throws if BadProperty

    // We now know that (ptmp, stmp) are consistent. We do not 
    // write them back to (P_, S_) before we are also sure that 
    // the properties to be set in the parent class are internally 
    // consistent.
    Archiving_Node::set_status(d);

    // if we get here, temporaries contain consistent set of properties
    P_ = ptmp;
    S_ = stmp;

    calibrate();
  }
EOF
  sysname sysname sysname 
  sysname sysname sysname)))

  (pp indent "}" ,nl)
  )

(define (output-update
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs defaults i-eqs pool-ions perm-ions 
	 synapse-info 
	 indent indent+)
  
  (let* ((vi (lookup-def 'v state-index-map))
	 (threshold (sprintf #<<EOF
	    // (threshold && maximum)
	    if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
	      {
		S_.r_ = V_.RefractoryCounts_;
		
		set_spiketime(Time::step(origin.get_steps()+lag+1));
		
		SpikeEvent se;
		network()->send(*this, se, lag);
	      }
EOF
        vi vi)))
    
    (pp indent  (,(sprintf #<<EOF
void ~A::update(Time const & origin, const long_t from, const long_t to)
EOF
  sysname)))

    (pp indent  #\{)

    (pp indent+ (,(sprintf #<<EOF
assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
    assert(from < to);

    double tout;
    long_t current_steps = origin.get_steps();

    for ( long_t lag = from ; lag < to ; ++lag )
      {
        double h = B_.step_;    
	double tt = 0.0 ; 
	V_.U_old_ = S_.y_[~A];
EOF
vi)))

    (case method
      ((leapfrog)
       (pp indent+ (,(sprintf #<<EOF

        tt = (current_steps+lag)*h;
        tout = (current_steps+lag+1)*h;

	while (tt < tout)
	{
	   B_.sys_ (tt, S_.y_, B_.dvdt, reinterpret_cast<void*>(this));

	   for (int i = 0; i < B_.N; i++) 
           {
              B_.u[i] = S_.y_[i] + 0.5 * h * B_.dvdt[i];
           }

	   B_.sys_ (tt+0.5*h, B_.u, B_.dvdt, reinterpret_cast<void*>(this));

           for (int i = 0; i < B_.N; i++) 
           {  
              B_.v[i] = S_.y_[i] + h * B_.dvdt[i];
           }

	   tt += h;
	}

        for (int i = 0; i < B_.N; i++)
        {
           S_.y_[i] = B_.v[i];
        }

EOF
 ))
))
      ((cvode)
       (pp indent+ (,(sprintf #<<EOF

        int N = NV_LENGTH_S (B_.y);
        tout = (current_steps+lag+1)*h;

	// adaptive step integration
	while (tt < tout)
	{
	  const int status = CVode (B_.sys_, tout, B_.y, &tt, CV_NORMAL);

	  if ( status != CV_SUCCESS )
	    throw GSLSolverFailure(get_name(), status);
	}

        for (int i = 0; i < N; i++)
        {
           S_.y_[i] = Ith(B_.y,i);
        }
EOF
 ))
))
      ((gsl)
       (pp indent+ (,(sprintf #<<EOF
			      
	// adaptive step integration
	while (tt < h)
	{
	  const int status = gsl_odeiv2_evolve_apply(B_.e_, B_.c_, B_.s_, 
				 &B_.sys_,              // system of ODE
				 &tt,                   // from t...
				  h,             // ...to t=t+h
				 &B_.IntegrationStep_ , // integration window (written on!)
				  S_.y_);	        // neuron state

	  if ( status != GSL_SUCCESS )
	    throw GSLSolverFailure(get_name(), status);
	}
EOF
))))
      )

    (let ((isyns (lookup-def 'i-synapses synapse-info))
	  (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))

      (for-each (lambda (isyn psc)
		  (let ((gi (lookup-def (nest-name (fourth isyn)) state-index-map)))
		    (if gi
			(pp indent+ (,(sprintf "      S_.y_[~A] += B_.spike_~A.get_value(lag);"
					       gi (second psc)))
			))))
		isyns pscs))

    (pp indent+ (,(sprintf #<<EOF
	// sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
	// refractory?
	if (S_.r_)
	  {
	    --S_.r_;
	  }
	else
	  {
           ~A
	  }
    
	// set new input current
	B_.I_stim_ = B_.currents_.get_value(lag);

	// log state data
	B_.logger_.record_data(current_steps + lag);

      }
EOF
  (if (lookup-def 'V_t defaults) threshold "")
  )))
       
  (pp indent  #\})

))

(define (output-buffers+nodes
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info
	 indent indent+)

  (let ((N (length state-index-map)))

    (case method


      ((leapfrog)

       (pp indent ,nl
	   (,(sprintf #<<EOF
~A::Buffers_::Buffers_(~A& n)
    : logger_(n),
      sys_(0),
      N(0),
      u(0),
      v(0),
      dvdt(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
  sysname sysname))

  (,(sprintf #<<EOF
~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
    : logger_(n),
      sys_(0),
      N(0),
      u(0),
      v(0),
      dvdt(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
  sysname sysname))
    ))

      ((cvode)

       (pp indent ,nl
	   (,(sprintf #<<EOF
~A::Buffers_::Buffers_(~A& n)
    : logger_(n),
      sys_(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
  sysname sysname))

  (,(sprintf #<<EOF
~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
    : logger_(n),
      sys_(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
  sysname sysname))
    ))

    ((gsl)

     (pp indent ,nl
	 (,(sprintf #<<EOF
~A::Buffers_::Buffers_(~A& n)
    : logger_(n),
      s_(0),
      c_(0),
      e_(0),
      N(0),
      u(0),
      jac(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
  sysname sysname))

  (,(sprintf #<<EOF
~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
    : logger_(n),
      s_(0),
      c_(0),
      e_(0),
      N(0),
      u(0),
      jac(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
  sysname sysname))
    ))


   )



  (pp indent  ,nl  (,(sprintf #<<EOF
~A::~A()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}
EOF
  sysname sysname)))

  (pp indent  ,nl  (,(sprintf #<<EOF
~A::~A(const ~A& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}
EOF
  sysname sysname sysname)))



  (pp indent (,(s+ sysname "::~" sysname "()")))

  (case method

   ((leapfrog)
    (pp indent #<<EOF
{
    if ( B_.u ) free(B_.u);
    if ( B_.v ) free(B_.v);
    if ( B_.dvdt ) free(B_.dvdt);
}
EOF
  ))


   ((gsl)
    (pp indent #<<EOF
{
    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
    if ( B_.s_ ) gsl_odeiv2_step_free(B_.s_);
    if ( B_.c_ ) gsl_odeiv2_control_free(B_.c_);
    if ( B_.e_ ) gsl_odeiv2_evolve_free(B_.e_);
    if ( B_.u ) free(B_.u);
    if ( B_.jac ) free(B_.jac);
}
EOF
  ))

   ((cvode)
    (pp indent #<<EOF
{

  if ( B_.sys_ )
  {
    /* Free y vector */
    N_VDestroy_Serial(B_.y);

    /* Free integrator memory */
//    CVodeFree(&B_.sys_);
  }
}
EOF
    ))
   )


  (pp indent ,nl (,(sprintf #<<EOF
  void ~A::init_node_(const Node& proto)
{
    const ~A& pr = downcast<~A>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}
EOF
  sysname sysname sysname)))

  (pp indent  ,nl (,(sprintf #<<EOF
void ~A::init_state_(const Node& proto)
{
    const ~A& pr = downcast<~A>(proto);
    S_ = State_(pr.P_);
}
EOF
  sysname sysname sysname)))


  (pp indent ,nl (,(sprintf #<<EOF
void ~A::init_buffers_()
{
EOF
 sysname)))

 (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
   (for-each (lambda (psc)
	       (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
       pscs))

       (pp indent+ (#<<EOF
    B_.currents_.clear();           
    Archiving_Node::clear_history();

    B_.logger_.reset();

    B_.step_ = Time::get_resolution().get_ms();
    B_.IntegrationStep_ = B_.step_;

    B_.I_stim_ = 0.0;
EOF
))


    (case method
     ((leapfrog)

      (pp indent+ (,(sprintf #<<EOF

    B_.N = ~A;
    B_.sys_ = ~A_dynamics;
    B_.u = (double *)malloc(sizeof(double) * B_.N);
    assert (B_.u);
    B_.v = (double *)malloc(sizeof(double) * B_.N);
    assert (B_.v);
    B_.dvdt = (double *)malloc(sizeof(double) * B_.N);
    assert (B_.dvdt);

EOF
N sysname))))

     ((cvode)

      (pp indent+ (,(sprintf #<<EOF

    int status, N;

    N = ~A;

    /* Creates serial vector of length N */
    B_.y = N_VNew_Serial(N);
    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) abort();

    for (int i = 0; i < N; i++)
    {
       Ith(B_.y,i) = S_.y_[i];
    }
 
    /* Calls CVodeCreate to create the solver memory and specify the 
     * Backward Differentiation Formula and the use of a Newton iteration */
    B_.sys_ = CVodeCreate(CV_BDF, CV_NEWTON);
    if (check_flag((void *)B_.sys_, "CVodeCreate", 0)) abort();

   /* Calls CVodeInit to initialize the integrator memory and specify the
    * right hand side function in y'=f(t,y), the initial time, and
    * the initial values. */
    status = CVodeInit (B_.sys_, ~A_dynamics, 0.0, B_.y);
    if (check_flag(&status, "CVodeInit", 1)) abort();

    /* Sets the relative and absolute error tolerances of CVode  */
    status = CVodeSStolerances (B_.sys_, 1e-7, 1e-7);
    if (check_flag(&status, "CVodeSStolerances", 1)) abort();

    /* Turn on CVode stability limit detection (only applicable for order 3 and greater) */
    status = CVodeSetStabLimDet (B_.sys_,true);
    if (check_flag(&status, "CVodeSetStabLimDet", 1)) abort();

    /* Sets the maximum order of CVode  */
    status = CVodeSetMaxOrd (B_.sys_,4);
    if (check_flag(&status, "CVodeSetMaxOrd", 1)) abort();

    /* Sets minimum and maximum step size. */
    //    status = CVodeSetMinStep (B_.sys_,1e-6);
    //    if (check_flag(&status, "CVodeSetMinStep", 1)) abort();
    status = CVodeSetMaxStep (B_.sys_,B_.step_);
    if (check_flag(&status, "CVodeSetMaxStep", 1)) abort();

   /* Calls CVodeSetUserData to configure the integrator to pass the 
    * params structure to the right-hand function */
   status = CVodeSetUserData(B_.sys_, reinterpret_cast<void*>(this));
   if (check_flag(&status, "CVodeSetUserData", 1)) abort();

    /* Initializes diagonal linear solver. */
    status = CVDiag (B_.sys_);
    if (check_flag(&status, "CVDiag", 1)) abort();
EOF
    N sysname))))

    ((gsl)

    (pp indent+ (,(sprintf 
#<<EOF
    static const gsl_odeiv2_step_type* T1 = gsl_odeiv2_step_rk2;
    B_.N = ~A;
  
    if ( B_.s_ == 0 )
      B_.s_ = gsl_odeiv2_step_alloc (T1, B_.N);
    else 
      gsl_odeiv2_step_reset(B_.s_);
    
    if ( B_.c_ == 0 )  
      B_.c_ = gsl_odeiv2_control_standard_new (1e-6, 1e-6, 1.0, 0.0);
    else
      gsl_odeiv2_control_init(B_.c_, 1e-6, 1e-6, 1.0, 0.0);
    
    if ( B_.e_ == 0 )  
      B_.e_ = gsl_odeiv2_evolve_alloc(B_.N);
    else 
      gsl_odeiv2_evolve_reset(B_.e_);
  
    B_.sys_.function  = ~A_dynamics; 
    B_.sys_.jacobian  = ~A_jacobian;
    B_.sys_.dimension = B_.N;
    B_.sys_.params    = reinterpret_cast<void*>(this);

    B_.u = (double *)malloc(sizeof(double) * B_.N);
    assert (B_.u);
    B_.jac = (double *)malloc(sizeof(double) * B_.N);
    assert (B_.jac);

EOF
  N sysname sysname))))
  )
  (pp indent ,nl #\})

  (let ((iv (lookup-def 'v state-index-map)))
    (if iv
	(pp indent ,nl (,(sprintf #<<EOF
void ~A::calibrate()
{
    B_.logger_.init();  
    V_.RefractoryCounts_ = 20;
    V_.U_old_ = S_.y_[~A];
}
EOF
  sysname iv)))
  ))
 ))



(define (nemo:nest-translator sys . rest)

  (define (cid x)  (second x))
  (define (cn x)   (first x))

  (let-optionals rest ((dirname ".") (method #f))

    (let ((method (or method 'gsl)))

      (if (not (member method '(gsl cvode leapfrog)))
	  (nemo:error 'nemo:nest-translator ": unknown method " method))

  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
    (let ((imports  ((dis 'imports)  sys))
	  (exports  ((dis 'exports)  sys)))
      (let* ((indent      0)
	     (indent+     (+ 2 indent ))

	     (sysname     (nest-name ((dis 'sysname) sys)))
	     (prefix      (->string sysname))
	     (deps*       ((dis 'depgraph*)  sys))
	     (consts      ((dis 'consts)     sys))
	     (asgns       ((dis 'asgns)      sys))
	     (states      ((dis 'states)     sys))
	     (reactions   ((dis 'reactions)  sys))
	     (defuns      ((dis 'defuns)     sys))
	     (components  ((dis 'components) sys))
	     
	     (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
	     (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))

	     (const-defs       (filter-map
				(lambda (nv)
				  (and (not (member (first nv) nest-builtin-consts))
				       (let ((v1 (canonicalize-expr/C++ (second nv))))
					 (list (nest-name (first nv)) v1))))
				consts))
	     
	     (gate-complex-info    (nemo:gate-complex-query sys))
	     (synapse-info         (nemo:post-synaptic-conductance-query sys))
	     (i-synapses           (lookup-def 'i-synapses synapse-info))
	     (defaults             (nemo:defaults-query sys))
             (geometry             (nemo:geometry-query sys))

	     (gate-complexes  (lookup-def 'gate-complexes gate-complex-info))
	     (perm-ions       (map (match-lambda ((comp i e erev val) `(,comp ,(nest-name i) ,(nest-name e) ,erev)))
				   (lookup-def 'perm-ions gate-complex-info)))
	     (acc-ions        (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out)))))
				   (lookup-def 'acc-ions gate-complex-info)))
	     (epools          (lookup-def 'pool-ions gate-complex-info))
	     (pool-ions       (pool-ion-name-map nest-name  epools))

	     (i-gates         (lookup-def 'i-gates gate-complex-info))

	     (comprc         (any (match-lambda ((name 'membrane-tau id) (list name id)) (else #f)) components))
	     (compcap        (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
	     (mrc            (or (and comprc (car ((dis 'component-exports) sys (cid comprc))))
				 (and compcap (car ((dis 'component-exports) sys (cid compcap))))))

	     (soma-geometry  (lookup-def 'soma geometry))
	     (marea          (and soma-geometry (third soma-geometry)))
		

	     (i-eqs (filter-map
		     (lambda (gate-complex) 
		       
		       (let* ((label             (first gate-complex))
			      (n                 (second gate-complex))
			      (subcomps          ((dis 'component-subcomps) sys n))
			      (acc               (lookup-def 'accumulating-substance subcomps))
			      (perm              (lookup-def 'permeating-ion subcomps))
			      (permqs            (and perm ((dis 'component-exports) sys (cid perm))))
			      (pore              (lookup-def 'pore subcomps))
			      (permeability      (lookup-def 'permeability subcomps))
			      (gates             (filter (lambda (x) (equal? (car x) 'gate)) subcomps))
			      (sts               (map (lambda (gate) ((dis 'component-exports) sys (cid gate))) gates))
			      )


                         (if (and pore (null? permqs))
                             (nemo:error 'nemo:nest-translator ": ion channel definition " label
                                         "permeating-ion component lacks exported quantities"))

			 (for-each 
			  (lambda (st)
			    (if (null? st)
				(nemo:error 'nemo:nest-translator: "ion channel definition " label
					    "gate component lacks exported quantities")))
			  sts)
			 
			 (if (not (or pore permeability))
			     (nemo:error 'nemo:nest-translator ": ion channel definition " label
					 "lacks any pore or permeability components"))


			 (cond ((and perm permeability (pair? gates))
				     (let* ((i     (nest-name (s+ 'i (cn perm))))
					    (pmax  (car ((dis 'component-exports) sys (cid permeability))))
					    (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
					    (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
					    (gion  `(* ,pmax ,(sum (map (lambda (gpwr) 
									  (match gpwr ((x)  x) (else `(* ,@gpwr))))
									gpwrs))))
					    )
					 (list i #f gion (nest-name (s+ 'i_ label) ))))

			       ((and perm pore (pair? gates))
				(case (cn perm)
				  ((non-specific)
				   (let* ((i     (nest-name 'i))
					  (e     (car permqs))
					  (gmax  (car ((dis 'component-exports) sys (cid pore))))
					  (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
					  (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
					  (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
									(match gpwr ((x)  x) (else `(* ,@gpwr))))
								      gpwrs))))
					  )
				     (list i e gion  (nest-name (s+ 'i_ label) ))))

				  (else
				   (let* ((i     (nest-name (s+ 'i (cn perm))))
					  (e     (car permqs))
					  (gmax  (car ((dis 'component-exports) sys (cid pore))))
					  (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
					  (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
					  (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
									(match gpwr ((x)  x) (else `(* ,@gpwr))))
								      gpwrs))))
					  )
				     (list i e gion (nest-name (s+ 'i_ label) ))))))
			       
			       ((and perm pore)
				(case (cn perm)
				  ((non-specific)
				   (let* ((i     (nest-name 'i))
					  (e     (car permqs))
					  (gmax  (car ((dis 'component-exports) sys (cid pore)))))
				     (list i e gmax (nest-name (s+ 'i_ label) ))))
				  (else
				   (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))))
			       
			       ((and acc pore (pair? gates))
				(let* ((i     (nest-name (s+ 'i (cn acc))))
				       (gmax  (car ((dis 'component-exports) sys (cid pore))))
				       (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
				       (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
				       (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
								     (match gpwr ((x)  x) (else `(* ,@gpwr))))
								   gpwrs))))
				       )
				  (list i #f gion  (nest-name (s+ 'i_ label) ))))
			       (else (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))
			       )))
		     gate-complexes))

	     (i-names (delete-duplicates (map first i-eqs)))
		
	     (i-eqs  (fold  (lambda (i-gate ax) 
			      (let ((i-gate-var (first i-gate)))
				(cons (list (nest-name 'i) #f i-gate-var (s+ 'i_ (second i-gate))) ax)))
			    i-eqs i-gates))

	     (i-eqs  (fold  (lambda (i-synapse ax) 
			      (cons (list (first i-synapse) (third i-synapse) (second i-synapse) (gensym 'i_syn )) ax))
			    i-eqs 
			    i-synapses))

	     (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
	     
	     (i-eqs  (fold (lambda (b ax) 
			     (match b 
				    ((and ps ((i e gion ii) . rst))  
				     (let loop ((ps ps) (summands (list)) (eqs (list)))
				       (if (null? ps)
					   
					   (let* ((sum0  (sum summands))
						  (sum1  (rhsexpr/C++ sum0))
						  (sum2  (add-params-to-fncall 
							  (canonicalize-expr/C++ sum1) builtin-fns)))
					     (append eqs (list (list i sum2)) ax))
					   
					   (match-let (((i e gion ii) (car ps)))
						      (loop (cdr ps) 
							    (cons ii summands) 
							    (let* ((expr0 (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
								   (expr1 (canonicalize-expr/C++ expr0)))
							      (cons (list ii expr1) eqs)))))))
				    
				    ((i e gion ii)
				     (let* ((expr0  (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
					    (expr1  (canonicalize-expr/C++ expr0)))
				       (cons (list i expr1) ax)))
				    
				    (else ax)))
			   (list) i-bkts))

	     (external-eq-defs   (external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++))

	     (asgn-eq-defs       (poset->asgn-eq-defs* poset sys nest-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
	     
	     (rate-eq-defs       (reverse (poset->rate-eq-defs* poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)))
	     
	     (reaction-eq-defs   (poset->reaction-eq-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++))
	     
	     (init-eq-defs       (poset->init-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
	     
	     (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
				      (poset->state-conserve-eq-defs poset sys nest-name nest-state-name)))
	     
	     (imports-sans-v (filter (lambda (x) (not (equal? 'v (first x)))) imports))

	     (v-eq    (and (not (null? i-names))
			   (let ((istim "(node.B_.I_stim_)" )) 
			     (cond
                              ((and mrc marea)
                               (list 'v (rhsexpr/C++ 
                                         `(* 1e3 (/ (+ (* ,istim (/ 100. ,marea)) 
                                                       (neg ,(sum i-names))) ,mrc)))))
                              (marea
                               (list 'v (rhsexpr/C++ 
                                         `(* 1e3 (+ (* ,istim (/ 100. ,marea))
                                                    (neg ,(sum i-names)))))))
                              (mrc
                               (list 'v (rhsexpr/C++ `(* 1e3 (/ (+ ,istim (neg ,(sum i-names))) ,mrc)))))
                              
                              (else
                               (list 'v (rhsexpr/C++ `(* 1e3 (+ ,istim (neg ,(sum i-names)))))))
                              ))
                           ))

	     
	     (state-index-map  (let ((acc (fold (lambda (def ax)
						  (let ((st-name (first def)))
						    (list (+ 1 (first ax)) 
							  (cons `(,st-name ,(first ax)) (second ax)))))
						(list 0 (list)) 
						(if (lookup-def 'v imports)
						    (cons (list 'v) rate-eq-defs)
						    rate-eq-defs)
						)))
				 (second acc)))
	     
	     (steady-state-index-map  (let ((acc (fold (lambda (def ax)
							 (let ((st-name (first def)))
							   (if (not (alist-ref st-name init-eq-defs))
							       (list (+ 1 (first ax)) 
								     (cons `(,st-name ,(first ax)) (second ax)))
							       ax)))
						       (list 0 (list)) 
						       rate-eq-defs)))
					(second acc)))
	     
	     (dfenv 
	      (map (lambda (x) (let ((n (first x))) (list n (nest-name (s+ "d_" n )))))
		   defuns))

	     )
	
	
	
	(for-each 
	 (lambda (a)
	   (let ((acc-ion   (car a)))
	     (if (assoc acc-ion perm-ions)
		 (nemo:error 'nemo:nest-translator 
			     ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
	 acc-ions)


	(let ((cpp-output  (open-output-file (make-output-fname dirname prefix ".cpp")))
	      (hpp-output  (open-output-file (make-output-fname dirname prefix ".h"))))
	  
	  ;; include files and other prelude definitions
	  (with-output-to-port cpp-output
	    (lambda ()
              (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
	      (output-prelude sysname steady-state-index-map indent)
	      ))
	  
	  (with-output-to-port cpp-output
	    (lambda () (pp indent "namespace nest {" ,nl)))
	  
	  ;; user-defined functions
	  (let ((define-fn  (make-define-fn sysname))
		(define-fn-header (make-define-fn-header sysname)))
	    
	    (for-each (lambda (fndef) 
			(if (not (member (car fndef) builtin-fns))
			    (with-output-to-port cpp-output
			      (lambda ()
				(apply define-fn-header (cons indent fndef))
				(pp indent ,nl)))
			    ))
		      defuns)

	    (for-each (lambda (fndef) 
			(if (not (member (car fndef) builtin-fns))
			    (with-output-to-port cpp-output
			      (lambda ()
				(apply define-fn (cons indent fndef))
				(pp indent ,nl)))
			    ))
		      defuns)
	    )
	  
	  
	  ;; derivative function
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-dy sysname method imports-sans-v const-defs state-index-map 
			 external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
			 pool-ions i-eqs v-eq
			 indent indent+)
	      ))

	  ;; Jacobian function
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-jac  sysname method state-index-map pool-ions i-eqs v-eq indent indent+)
	      ))

	  ;;  system recordables
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-recordables
	       sysname state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; system parameters
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-parameters sysname imports-sans-v const-defs 
				 external-eq-defs defaults indent indent+)
	      ))
	  
	  ;; initial values function
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-init sysname state-index-map steady-state-index-map 
			      imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
			      reaction-eq-defs i-eqs pool-ions perm-ions defaults indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; accessors and modifiers for states and parameters
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-accessors+modifiers
	       sysname imports-sans-v state-index-map 
	       const-defs asgn-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions 
	       indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  
	  ;; buffer and node initialization
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-buffers+nodes
	       sysname method state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions 
	       synapse-info
	       indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; state update
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-update
	       sysname method state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs defaults i-eqs pool-ions perm-ions 
	       synapse-info indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; spike handler
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-handle
	       sysname state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions 
	       synapse-info
	       indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  (with-output-to-port cpp-output
	    (lambda () (pp indent "}" ,nl)))


	     (with-output-to-port hpp-output
	       (lambda ()
                 (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
		 (output-header
		  sysname method state-index-map steady-state-index-map defaults
		  external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
		  reaction-eq-defs i-eqs pool-ions perm-ions 
		  synapse-info
		  indent indent+)
		 (pp indent ,nl)
		 ))
	     
	     (close-output-port cpp-output)
	     (close-output-port hpp-output)
	   
	   ))
      ))
  ))
)
)