;;       
;; 
;; An extension for translating NEMO models to NMODL descriptions.
;;
;; Copyright 2008-2010 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-nmodl

	(nemo:nmodl-translator)

	(import scheme chicken utils data-structures lolevel srfi-1 srfi-13 )

	(require-extension lolevel datatype matchable strictly-pretty 
			   environments varsubst datatype 
			   nemo-core nemo-utils nemo-ionch)
(declare (lambda-lift))

(define (nmodl-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 (nmodl-state-name n s)
  (nmodl-name (if n (s+ n s) s)))

(define (rhsvars rhs)
  (enum-freevars rhs (list) (list)))

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


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


(define (letblk/NMODL 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 e1))
			   (doc:group (doc:nest 2 e2))))))
	
(define ifthen/NMODL  (doc:ifthen 0 (doc:text "if") (doc:text "") (doc:text "else")))
(define group/NMODL   (doc:block 2 (doc:text "(") (doc:text ")")))
(define block/NMODL   (doc:block 2 (doc:text "{") (doc:text "}")))
(define binop/NMODL   (doc:binop 2))

(define (format-op/NMODL indent op args)
  (let ((op1 (doc:text (->string op))))
    (let ((res
	   (if (null? args) op1
	       (match args
		      ((x)           (doc:connect op1 x))
		      ((x y)         (binop/NMODL x op1 y))
		      ((x y z)       (binop/NMODL x op1 (binop/NMODL y op1 z)))
		      (lst           (let* ((n   (length lst))
					    (n/2 (inexact->exact (round (/ n 2)))))
				       (binop/NMODL (format-op/NMODL indent op (take lst n/2 )) op1 
						    (format-op/NMODL indent op (drop lst n/2 )))))))))
      res)))


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

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

(define nmodl-builtin-consts
  `(celsius diam))

(define nmodl-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 (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? ) (nmodl-name expr))
	 ((? atom? ) expr)))


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

       ((op . rest)  
       (let ((op (case op ((pow)  '^) ((abs) 'fabs) ((ln) 'log) (else op))))
	 (let ((fe
		(if (member op nmodl-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/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/NMODL indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if (or mul? plmin?) (group/NMODL fx) fx)))) rest)))
			((*)  
			 (format-op/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/NMODL indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if plmin? (group/NMODL fx) fx)))) rest)))
			
			((^)  
			 (format-op/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/NMODL indent+ x)))
						   (if (or (symbol? x)  (number? x)) fx
						       (if (or mdiv? plmin?) (group/NMODL fx) fx)))) rest)))
			
			(else
			 (format-op/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/NMODL indent+ x))) fx)) rest)))))
		    
		    (let ((op (case op ((neg) '-) (else op))))
		      (format-fncall/NMODL indent op (map (lambda (x) (format-expr/NMODL indent+ x)) rest))))))
	   (if rv 
	       (format-op/NMODL indent " = " (list (format-expr/NMODL indent+ rv ) fe))
	       fe))))
      
      (else  (let ((fe (doc:text (->string expr))))
	       (if rv 
		   (format-op/NMODL indent " = " (list (format-expr/NMODL indent+ rv ) fe))
		   fe)))))))
	       

	  
(define (expr->string/NMODL x . rest)
  (let-optionals rest ((rv #f) (width 72))
    (sdoc->string (doc:format width (format-expr/NMODL 2 x rv)))))
  

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

       ((op . rest)  
	(let ((op (case op ((pow)  '^) ((abs) 'fabs) (else op))))
	  (let ((fe
		(if (member op nmodl-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-lineq-op/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-lineq/NMODL indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if (or mul? plmin?) (group/NMODL fx) fx)))) rest)))
			((*)  
			 (format-lineq-op/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-lineq/NMODL indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if plmin? (group/NMODL fx) fx)))) rest)))
			
			((^)  
			 (format-lineq-op/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-lineq/NMODL indent+ x)))
						   (if (or (symbol? x)  (number? x)) fx
						       (if (or mdiv? plmin?) (group/NMODL fx) fx)))) rest)))
			
			(else
			 (format-lineq-op/NMODL indent op 
					  (map (lambda (x) 
						 (let ((fx (format-lineq/NMODL indent+ x))) fx)) rest)))))
		    
		    (case op
		      ((neg) (format-lineq-op/NMODL indent '* (map (lambda (x) (format-lineq/NMODL indent+ x)) 
							     (cons "(-1)" rest))))
		      (else  (format-fncall/NMODL indent op (map (lambda (x) (format-lineq/NMODL indent+ x)) 
								 rest)))))))

	   (if rv (format-lineq-op/NMODL indent " = " (list (format-lineq/NMODL indent+ rv ) fe)) fe))))
      
      (else  (let ((fe (doc:text (->string expr))))
	       (if rv 
		   (format-lineq-op/NMODL indent " = " (list (format-lineq/NMODL indent+ rv ) fe))
		   fe)))))))
	       

(define (conserve-lineq->string/NMODL x val . rest)
  (let-optionals rest ((width 72))
    (s+ "CONSERVE " (sdoc->string (doc:format width (format-lineq/NMODL 2 x #f))) 
	" = " (number->string val))))
  

(define (make-define-fn table? min-v max-v with depend)
  (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)))
	(pp indent ,nl (FUNCTION ,n (,(sl\ ", " vars)) "{" ))
	(let* ((body1 (canonicalize-expr/NMODL (rhsexpr/NMODL body)))
	       (lbs   (enum-bnds body1 (list))))
	  (if (not (null? lbs)) (pp indent+ (LOCAL ,(sl\ ", " lbs))))
	  (if (and table? min-v max-v with)
	      (match vars
		     (('v)  (pp indent+ (TABLE ,@(if depend `(DEPEND ,depend) `(""))
					       FROM ,min-v TO ,max-v WITH ,with)))
		     (else  (void))))
	  (pp indent+ ,(expr->string/NMODL body1 (nmodl-name n))))
	(pp indent "}")))  
    ))

(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 (reaction-transition-eqs n initial open transitions power method)
  (match-let (((g  node-subs)  (transitions-graph n open transitions nmodl-state-name)))
     (let* ((out-edges  (g 'out-edges))
	    (in-edges   (g 'in-edges))
	    (nodes   ((g 'nodes)))
	    (snode   (find (lambda (s) (not (eq? (second s) open))) nodes)))
       ;; generate differential equations for each state in the transitions system
      (let ((eqs    (fold (lambda (s ax) 
			    (if (= (first snode) (first s) ) ax
				(let* ((out   (out-edges (first s)))
				       (in    (in-edges (first s)))
				       (open? (eq? (second s) open))
				       (name  (nmodl-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/NMODL rhs1)))
				    (case method
				      ((expeuler)  (cons (list name (canonicalize-expr/NMODL (expeuler 'dt name fbody0))) 
							 ax))
				      (else        (cons (list name (canonicalize-expr/NMODL fbody0))
							 ax))
				      )))))
			  (list) nodes)))
	eqs))))
	    
	

(define (reaction-keqs n initial open transitions power)
  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) 
				       nemo:binding? identity nemo:bind nemo:subst-term))
	 (state-list     (let loop ((lst (list)) (tlst transitions))
			   (if (null? tlst)  (delete-duplicates lst eq?)
			       (match (car tlst) 
				      (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
				       (loop (cons* s0 s1 lst) (cdr tlst)))
				      (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
				       (loop (cons* s0 s1 lst) (cdr tlst)))
				      (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
				       (loop (cons* s0 s1 lst) (cdr tlst)))
				      (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
				       (loop (cons* s0 s1 lst) (cdr tlst)))
				      (else
				       (nemo:error 'nemo:nmodl-reaction-keqs ": invalid transition equation " 
						   (car tlst) " in state complex " n))
				      (else (loop lst (cdr tlst)))))))
	 (state-subs     (fold (lambda (s ax) (subst-extend s (nmodl-state-name n s) ax)) subst-empty state-list)))
    ;; generate kinetic equations for each edge in the transitions system
    (list n
	  (map
	   (lambda (e) 
	     (match e
		    (('-> s0 s1 rexpr)
		     (let ((i  (lookup-def s0 state-subs))
			   (j  (lookup-def s1 state-subs)))
		       `(-> ,i ,j ,(canonicalize-expr/NMODL 
				    (subst-convert rexpr state-subs)))))
		    
		    ((s0 '-> s1 rexpr)
		     (let ((i  (lookup-def s0 state-subs))
			   (j  (lookup-def s1 state-subs)))
		       `(-> ,i ,j ,(canonicalize-expr/NMODL 
				    (subst-convert rexpr state-subs)))))
		    
		    (('<-> s0 s1 rexpr1 rexpr2)
		     (let ((i  (lookup-def s0 state-subs))
			   (j  (lookup-def s1 state-subs)))
		       `(<-> ,i ,j 
			     ,(canonicalize-expr/NMODL (subst-convert rexpr1 state-subs))
			     ,(canonicalize-expr/NMODL (subst-convert rexpr2 state-subs)))))
		    
		    ((s0 '<-> s1 rexpr1 rexpr2)
		     (let ((i  (lookup-def s0 state-subs))
			   (j  (lookup-def s1 state-subs)))
		       `(<-> ,i ,j 
			     ,(canonicalize-expr/NMODL (subst-convert rexpr1 state-subs))
			     ,(canonicalize-expr/NMODL (subst-convert rexpr2 state-subs)))))
		    
		    
		    (else (nemo:error 'nemo:nmodl-reaction-keqs ": invalid transition equation " 
				      e " in state complex " n))))
	   transitions))))
	


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


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


(define (reaction-eq n open transitions)
  (list (nmodl-name n) (nmodl-name (nmodl-state-name n open))))

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

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


(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 kinetic method)
  (fold-right
   (lambda (lst ax)
     (fold  (lambda (x ax) 
	      (match-let (((i . n)  x))
			 (let ((en (environment-ref sys n)))
			   (if (and (not (member n kinetic)) (nemo:quantity? en))
			       (cases nemo:quantity en
				      (REACTION  (name initial open transitions conserve power) 
						 (append (reaction-transition-eqs name initial open transitions 
										  power method) ax))
				      
				      (RATE (name initial rhs)
					    (let ((fbody0  (rhsexpr/NMODL rhs))
						  (dy      name ))
					      (case method
						((expeuler)  
						 (cons (list dy (canonicalize-expr/NMODL (expeuler 'dt name fbody0))) 
						       ax))
						(else
						 (cons (list dy (canonicalize-expr/NMODL fbody0)) ax)))))

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


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


(define (poset->state-init-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 (nemo:rhs? initial)
						     (cons* (state-init name initial) 
							    (state-init (nmodl-state-name name open) name) ax) 
						     ax))
				      (RATE  (name initial rhs)
					     (if (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:lineq? conserve))
						    (cons (state-lineqs (nmodl-name name) transitions conserve
									nmodl-state-name) ax) 
						    ax))
				      (else  ax))
			       ax))))
	    ax lst))
   (list) poset))


(define (find-locals defs)
  (concatenate (map (lambda (def) (match def (('let bnds _) (map first bnds)) (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 (nemo:nmodl-translator sys . rest)
  (define (cid x)  (second x))
  (define (cn x)   (first x))
  (let-optionals rest ((method 'cnexp) (table? #f) (min-v -100) (max-v 100) (step 0.5) 
		       (depend #f)  (kinetic (list)) (linear? #f))
  (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 ))
	     (table-with    (and table? (inexact->exact (round (/ (abs (- max-v min-v)) step)))))
	     (eval-const    (dis 'eval-const))
	     (sysname       (nmodl-name ((dis 'sysname) sys)))
	     (consts        ((dis 'consts)  sys))
	     (asgns         ((dis 'asgns)   sys))
	     (states        ((dis 'states)  sys))
	     (kinetic       (cond ((eq? kinetic 'all) (delete-duplicates (map first states)))
				  ((symbol? kinetic)  (list kinetic))
				  (else kinetic)))
	     (reactions     ((dis 'reactions) sys))
	     (rates         ((dis 'rates) 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)))
	     (ionch-info    (nemo:ionch-query sys))
	     (ionchs        (lookup-def 'ion-channels ionch-info))
	     (perm-ions     (map (match-lambda ((comp i e erev) `(,comp ,(nmodl-name i) ,(nmodl-name e) ,erev)))
				 (lookup-def 'perm-ions ionch-info)))
	     (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map nmodl-name (list i in out)))))
				 (lookup-def 'acc-ions ionch-info)))
	     (epools        (lookup-def 'pool-ions ionch-info))
	     (pool-ions     (map (lambda (lst) (map nmodl-name lst)) epools))

	     (i-gates       (lookup-def 'i-gates ionch-info))

	     (has-kinetic?  (or (not (null? (filter (lambda (x) (member (car x) kinetic)) states)))))
	     (has-ode?      (or (not (null? (filter (lambda (x) (not (member (car x) kinetic))) states)))
				(not (null? pool-ions))))

	     (asgn-eq-defs       (poset->asgn-eq-defs poset sys))
	     (reaction-eq-defs   (poset->reaction-eq-defs poset sys kinetic)) 
	     (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys kinetic method)))
	     (kstate-eq-defs     (poset->kinetic-eq-defs poset sys kinetic))
	     (conserve-eq-defs   (poset->state-conserve-eq-defs poset sys))
	     (state-init-defs    (poset->state-init-defs poset sys))

	     )
	
	   (pp indent ,nl (TITLE ,sysname))
	   
	   (pp indent ,nl (NEURON "{"))
	   (if (not (null? exports)) (pp indent+ (RANGE ,(sl\ ", " (map nmodl-name exports)))))

	   (for-each (lambda (x)
		       (case (first x)
			 ((non-specific) 
			  (pp indent+ (RANGE ,(third x))
			      (NONSPECIFIC_CURRENT ,(second x))))
			 (else 
			  (pp indent+ (RANGE ,(second x))
			      (USEION ,(first x) READ ,(third x) WRITE ,(second x))))))
		     (delete-duplicates perm-ions (lambda (x y) (eq? (car x) (car y)))))


	   (for-each (lambda (acc-ion)
		       (let ((pool-ion (assoc (first acc-ion) pool-ions)))
			 (if pool-ion
			     (pp indent+ (RANGE ,(second acc-ion))
				 (USEION ,(first acc-ion) 
					 READ  ,(sl\ ", " (list (third acc-ion) (fourth acc-ion) (second pool-ion)))
					 WRITE ,(sl\ ", " (list (second acc-ion) (third pool-ion )))))
			     (pp indent+ (RANGE ,(second acc-ion))
				 (USEION ,(first acc-ion) 
					 READ ,(sl\ ", "  (list (third acc-ion) (fourth acc-ion) ))
					 WRITE ,(second acc-ion))))))
		       (delete-duplicates acc-ions (lambda (x y) (eq? (car x) (car y)))))

	   (let* ((const-names   (map first consts))
		  (is-const?     (lambda (x) (member x const-names)))
		  (range-consts  (delete-duplicates 
				  (fold (lambda (def ax) 
					  (let* ((rhs   (second def))
						 (vars  (rhsvars rhs)))
					    (append (filter is-const? vars) ax)))
					(list) asgn-eq-defs ))))
	     (if (not (null? range-consts)) (pp indent+ (RANGE ,(sl\ ", " range-consts)))))
	   
	   
	   (pp indent "}")
	   
	   (let* ((define-fn (make-define-fn table? min-v max-v table-with depend)))
	     (for-each (lambda (fndef) 
			 (if (not (member (car fndef) builtin-fns))
			     (apply define-fn (cons indent fndef)))) 
		       defuns))
	   
	   (let* ((parameter-defs 
		   (filter-map
		    (lambda (nv)
		      (and (not (member (first nv) nmodl-builtin-consts))
			   (let ((v1 (canonicalize-expr/NMODL (second nv))))
			     (list (first nv) v1))))
		    consts))
		  (parameter-locals  (find-locals (map second parameter-defs)))
		  (state-defs 
		   (append
		    (map (lambda (st)
			   (if (pair? st) (nmodl-state-name (first st) (second st)) 
			       (nmodl-name st)))
			 states)
		    (map nmodl-name reactions)))
		  (assigned-defs
		   (filter-map
		    (lambda (x) 
		      (let ((x1 (nmodl-name x)))
			(and (not (or (member x1 state-defs) (assoc x1 parameter-defs)))
			     x1)))
		    (delete-duplicates
		     (append asgns
			     (map first imports) 
			     (map second perm-ions) (map third perm-ions)
			     (map second acc-ions)  (map fourth acc-ions)
			     (map second pool-ions) (map third pool-ions)
			     )))))
		  
	     (pp indent ,nl (PARAMETER "{"))
	     (if (not (null? parameter-locals)) (pp indent+ (LOCAL ,(sl\ ", " parameter-locals))))
	     (for-each (lambda (def)
			 (let ((n (nmodl-name (first def))) (b (second def)))
			   (pp indent+ ,(expr->string/NMODL b n)))) parameter-defs)
	     (case method  ((expeuler)  (pp indent+ dt)))
	     (pp indent "}")
	     (pp indent ,nl (STATE "{"))
	     (for-each (lambda (x) (pp indent+ ,x)) state-defs)
	     (pp indent "}")
	     (pp indent ,nl (ASSIGNED "{"))
	     (for-each (lambda (x) (pp indent+ ,x)) assigned-defs)
	     (pp indent "}"))

	   (if (not (null? asgns))
	       (begin
		 (pp indent ,nl (PROCEDURE asgns () "{"))
		 (let ((locals    (find-locals (map second asgn-eq-defs))) )
		   (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals)))))
#|
                 This seems to cause a segmentation fault in nrnoc:

		 (if (and table? min-v max-v table-with)
		     (pp indent+ (TABLE ,(sl\ ", " (map first asgn-eq-defs))
				  ,@(if depend `(DEPEND ,depend) `(""))
				  FROM ,min-v TO ,max-v WITH ,table-with)))
|#

		 (for-each (lambda (def)
			     (let ((n (nmodl-name (first def)) )
				   (b (second def)))
			       (pp indent+ ,(expr->string/NMODL b n)))) asgn-eq-defs)
		 (pp indent "}")))
	     
	   (if (not (null? reactions))
	       (begin
		 (pp indent ,nl (PROCEDURE reactions () "{"))
		 (let ((locals    (find-locals (map second reaction-eq-defs))) )
		   (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
		   (for-each (lambda (def)
			       (let ((n (nmodl-name (first def))) (b (second def)))
				 (pp indent+ ,(expr->string/NMODL b n)))) reaction-eq-defs))
		 
		 (pp indent "}")))

	   (if (not (null? pool-ions))
	       (begin
		 (pp indent ,nl (PROCEDURE pools () "{"))
		 (for-each (lambda (pool-ion)
			     (pp indent+ (,(third pool-ion) = ,(first pool-ion))))
			   pool-ions)
		 (pp indent "}")))

	   (pp indent ,nl (BREAKPOINT "{"))
	   (let* ((i-eqs (filter-map
			  (lambda (ionch) 

			    (let* ((label     (first ionch))
				   (n         (second ionch))
				   (subcomps  ((dis 'component-subcomps) sys n))
				   (acc       (lookup-def 'accumulating-substance subcomps))
				   (perm      (lookup-def 'permeating-substance subcomps))
				   (permqs    (and perm ((dis 'component-exports) sys (cid perm))))
				   (pore      (lookup-def 'pore subcomps))
				   (gate      (lookup-def 'gate subcomps))
				   (sts       (and gate ((dis 'component-exports) sys (cid gate)))))

			      (if (not pore)
				  (nemo:error 'nemo:nmodl-translator ": ion channel definition " label
					      "lacks any pore components"))

			      (cond ((and perm pore gate)
				     (case (cn perm)
				       ((non-specific)
					(let* ((i     (nmodl-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)))
				       (else
					(let* ((i     (nmodl-name (s+ 'i (cn perm))))
					       (e     (nmodl-name (s+ 'e (cn perm))))
					       (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)))))

				    ((and perm pore)
				     (case (cn perm)
				       ((non-specific)
					(let* ((i     (nmodl-name 'i))
					       (e     (car permqs))
					       (gmax  (car ((dis 'component-exports) sys (cid pore)))))
					  (list i e gmax)))
				       (else
					(nemo:error 'nemo:nmodl-translator ": ion channel definition " label
						    (s+ "(" n ")")
						    "lacks gate component"))))

				      ((and acc pore gate)
				       (let* ((i     (nmodl-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)))

				      (else (nemo:error 'nemo:nmodl-translator ": invalid ion channel definition " 
							label))
				      )))
			  ionchs))
		  (i-eqs  (fold  (lambda (i-gate ax) 
				   (let ((i-gate-var (first i-gate)))
				     (cons (list (nmodl-name 'i) #f i-gate-var) 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) . rst))
					  (let loop ((ps ps) (summands (list)))
					    (if (null? ps)
						(let* ((sum0  (sum summands))
						       (sum1  (rhsexpr/NMODL sum0))
						       (sum2  (canonicalize-expr/NMODL sum1)))
						  (cons (list i sum2) ax))
						(match-let (((i e gion) (car ps)))
						  (loop (cdr ps) (cons (if e `(* ,gion (- v ,e)) gion) summands))))))
					   
					 ((i e gion)
					  (let* ((expr0  (rhsexpr/NMODL (if e `(* ,gion (- v ,e)) gion)))
						 (expr1  (canonicalize-expr/NMODL expr0)))
					    (cons (list i expr1) ax)))
					 
					 (else ax)))
				(list) i-bkts))
		  (locals (find-locals (map second i-eqs))))
	     (if (not (null? locals))    (pp indent+ (LOCAL ,(sl\ ", " locals))))
	     (if (not (null? asgns))     (pp indent+ (asgns ())))
	     (if has-ode?
		 (case method
		   ((#f expeuler)  (pp indent+ (SOLVE states)))
		   (else           (pp indent+ (SOLVE states METHOD ,method)))))
	     (if has-kinetic?   (pp indent+  (SOLVE kstates METHOD sparse)))
	     (if (not (null? reactions))   (pp indent+ (reactions ())))
	     (if (not (null? pool-ions)) (pp indent+ (pools ())))
	     (for-each (lambda (p) (pp indent+ ,(expr->string/NMODL (second p) (first p)))) i-eqs)
	     (pp indent "}"))
	   
	   (if has-ode?
	       (let ((locals   (find-locals (map second rate-eq-defs))))
		 (case method
		   ((expeuler) (pp indent ,nl (PROCEDURE states () "{")))
		   (else       (pp indent ,nl (DERIVATIVE states "{"))))
		 (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
		 (let ((prime (case method 
				((expeuler) identity)
				(else  (lambda (x) (s+ x "'"))))))
		   (for-each (lambda (def)
			       (let ((n (prime (first def)))
				     (b (second def)))
				 (pp indent+ ,(expr->string/NMODL b n))))
			     rate-eq-defs))
		 (pp indent "}")))

	   (if has-kinetic?
	       (begin
		 (pp indent ,nl (KINETIC kstates "{"))
		 (let ((locals            (concatenate (find-locals (map third (map second kstate-eq-defs))))))
		   (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
		   (for-each
		    (lambda (def)
		      (let* ((n             (first def))
			     (eqs           (second def))
			     (conserve-eqs  (lookup-def n conserve-eq-defs)))
			(for-each
			 (lambda (eq)
			   (match eq
				  (('-> s0 s1 rexpr) 
				   (pp indent+ (~ ,s0 -> ,s1 (,(expr->string/NMODL rexpr)))))
				  (('<-> s0 s1 rexpr1 rexpr2)  
				   (pp indent+ (~ ,s0 <-> ,s1 (,(expr->string/NMODL rexpr1) #\,
							       ,(expr->string/NMODL rexpr2) 
							       ))))
				  ))
			   eqs)
			(if conserve-eqs
			    (for-each (lambda (eq) 
					(let ((val  (first eq))
					      (expr (third eq)))
					  (pp indent+ ,(conserve-lineq->string/NMODL expr val))))
				    conserve-eqs))
			))
		    kstate-eq-defs))
		 (pp indent "}")))
	   
	   
	   (let ((locals (concatenate (find-locals (map second state-init-defs)))) )
	       (pp indent ,nl (INITIAL "{"))
	       (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
	       (if (not (null? asgns))  (pp indent+ (asgns ())))
	       (for-each (lambda (def)
			   (let ((n (first def)) (b (second def)))
			     (pp indent+ ,(expr->string/NMODL b n)))) state-init-defs)

	       (if has-kinetic?
		   (pp indent+ (SOLVE kstates STEADYSTATE sparse)))
	       
	       (pp indent "}")

	       (pp indent ,nl (PROCEDURE print_state () "{"))

	       (let ((lst (sort (map (compose ->string first) rate-eq-defs) string<?)))
		 (for-each (lambda (x)
			     (pp indent+ (printf (,(s+ #\" x " = %g\\n"  #\") ", " ,x ))))
			   lst))

	       (pp indent "}")

	       ))))))
)