;;       
;; 
;; An extension for generating plain equational representation from NEMO models.
;;
;; Copyright 2008-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-eqn

	(nemo:eqn-translator)

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

(define eqn-builtin-consts
  `())

(define eqn-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 (eqn-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/EQN expr)
  (match expr 
	 (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/EQN 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 (eqn-name s)) (map (lambda (x) (rhsexpr/EQN x)) es)) expr))
	 (id          (if (symbol? id) (eqn-name id) id))))

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

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


(define group/EQN   (doc:block 2 (doc:text "(") (doc:text ")")))
(define block/EQN   (doc:block 2 (doc:empty) (doc:empty)))
(define (stmt/EQN x) x)


(define (ifthen/EQN 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/EQN 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/EQN e1)))
			   (doc:group (doc:nest 2 e2))))))
	

(define (format-op/EQN 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/EQN indent op (take lst n/2 )) op1 
				   (format-op/EQN indent op (drop lst n/2 )))
			     (doc:space)))))))))

(define (format-fncall/EQN indent op args)
  (let ((op1 (doc:text (->string op))))
    (doc:cons op1 (group/EQN ((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? ) (eqn-name expr))
	 ((? atom? ) expr)))


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

       ((op . rest)  
       (let ((op (case op ((pow)  '^)  ((ln) 'log) (else op))))
	 (let ((fe
		(if (member op eqn-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/EQN indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/EQN indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if (or mul? plmin?) (group/EQN fx) fx)))) rest)))
			((*)  
			 (format-op/EQN indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/EQN indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if plmin? (group/EQN fx) fx)))) rest)))
			
			((^)  
			 (format-op/EQN indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/EQN indent+ x)))
						   (if (or (symbol? x)  (number? x)) fx
						       (if (or mdiv? plmin?) (group/EQN fx) fx)))) rest)))
			
			(else
			 (format-op/EQN indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/EQN indent+ x))) fx)) rest)))))
		    
		    (let ((op (case op ((neg) '-) (else op))))
		      (format-fncall/EQN indent op (map (lambda (x) (format-expr/EQN indent+ x)) rest))))))
	   (if rv 
	       (stmt/EQN (format-op/EQN indent " = " (list (format-expr/EQN indent+ rv ) fe)))
	       fe))))
      
      (else  (let ((fe (doc:text (->string expr))))
	       (if rv 
		   (stmt/EQN (format-op/EQN indent " = " (list (format-expr/EQN indent+ rv ) fe)))
		   fe)))))))
	       
	  
(define (expr->string/EQN x . rest)
  (let-optionals rest ((rv #f) (width 72))
    (sdoc->string (doc:format width (format-expr/EQN 2 x rv)))))


(define (expeuler dt name rhs)
  (define (isname? x) (equal? x name))
  (let ((res 
	 (match rhs
		((or ('- A ('* B (and x (? isname?))))
		     ('+ ('neg ('* B (and x (? isname?)))) A))
		 (let ((xexp (string->symbol (s+ x 'exp))))
		   `(let ((,xexp (exp (* (neg ,B) ,dt))))
		      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B))))))

		((or ('- A ('* (and x (? isname?)) . B))
		     ('+ ('neg ('* (and x (? isname?)) . B)) A))
		 (let ((xexp (string->symbol (s+ x 'exp)))
		       (B1   (if (null? (cdr B)) (car B) `(* ,@B))))
		   `(let ((,xexp (exp (* (neg ,B1) ,dt))))
		      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B1))))))
		
		(('+ ('neg ('* (and x1 (? isname?)) Alpha))
		     ('* ('- 1 (and x2 (? isname?))) Beta))
		 (let ((A  Alpha)
		       (B  `(+ ,Alpha ,Beta)))
		   (let ((xexp (string->symbol (s+ x1 'exp))))
		     `(let ((,xexp (exp (* (neg ,B) ,dt))))
			(+ (* ,x1 ,xexp) (* (- 1 ,xexp) (/ ,A ,B)))))))
		
		(('let bnds body)
		 `(let ,bnds ,(expeuler dt name body)))
		
		(else (nemo:error 'nemo:expeuler ": unable to rewrite equation " rhs 
				  "in exponential Euler form")))))

    res))



(define (state-init n init)
  (let* ((init  (rhsexpr/EQN init))
	 (init1 (canonicalize-expr/EQN init)))
    (list (eqn-name n) init1)))


(define (asgn-eq n rhs)
  (let* ((fbody   (rhsexpr/EQN rhs))
	 (fbody1  (canonicalize-expr/EQN fbody)))
    (list (eqn-name n) fbody1)))


(define (reaction-eq n open transitions conserve)
  (if (symbol? open)
      (list (eqn-name n) (eqn-state-name n open))
      (list (eqn-name n) (sum (map (lambda (x) (eqn-state-name n x)) open)))
      ))


(define (reaction-transition-eqs n initial open transitions conserve power )
  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve eqn-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 (and cnode (= (first cnode) (first s) )) ax
				(let* ((out   (out-edges (first s)))
				       (in    (in-edges (first s)))
				       (open? (eq? (second s) open))
				       (name  (eqn-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/EQN rhs1))
					 (fbody1 (canonicalize-expr/EQN fbody0)))
				    (cons (list name fbody1) 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/EQN rhs))
						  (dy      (eqn-name name) ))
					      (cons (list dy (canonicalize-expr/EQN 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 (eqn-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 (eqn-name name) transitions conserve
									eqn-state-name) ax) 
						    ax))
				      (else  ax))
			       ax))))
	    ax lst))
   (list) poset))


(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)
  (lambda (globals 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)))
	(pp indent ,nl (function ,(eqn-name n) (,(slp ", " vars)) ))
	(let* ((body0 (rhsexpr/EQN body))
	       (body1 (canonicalize-expr/EQN body0))
	       (lbs   (enum-bnds body1 (list))))
	  (pp indent+ ,(expr->string/EQN body1))
	  (pp indent end))
	  ))
    ))


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

  (pp indent ,nl ("State equations: "))

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

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

  	   
  (for-each (lambda (def)
	      (let ((n (first def)) )
		(pp indent+
		    ,(expr->string/EQN 
		      (second def)
		      (s+ "d(" n ")") 
		      ))))
	    rate-eq-defs)


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

)


(define (output-init sysname globals 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 ("Initial values: "))
    
    (let* ((init-eqs 
	    (append 
	     
	     const-defs
	     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)))))))

      (for-each (lambda (n)
		  (let ((b  (lookup-def n init-eqs)))
		    (if b (pp indent+ ,(expr->string/EQN b (eqn-name n))))))
		init-order))

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

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

(define (nemo:eqn-translator sys filename)
  (define (cid x)  (second x))
  (define (cn x)   (first x))
  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
    (let ((imports  ((dis 'imports)  sys))
	  (exports  ((dis 'exports)  sys)))
      (let* ((indent      0)
	     (indent+     (+ 2 indent ))
	     (sysname     (eqn-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) eqn-builtin-consts))
				       (let ((v1 (canonicalize-expr/EQN (second nv))))
					 (list (eqn-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 ,(eqn-name i) ,(eqn-name e) ,erev)))
				 (lookup-def 'perm-ions gate-complex-info)))
	     (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map eqn-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 eqn-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))
			      (gates             (filter (lambda (x) (equal? (car x) 'gate)) subcomps))
			      (sts               (map (lambda (gate) ((dis 'component-exports) sys (cid gate))) gates))
			      )


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

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


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

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

				  (else
				   (let* ((i     (eqn-name (s+ 'i (cn perm))))
					  (e     (car permqs))
					  (gmax  (car ((dis 'component-exports) sys (cid pore))))
					  (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
					  (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
					  (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
									(match gpwr ((x)  x) (else `(* ,@gpwr))))
								      gpwrs))))
					  )
				     (list i e gion (eqn-name (s+ 'i_ label) ))))))
			       
			       ((and perm pore)
				(case (cn perm)
				  ((non-specific)
				   (let* ((i     (eqn-name 'i))
					  (e     (car permqs))
					  (gmax  (car ((dis 'component-exports) sys (cid pore)))))
				     (list i e gmax (eqn-name (s+ 'i_ label) ))))
				  (else
				   (nemo:error 'nemo:eqn-translator ": invalid ion channel definition " label))))
			       
			       ((and acc pore (pair? gates))
				(let* ((i     (eqn-name (s+ 'i (cn acc))))
				       (gmax  (car ((dis 'component-exports) sys (cid pore))))
				       (pwrs  (map (lambda (st) (map (lambda (n) (rate/reaction-power sys n)) st)) sts))
				       (gpwrs (map (lambda (st pwr) (map (lambda (s p) (if p `(pow ,s ,p) s)) st pwr)) sts pwrs))
				       (gion  `(* ,gmax ,(sum (map (lambda (gpwr) 
								     (match gpwr ((x)  x) (else `(* ,@gpwr))))
								   gpwrs))))
				       )
				  (list i #f gion  (eqn-name (s+ 'i_ label) ))))
			       (else (nemo:error 'nemo:eqn-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 (eqn-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/EQN sum0))
						  (sum2  (canonicalize-expr/EQN 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/EQN (if e `(* ,gion (- v ,e)) gion)))
								   (expr1 (canonicalize-expr/EQN expr0)))
							      (cons (list ii expr1) eqs)))))))
				    
				    ((i e gion ii)
				     (let* ((expr0  (rhsexpr/EQN (if e `(* ,gion (- v ,e)) gion)))
					    (expr1  (canonicalize-expr/EQN 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 eqn-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) (eqn-name (s+ 'i_ (first gate-complex)))) gate-complexes )
							(map (lambda (i-gate) (eqn-name (s+ 'i_ (second i-gate)))) i-gates )

							))))
	     (v-eq    (if (and mcap (member 'v globals))
			  (list 'v (rhsexpr/EQN `(/ (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 1 (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 1 (list)) 
						       rate-eq-defs)))
					(second acc)))
	     
	     (dfenv 
	      (map (lambda (x) (let ((n (first x)))
				 (list n (eqn-name (s+ "d_" n )))))
		   defuns))


	     )
	
	
	
	(for-each 
	 (lambda (a)
	   (let ((acc-ion   (car a)))
	     (if (assoc acc-ion perm-ions)
		 (nemo:error 'nemo:eqn-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:eqn-translator 
			     ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
	 pool-ions)

	(let ((output  (open-output-file filename)))
	  
	  (with-output-to-port output
	    (lambda ()
	      (if (not (null? globals)) 
		  (pp indent ("model global quantities:"  ,(slp " " globals)) ,nl))))

	  (with-output-to-port output
	    (lambda ()
	      (if (not (null? perm-ions)) 
		  (pp indent ("imported quantities:"  ,(slp " " (map car imports)) ,nl)))))

	  (with-output-to-port output
	    (lambda ()
	      (if (not (null? perm-ions)) 
		  (pp indent ("permeating ions:"  ,(slp " " (delete-duplicates (map car perm-ions))) ,nl)))))

	  (with-output-to-port output
	    (lambda ()
	      (if (not (null? acc-ions)) 
		  (pp indent ("accumulating ions:"  ,(slp " " (delete-duplicates (map car acc-ions)))) ,nl))))

	  (with-output-to-port output
	    (lambda ()
	      (if (not (null? pool-ions)) 
		  (pp indent ("pool ions:"  ,(slp " " (map cadddr epools))) ,nl))))

	  (with-output-to-port output
	    (lambda ()
	      (if (not (null? i-names)) 
		  (pp indent ("ionic currents:"  ,(slp " " (map car i-eqs))) ,nl))))

	  ;; initial values function
	  (with-output-to-port output
	    (lambda ()
	      (output-init sysname globals 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 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+)))

	   ;; user-defined functions
	   (let* (;;(with       (inexact->exact (round (/ (abs (- max-v min-v)) step))))
		  (define-fn  (make-define-fn)))
	     (for-each (lambda (fndef) 
			 (if (not (member (car fndef) builtin-fns))
			     (with-output-to-port output
			       (lambda ()
				 (apply define-fn (cons globals (cons indent fndef)))
				 (pp indent ,nl)))
			 ))
		       defuns))
	 
	   (close-output-port output))

	))
    ))

)