;;       
;; 
;; An extension for translating NEMO models to NEST code.
;;
;; Copyright 2011-2012 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 lolevel ports extras srfi-1 srfi-13)
	
	(require-extension lolevel matchable strictly-pretty environments
			   varsubst datatype nemo-core nemo-utils
			   nemo-gate-complex)

(define nest-builtin-consts
  `())

(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
      fpvector-ref))

(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)
			    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:empty) (doc:empty)))
(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:connect (doc:group (doc:connect (doc:text "if") c))
		   (doc:connect (doc:nest 2 e1)
				(doc:nest 2 (doc:connect (doc:text "else") e2))))
      (doc:text "end"))))

(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 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 ((pow)  '^)  ((ln) 'log) (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)))
			
			((^)  
			 (format-op/C++ indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/C++ indent+ x)))
						   (if (or (symbol? x)  (number? x)) fx
						       (if (or mdiv? 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 (state-init n init)
  (let* ((init  (rhsexpr/C++ init))
	 (init1 (canonicalize-expr/C++ init)))
    (list (nest-name n) init1)))


(define (asgn-eq n rhs)
  (let* ((fbody   (rhsexpr/C++ rhs))
	 (fbody1  (canonicalize-expr/C++ fbody)))
    (list (nest-name n) fbody1)))


(define (reaction-eq n open transitions conserve)
  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nest-state-name)))
    (let ((nodes ((g 'nodes))))
      (list (nest-name n) (nest-state-name n open)))))


(define (reaction-transition-eqs n initial open transitions conserve power)
  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nest-state-name)))
     (let* ((out-edges  (g 'out-edges))
	    (in-edges   (g 'in-edges))
	    (nodes      ((g 'nodes))))
       ;; generate differential equations for each state in the transitions system
      (let ((eqs    (fold (lambda (s ax) 
			    (if (= (first cnode) (first s) ) ax
				(let* ((out   (out-edges (first s)))
				       (in    (in-edges (first s)))
				       (open? (eq? (second s) open))
				       (name  (nest-name (lookup-def (second s) node-subs))))
				  (let* ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
						       `(+ (neg ,(sum (map third out)))
							   ,(sum (map third in))))
						      ((and (not (null? out)) (null? in))
						       `(neg ,(sum (map third out))))
						      ((and (null? out) (not (null? in)))
						       (sum (map third in)))))
					 (fbody0 (rhsexpr/C++ rhs1)))
				    
				    (cons (list name (canonicalize-expr/C++ fbody0)) ax)
				    ))))
			  (list) nodes)))
	eqs))))
	    
	

(define (poset->asgn-eq-defs poset sys)
  (fold-right
   (lambda (lst ax)
     (fold  (lambda (x ax) 
	      (match-let (((i . n)  x))
			 (let ((en (environment-ref sys n)))
			   (if (nemo:quantity? en)
			       (cases nemo:quantity en
				      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
				      (else  ax))
			       ax))))
	    ax lst))
   (list) poset))

(define (poset->rate-eq-defs poset sys)
  (fold-right
   (lambda (lst ax)
     (fold  (lambda (x ax) 
	      (match-let (((i . n)  x))
			 (let ((en (environment-ref sys n)))
			   (if (nemo:quantity? en)
			       (cases nemo:quantity en

				      (REACTION  (name initial open transitions conserve power) 
						 (append (reaction-transition-eqs name initial open transitions 
										  conserve power) ax))
				      
				      (RATE (name initial rhs power)
					    (let ((fbody0  (rhsexpr/C++ rhs))
						  (dy      (nest-name name) ))
					      (cons (list dy (canonicalize-expr/C++ fbody0)) ax)))

				      (else  ax))
			       ax))))
	    ax lst))
   (list) poset))


(define (poset->reaction-eq-defs poset sys)
  (fold-right
   (lambda (lst ax)
     (fold  (lambda (x ax) 
	      (match-let (((i . n)  x))
			 (let ((en (environment-ref sys n)))
			   (if (nemo:quantity? en)
			       (cases nemo:quantity en
				      (REACTION  (name initial open transitions conserve power) 
						 (cons (reaction-eq name open transitions conserve) ax))
				      (else  ax))
			       ax))))
	    ax lst))
   (list) poset))


(define (poset->init-defs poset sys)
  (fold-right
   (lambda (lst ax)
     (fold-right
      (lambda (x ax) 
	      (match-let (((i . n)  x))
			 (let ((en (environment-ref sys n)))
			   (if (nemo:quantity? en)
			       (cases nemo:quantity en
				      (REACTION  (name initial open transitions conserve power)
						 (if (nemo:rhs? initial)
						     (cons* (state-init name initial) 
							    (state-init (nest-state-name name open) name) ax) 
						     ax))

				      (RATE  (name initial rhs power)
					     (if (and initial (nemo:rhs? initial))
						 (cons (state-init name initial) ax) 
						 ax))

				      (else  ax))
			       ax))))
	    ax lst))
   (list) poset))


(define (poset->state-conserve-eq-defs poset sys)
  (fold-right
   (lambda (lst ax)
     (fold  (lambda (x ax) 
	      (match-let (((i . n)  x))
			 (let ((en (environment-ref sys n)))
			   (if (nemo:quantity? en)
			       (cases nemo:quantity en
				      (REACTION (name initial open transitions conserve power)
						(if (and (list? conserve) (every nemo:conseq? conserve))
						    (cons (state-conseqs (nest-name name) transitions conserve
									nest-state-name) ax) 
						    ax))
				      (else  ax))
			       ax))))
	    ax lst))
   (list) poset))


(define (find-locals defs)
  (concatenate
   (map (lambda (def) 
	  (match def 
		 (('let bnds body) 
		  (let ((bexprs (map second bnds)))
		    (concatenate (list (map first bnds) 
				       (find-locals bexprs )
				       (find-locals (list body))))))
		 (('if c t e)      (append (find-locals (list t)) (find-locals (list e))))
		 ((s . rest)       (find-locals rest))
		 (else (list)))) 
	defs)))


(define (reaction-power sys n)
  (let ((en (environment-ref sys n)))
    (if (nemo:quantity? en)
	(cases nemo:quantity en
	       (REACTION  (name initial open transitions conserve power)  power)
	       (else  #f))  #f)))


(define (bucket-partition p lst)
  (let loop ((lst lst) (ax (list)))
    (if (null? lst) ax
	(let ((x (car lst)))
	  (let bkt-loop ((old-bkts ax) (new-bkts (list)))
	    (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
		(if (p x (caar old-bkts ))
		    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
		    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))


(define (make-define-fn )
  (lambda (indent n proc)
    (let ((lst (procedure-data proc))
	  (indent+ (+ 2 indent)))
      (let ((rt       (lookup-def 'rt lst))
	    (formals  (lookup-def 'formals lst))
	    (vars     (lookup-def 'vars lst))
	    (body     (lookup-def 'body lst))
	    (rv       (gensym 'rv)))
	(pp indent ,nl (,rt ,(nest-name n) (,(slp ", " vars)) "{" ))
	(let* ((body0 (rhsexpr/C++ body))
	       (body1 (canonicalize-expr/C++ body0))
	       (lbs   (enum-bnds body1 (list))))
	  (pp indent+ (double ,rv ";"))
	  (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs) ";")))
	  (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
	  (pp indent+ ,(s+ "return " rv ";"))
	  )
	(pp indent "}")))  
    ))


(define (output-dy sysname parameters state-index-map 
		   rate-eq-defs reaction-eq-defs asgn-eq-defs
		   pool-ions mcap i-eqs v-eq indent indent+)


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

  (pp indent+ ,nl
      ("// S is shorthand for the type that describes the model state ")
      (typedef ,(s+ sysname "::state_t") "S;")

      ,nl

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

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

  ;; TODO parameters, time
	   

  (let* (
	 (asgn-eq-defs 
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   asgn-eq-defs))

	 (rate-eq-defs 
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   rate-eq-defs))

	 (reaction-eq-defs 
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   reaction-eq-defs))

	 (i-eqs 
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   i-eqs))

	 (v-eq
	  (and v-eq
	       (list (first v-eq) (canonicalize-expr/C++ (second v-eq)))))

	 (eqs 
	  (append 
	   
	   asgn-eq-defs
	   reaction-eq-defs
		 
	   (map (lambda (pool-ion)
		  (let ((n (third pool-ion))
			(b (first pool-ion)))
		    (list n b)))
		pool-ions)))
	 
	 (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)))))
	 )

    (if (not (null? eq-locals)) 
	(pp indent+ (double ,(slp ", " eq-locals) ";")))
    
    (for-each (lambda (def)
		(let ((n (first def)) )
		  (pp indent+
		      ,(s+ "// rate equation " n)
		      ,(expr->string/C++ 
			(s+ "y[" (lookup-def n state-index-map) "]") n))))
	      rate-eq-defs)

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

    (for-each (lambda (def)
		(let ((n (s+ "f[" (lookup-def (first def) state-index-map) "]")) 
		      (b (second def)))
		  (pp indent+ ,(s+ "// state " (first def))
		      ,(expr->string/C++ b n))))
	      rate-eq-defs)
    
    (for-each (lambda (def) 
		(pp indent+ ,(expr->string/C++ (second def) (first def)))) 
	    i-eqs)
    
    (if v-eq
	(let ((vi (lookup-def 'v state-index-map)))
	  (if vi
	      (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "f[" vi "]"))))))
    
    (pp indent+ ,nl ("return GSL_SUCCESS;"))
    (pp indent  #\})

    ))



(define (output-parameters sysname const-defs indent indent+)

  (let* ((parameter-defs
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   const-defs))
	 (parameter-locals  (find-locals (map second parameter-defs))))

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

      (pp indent+ ,const-exprs)
      
      (pp indent ("{}"))
      
      )))


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

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

	 (asgn-eq-defs
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   asgn-eq-defs))
	 
	 (init-eq-defs
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   init-eq-defs))

	 (reaction-eq-defs 
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   reaction-eq-defs))

	 (i-eqs 
	  (map 
	   (lambda (def) (list (first def) (canonicalize-expr/C++ (second def))))
	   i-eqs))

	 (v-eq
	  (and v-eq
	       (list (first v-eq) (canonicalize-expr/C++ (second v-eq)))))
	 
	 (init-eqs 
	  (append 
	   
	   asgn-eq-defs
	   init-eq-defs
	   
	   (map (lambda (pool-ion)
		  (let ((n (third pool-ion))
			(b (first 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  #\{)

    (if (not (null? init-locals)) 
	(pp indent+ (double ,(slp ", " init-locals) ";")))

    (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)

    (if v-eq
	(let ((vi (lookup-def 'v state-index-map)))
	  (if vi
	      (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "y_[" vi "]"))))))
      
    (pp indent  #\})
    
    (pp indent ,nl (,(s+ sysname "::State_::State_") (const Parameters_& p) ": r_(s.r_)"))
    (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];
     r_ = s.r_;
     return *this;
EOF
N)))

     (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_f;
  F.df  = &~A_df;
  F.fdf = &~A_fdf;
  F.n   = ~A;
  F.params = &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 sysname N N N)))
)



(define (output-accessors+modifiers
	 sysname state-index-map const-defs asgn-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs v-eq 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, names::~A, ~A);" n n)))))
   const-defs)

  (pp indent  #\})

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

  (for-each
   (lambda (def)
     (let ((n (first def)))
       (pp indent+ (,(sprintf "updateValue<double_t>(d, names::~A, ~A);" 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 (lookup-def (second def) state-index-map)))
       (pp indent+ (,(sprintf "def<double_t>(d, names::~A, y_[~A]);" n i)))))
   reaction-eq-defs)

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

  (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 (lookup-def (second def) state-index-map)))
       (pp indent+ (,(sprintf "updateValue<double_t>(d, names::~A, y_[~A]);" n i)))))
   reaction-eq-defs)

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

  (pp indent  #\})

)

(define (output-buffers+nodes
	 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+)

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

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

  (pp indent  ,nl  (,(sprintf #<<EOF
~A::Buffers_::Buffers_(const Buffers_&, ~A& n)
    : logger_(n),
      s_(0),
      c_(0),
      e_(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 "()")))
  (pp indent #<<EOF
{
    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
    if ( B_.s_ ) gsl_odeiv_step_free(B_.s_);
    if ( B_.c_ ) gsl_odeiv_control_free(B_.c_);
    if ( B_.e_ ) gsl_odeiv_evolve_free(B_.e_);
}
EOF
  )

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

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

  (pp indent ,nl (,(sprintf #<<EOF
void ~A::init_buffers_()
{
    B_.spike_exc_.clear();          
    B_.spike_inh_.clear();          
    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;

    static const gsl_odeiv_step_type* T1 = gsl_odeiv_step_rkf45;
  
    if ( B_.s_ == 0 )
      B_.s_ = gsl_odeiv_step_alloc (T1, ~A);
    else 
      gsl_odeiv_step_reset(B_.s_);
    
    if ( B_.c_ == 0 )  
      B_.c_ = gsl_odeiv_control_y_new (1e-3, 0.0);
    else
      gsl_odeiv_control_init(B_.c_, 1e-3, 0.0, 1.0, 0.0);
    
    if ( B_.e_ == 0 )  
      B_.e_ = gsl_odeiv_evolve_alloc(~A);
    else 
      gsl_odeiv_evolve_reset(B_.e_);
  
    B_.sys_.function  = ~A_dynamics; 
    B_.sys_.jacobian  = 0;
    B_.sys_.dimension = ~A;
    B_.sys_.params    = reinterpret_cast<void*>(this);
}
EOF
  sysname N N sysname N)))

  (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 (output-update
	 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+)
  
  (let ((vi (lookup-def 'v state-index-map)))
    
    (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);

    for ( long_t lag = from ; lag < to ; ++lag )
      {
    
	double tt = 0.0 ; //it's all relative!
	V_.U_old_ = S_.y_[~A];

   
	// adaptive step integration
	while (tt < B_.step_)
	{
	  const int status = gsl_odeiv_evolve_apply(B_.e_, B_.c_, B_.s_, 
				 &B_.sys_,              // system of ODE
				 &tt,                   // from t...
				  B_.step_,             // ...to t=t+h
				 &B_.IntegrationStep_ , // integration window (written on!)
				  S_.y_);	        // neuron state

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

	// sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
	// refractory?
	if (S_.r_)
	  {
	    --S_.r_;
	  }
	else
	  {
	    // (threshold && maximum)
	    if (S_.y_[~A] >= P_.V_T + 30. && 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);
	      }
	  }
    
	// set new input current
	B_.I_stim_ = B_.currents_.get_value(lag);

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

      }
EOF
  vi vi 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 indent indent+)
  
  (pp indent  ,nl (,(sprintf #<<EOF
void ~A::handle(SpikeEvent & e)
  {
    assert(e.get_delay() > 0);

    if(e.get_weight() > 0.0)
      {
	B_.spike_exc_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
				e.get_weight() * e.get_multiplicity());
      }
    else
      {
	// add with negative weight, ie positive value, since we are changing a
	// conductance
	B_.spike_inh_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
				-e.get_weight() * e.get_multiplicity());
      }
  }
EOF
  sysname)))

  (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_(names::~A, &~A::get_y_elem_<~A>);" 
			      n sysname i)))))
   state-index-map)
  (pp indent #\})
  )


(define (output-prelude sysname 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 <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_exp.h>
"  sysname))))


(define (output-header 
	 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 ("#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\"
#include <gsl/gsl_odeiv.h>
"))
  
  (pp indent (extern "\"C\""
		     int ,(s+ sysname "_dynamics")
		     "(double, const double*, double*, void*)" #\;) ,nl)

  (pp indent 
      (,(sprintf "class ~A:" sysname)
       "public Archiving_Node { "
       "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
))

  (pp indent (,(sprintf #<<EOF
    // make dynamics function quasi-member
    friend int ~A_dynamics(double, const double*, double*, void*);

    // 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 sysname)))

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

  (for-each 
   (lambda (x) (pp indent+ (,(sprintf "double ~A;" (first x)))))
   const-defs)

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

  (pp indent+ "};")

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

  (pp indent+ ,nl 
      (,(sprintf "double y_[~A];" (length state-index-map)))
      "int_t     r_;"
      "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))
#<<EOF
UniversalDataLogger<hh_cond_exp_traub> logger_;

  RingBuffer spike_exc_;
  RingBuffer spike_inh_;
  RingBuffer currents_;

  gsl_odeiv_step*    s_;    //!< stepping function
  gsl_odeiv_control* c_;    //!< adaptive stepsize control function
  gsl_odeiv_evolve*  e_;    //!< evolution function
  gsl_odeiv_system   sys_;  //!< struct describing system

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

  /** 
   * 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))
      "}; ")

  )

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

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

  (let-optionals rest ((dirname "."))
  (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
    (let ((imports  ((dis 'imports)  sys))
	  (exports  ((dis 'exports)  sys)))
      (let* ((indent      0)
	     (indent+     (+ 2 indent ))
	     (eval-const  (dis 'eval-const))
	     (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))

	     (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
	     (perm-ions     (map (match-lambda ((comp i e erev) `(,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     (map (lambda (lst) (map nest-name lst)) epools))

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

	     (capcomp       (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
	     (mcap          (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
		

	     (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))
			      (gate              (lookup-def 'gate subcomps))
			      (sts               (and gate ((dis 'component-exports) sys (cid gate)))))
			 
			 (if (not (or pore permeability))
			     (nemo:error 'nemo:nest-translator ": ion channel definition " label
					 "lacks any pore or permeability components"))


			 (cond ((and perm permeability gate)
				     (let* ((i     (nest-name (s+ 'i (cn perm))))
					    (pmax  (car ((dis 'component-exports) sys (cid permeability))))
					    (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
					    (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
					    (gion  `(* ,pmax ,@sptms)))
					 (list i #f gion (nest-name (s+ 'i_ label) ))))

			       ((and perm pore gate)
				(case (cn perm)
				  ((non-specific)
				   (let* ((i     (nest-name 'i))
					  (e     (car permqs))
					  (gmax  (car ((dis 'component-exports) sys (cid pore))))
					  (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
					  (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
					  (gion  `(* ,gmax ,@sptms)))
				     (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 (n) (reaction-power sys n)) sts))
					  (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
					  (gion  `(* ,gmax ,@sptms)))
				     (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 gate)
				(let* ((i     (nest-name (s+ 'i (cn acc))))
				       (gmax  (car ((dis 'component-exports) sys (cid pore))))
				       (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
				       (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
				       (gion  `(* ,gmax ,@sptms)))
				  (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-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  (canonicalize-expr/C++ sum1)))
					     (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))

	     (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
	     
	     (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys)))
	     
	     (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
	     
	     (init-eq-defs       (poset->init-defs poset sys))
	     
	     (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
				      (poset->state-conserve-eq-defs poset sys)))
	     
	     (globals          (map nest-name 
				    (delete-duplicates (append 
							exports
							(map second  perm-ions)
							(map third   perm-ions)
							(map second  acc-ions)
							(map third   acc-ions)
							(map fourth  acc-ions)
							(map second  pool-ions)
							(map third   pool-ions)
							(map first   imports) 
							(map first   const-defs)
							(map first   asgn-eq-defs)
							(map (lambda (gate-complex) (nest-name (s+ 'i_ (first gate-complex)))) gate-complexes )
							(map (lambda (i-gate) (nest-name (s+ 'i_ (second i-gate)))) i-gates )
							))))

	     (parameters (map nest-name (map first const-defs)))

	     (v-eq    (if (and mcap (member 'v globals))
			  (list 'v (rhsexpr/C++ `(/ (neg ,(sum i-names)) ,mcap)))
			  (list 'v 0.0)))
	     
	     (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 (member 'v globals)
						    (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)

	   (for-each 
	    (lambda (p)
	      (let ((pool-ion  (car p)))
		(if (assoc pool-ion perm-ions)
		    (nemo:error 'nemo:nest-translator 
				": ion species " pool-ion " cannot be declared as both pool and permeating"))))
	    pool-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 ()
		 (output-prelude sysname indent)
		 ))

	     ;; user-defined functions
	     (let ((define-fn  (make-define-fn)))
	       (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))

	     ;;  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)
		 ))

	     
	     ;; derivative function
	     (with-output-to-port cpp-output
	       (lambda ()
		 (output-dy sysname globals state-index-map 
			    rate-eq-defs reaction-eq-defs asgn-eq-defs
			    pool-ions mcap i-eqs v-eq
			    indent indent+)
		 ))

	     ;; system parameters
	     (with-output-to-port cpp-output
	       (lambda ()
		 (output-parameters sysname const-defs
				    indent indent+)
		 ))

	     ;; initial values function
	     (with-output-to-port cpp-output
	       (lambda ()
		 (output-init sysname state-index-map steady-state-index-map 
			      const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
			      reaction-eq-defs i-eqs v-eq
			      pool-ions perm-ions indent indent+)
		 (pp indent ,nl)
		 ))

	     ;; accessors and modifiers for states and parameters
	     (with-output-to-port cpp-output
	       (lambda ()
		 (output-accessors+modifiers
		  sysname state-index-map 
		  const-defs asgn-eq-defs rate-eq-defs 
		  reaction-eq-defs i-eqs v-eq pool-ions perm-ions 
		  indent indent+)
		 (pp indent ,nl)
		 ))

	     
	     ;; buffer and node initialization
	     (with-output-to-port cpp-output
	       (lambda ()
		 (output-buffers+nodes
		  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)
		 ))

	     ;; state update
	     (with-output-to-port cpp-output
	       (lambda ()
		 (output-update
		  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)
		 ))

	     ;; 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 indent indent+)
		 (pp indent ,nl)
		 ))


	     (with-output-to-port hpp-output
	       (lambda ()
		 (output-header
		  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)
		 ))

	     
	     (close-output-port cpp-output)
	     (close-output-port hpp-output)
	   
	   ))
      ))
  ))


)