;;       
;; TODO: implement a new component type to be used for thresholds
;; TODO: finish gsl-fminimizer
;; 
;; 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 srfi-69)
	
	(require-extension lolevel matchable strictly-pretty 
			   varsubst datatype nemo-core nemo-utils
			   nemo-gate-complex nemo-synapse)

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

(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 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) (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 (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 (hash-table-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 (hash-table-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 (hash-table-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 (hash-table-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 (hash-table-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 (rate/reaction-power sys n)
  (let ((en (hash-table-ref sys n)))
    (if (nemo:quantity? en)
	(cases nemo:quantity en
	       (REACTION  (name initial open transitions conserve power)  
			  power)
	       (RATE    (name initial rhs 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 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   (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))
		'("void* pnode"))))
	  (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) "{" ))
	  (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) ";")))
	    (if (not (null? consts)) 
		(begin (pp indent+ (double ,(slp ", " (map nest-name consts)) ";")
			   ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));"))
		       (for-each (lambda (x) (pp indent+ (,(nest-name x) = ,(s+ "node.P_." (nest-name x) ";")))) consts)
		       ))
	    (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
	    (pp indent+ ,(s+ "return " rv ";"))
	    ))
	(pp indent "}"))
    )))


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


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


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



    (pp indent+ (double ,(slp ", " (delete-duplicates 
				    (map ->string
					 (append 
					  eq-locals
					  eq-order
					  (map first state-index-map)
					  (map first i-eqs)
					  (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));")
	
	,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+ "node.P_." n) n))))
	      const-defs)

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

    (for-each (lambda (def)
		(let ((n (first def)) )
		  (pp indent+
		      ,(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+ 
			    ,(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+ "// rate equation " (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  #\{)

    (pp indent+ (double ,(slp ", " (delete-duplicates
				    (map ->string 
					 (append 
					  init-locals
					  init-order
					  (map first state-index-map) 
					  (map first i-eqs)
					  (map first const-defs)
					  (map first reaction-eq-defs)
					  ))
				    string=?))
			      ";"))

    (for-each (lambda (def)
		(let ((n (first def)) )
		  (pp indent+
		      ,(expr->string/C++ 
			(s+ "p." n) n))))
	      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)

    (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 State_& s) ": 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_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 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, ~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)

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


  (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 
	 synapse-info
	 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 #<<EOF
void ~A::init_buffers_()
{
EOF
)
     (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+ (,(sprintf #<<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;

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

    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);
	}
EOF
vi)))

    (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(origin.get_steps() + lag);

      }
EOF
  (if (lookup-def 'V_t const-defs) threshold "")
  )))
  
  (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 <cmath>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_exp.h>
"  sysname)))
  (if (not (null? steady-state-index-map))
      (pp indent ("#include <gsl/gsl_multimin.h>")))
  )


(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 
	 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\"
#include <gsl/gsl_odeiv.h>
"))
  
  (pp indent "namespace nest {" ,nl)

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

  (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 
      "enum StateVecElems {"
      (,(slp ", " (map (lambda (x) 
			 (let ((n (string-upcase (->string (first x))))
			       (ni (second x)))
			   (s+ n " = " ni)))
		       state-index-map)))
      "};"
      (,(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))
      (,(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
#<<EOF

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

  (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 (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) (hash-table-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))
	     (synapse-info         (nemo:post-synaptic-conductance-query sys))
	     (i-synapses           (lookup-def 'i-synapses synapse-info))

	     (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) (rate/reaction-power sys n)) sts))
					    (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) 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) (rate/reaction-power sys n)) sts))
					  (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) 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) (rate/reaction-power sys n)) sts))
					  (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) 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) (rate/reaction-power sys n)) sts))
				       (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) 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-eqs  (fold  (lambda (i-synapse ax) 
			      (cons (list (first i-synapse) (third i-synapse) (fourth 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  (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 (cons "(-node.B_.I_stim_)" i-names))) ,mcap)))
			  (list 'v 0.0)))
	     
	     (v-init-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 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)))
	       (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 const-defs state-index-map 
			    rate-eq-defs reaction-eq-defs asgn-eq-defs
			    pool-ions mcap 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 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-init-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 
		  synapse-info
		  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 
		  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 ()
		 (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 
		  synapse-info
		  indent indent+)
		 (pp indent ,nl)
		 ))
	     
	     (close-output-port cpp-output)
	     (close-output-port hpp-output)
	   
	   ))
      ))
  ))


)