;;       
;; 
;; An extension for translating NEMO models to Matlab/Octave code.
;;
;; Copyright 2008-2011 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-matlab

	(nemo:matlab-translator
	 nemo:octave-translator)

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

(define matlab-builtin-consts
  `())

(define matlab-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 (matlab-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/MATLAB expr)
  (match expr 
	 (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/MATLAB x)) es)))
	 (('pow x y)  (if (and (integer? y)  (positive? y))
			  (if (> y 1)  (let ((tmp (gensym "x")))
					 `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
			      x)
			    expr))
	 ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (matlab-name s)) (map (lambda (x) (rhsexpr/MATLAB x)) es)) expr))
	 (id          (if (symbol? id) (matlab-name id) id))))

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

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


(define group/MATLAB   (doc:block 2 (doc:text "(") (doc:text ")")))
(define block/MATLAB   (doc:block 2 (doc:empty) (doc:empty)))
(define (stmt/MATLAB x) 
  (match x
	 (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
	 (else  (doc:cons x (doc:text ";")))))


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

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

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


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

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


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


(define (reaction-eq n open transitions conserve)
  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve matlab-state-name)))
    (let ((nodes ((g 'nodes))))
;      (if (< 2 (length nodes))
;	  (list (matlab-name n) `(abs (/ ,(matlab-state-name n open) 
;					 ,(sum (filter-map (lambda (x) (and (not (= (first x) (first cnode)))
;									    (matlab-state-name n (second x)))) nodes)))))
	  (list (matlab-name n) (matlab-state-name n open)))))


(define (reaction-transition-eqs n initial open transitions conserve power method)
  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve matlab-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  (matlab-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/MATLAB rhs1)))
				    (case method
				      ((expeuler)  (cons (list name (canonicalize-expr/MATLAB 
								     (expeuler 'dt name fbody0))) 
							 ax))
				      (else        (cons (list name (canonicalize-expr/MATLAB fbody0))
							 ax))
				      )))))
			  (list) nodes)))
	eqs))))
	    
	

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

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

				      (REACTION  (name initial open transitions conserve power) 
						 (append (reaction-transition-eqs name initial open transitions 
										  conserve power method) ax))
				      
				      (RATE (name initial rhs)
					    (let ((fbody0  (rhsexpr/MATLAB rhs))
						  (dy      (matlab-name name) ))
					      (case method
						((expeuler)  
						 (cons (list dy (canonicalize-expr/MATLAB (expeuler 'dt name fbody0))) 
						       ax))
						(else
						 (cons (list dy (canonicalize-expr/MATLAB fbody0)) ax)))))

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


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


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

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

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


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


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


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

  
(define (make-define-fn)
  (lambda (indent n proc)
    (let ((lst (procedure-data proc))
	  (indent+ (+ 2 indent)))
      (let ((retval   (matlab-name (gensym "retval")))
	    (rt       (lookup-def 'rt lst))
	    (formals  (lookup-def 'formals lst))
	    (vars     (lookup-def 'vars lst))
	    (body     (lookup-def 'body lst)))
	(pp indent ,nl (function ,retval = ,(matlab-name n) (,(slp ", " vars)) ))
	(let* ((body0 (rhsexpr/MATLAB body))
	       (body1 (canonicalize-expr/MATLAB body0))
	       (lbs   (enum-bnds body1 (list))))
	  (pp indent+ ,(expr->string/MATLAB body1 retval))
	  (pp indent end))
	  ))))


(define (output-dy sysname method 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 (function dy = ,sysname (,(slp ", " (case method
						       ((lsode)  `(y t))
						       ((odepkg) `(t y))
						       (else     `(y t)))))))
  
  (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
	   
  (if (member 'v globals)
      (let ((vi (lookup-def 'v state-index-map)))
	(if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
	   
  (for-each (lambda (def)
	      (let ((n (first def)) )
		(pp indent+
		    ,(s+ "% rate eq " n)
		    ,(expr->string/MATLAB 
		      (s+ "y(" (lookup-def n state-index-map) ")") n))))
	    rate-eq-defs)

  (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+ 
				  ,(s+ "% equation for " n)
				  ,(expr->string/MATLAB b (matlab-name n))))))
		    eq-order))

  
  
  (pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'dy))
  
  (for-each (lambda (def)
	      (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")")) 
		    (b (second def)))
		(pp indent+ ,(s+ "% state " (first def))
		    ,(expr->string/MATLAB b n))))
	    rate-eq-defs)
  
  (for-each (lambda (def) 
	      (pp indent+ ,(expr->string/MATLAB (second def) (first def)))) 
	    i-eqs)

  (if v-eq
      (let ((vi (lookup-def 'v state-index-map)))
	(if vi
	    (pp indent+ ,(expr->string/MATLAB (second v-eq) (s+ "dy(" vi ")"))))))
  
  (pp indent ,nl (end))
)

(define (output-steadystate  sysname method globals steady-state-index-map pool-ions
			     const-defs init-eq-defs asgn-eq-defs rate-eq-defs reaction-eq-defs 
			     indent indent+)
  
  (if (not (null? steady-state-index-map))
      (begin
	(pp indent ,nl (function y = ,(s+ sysname "_steadystate") (x)))
	
	(if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
	
	(for-each 
	 (lambda (def)
	   (let* ((n   (first def)) 
		  (ni  (lookup-def n steady-state-index-map)))
	     (if ni (pp indent+ ,(expr->string/MATLAB (s+ "x(" ni ")") n)))
	     ))
	 rate-eq-defs)
	
	(pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'y))


	(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/MATLAB b (matlab-name n))))))
		    init-order))
	
	(for-each 
	 (lambda (def)
	   (let* ((n   (first def)) 
		  (ni  (lookup-def n steady-state-index-map))
		  (b   (second def)))
	     (if ni (pp indent+ ,(s+ "% state " n)
			,(expr->string/MATLAB b (s+ "y(" ni ")"))))
	     ))
	 rate-eq-defs)
	
	(pp indent ,nl (end))))
)


(define (output-print-state sysname state-index-map indent indent+)
  (pp indent ,nl (function  ,(s+ sysname "_print_state") (y)))
  
  (let ((lst (sort (map (lambda (x) (cons (->string (car x)) (cdr x))) state-index-map)
		   (lambda (x y) (string<? (car x) (car y))))))
    (for-each (lambda (p)
		(let ((n (first p)) (i (second p)))
		  (pp indent+ (,n " = " "y(" ,i ")"))))
	      lst))
  
  (pp indent ,nl (end))
)


(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 (function y0 = ,(s+ sysname "_init") (Vinit)))
    
    (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
    
    (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
    
    (if (member 'v globals)
	(let ((vi (lookup-def 'v state-index-map)))
	  (pp indent+ ,(expr->string/MATLAB `Vinit 'v))
	  (if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))


    (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/MATLAB b (matlab-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/MATLAB n  (s+ "y0(" ni ")"))))))
	      init-eq-defs)

    (if (not (null? steady-state-index-map))
	(begin
	  (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'xs)
	      ,(expr->string/MATLAB 
		`(fsolve ,(s+ #\' sysname "_steadystate" #\' ) xs) 
		"[ys,fval,info]"))
	  
	  
	  (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/MATLAB (s+ "ys(" si ")") n))
			(pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") (s+ "y0(" ni ")")))))))
	   rate-eq-defs)


	  ))

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

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

    )

(define (matlab-translator1 sys . rest)
  (define (cid x)  (second x))
  (define (cn x)   (first x))
  (let-optionals rest ((mode 'multiple) (method 'lsode) (filename #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 ))
	     (eval-const  (dis 'eval-const))
	     (sysname     (matlab-name ((dis 'sysname) sys)))
	     (prefix      sysname)
	     (filename    (or filename (s+ sysname ".m")))
	     (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) matlab-builtin-consts))
				       (let ((v1 (canonicalize-expr/MATLAB (second nv))))
					 (list (matlab-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 ,(matlab-name i) ,(matlab-name e) ,erev)))
				 (lookup-def 'perm-ions gate-complex-info)))
	     (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map matlab-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 matlab-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:matlab-translator ": ion channel definition " label
					 "lacks any pore or permeability components"))


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

			       ((and perm pore gate)
				(case (cn perm)
				  ((non-specific)
				   (let* ((i     (matlab-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  (matlab-name (s+ 'i_ label) ))))

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

							))))
	     (v-eq    (if (and mcap (member 'v globals))
			  (list 'v (rhsexpr/MATLAB `(/ (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 (matlab-name (s+ "d_" n )))))
		   defuns))

	     (output-globals 
	      (lambda () (if (not (null? globals)) (pp indent (global ,(slp " " globals))))))
	     )
	
	
	
	(for-each 
	 (lambda (a)
	   (let ((acc-ion   (car a)))
	     (if (assoc acc-ion perm-ions)
		 (nemo:error 'nemo:matlab-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:matlab-translator 
				": ion species " pool-ion " cannot be declared as both pool and permeating"))))
	    pool-ions)

	   (let ((output  (case mode
			    ((single)  (open-output-file filename))
			    (else #f))))

	     (if output (with-output-to-port output (lambda () (output-globals))))
	     
	     ;; derivative function
	     (let ((output1 (or output (open-output-file (s+ prefix "_dy.m")))))
	       (with-output-to-port output1
		 (lambda ()
		   (output-dy sysname method globals state-index-map 
			      rate-eq-defs reaction-eq-defs asgn-eq-defs
			      pool-ions mcap i-eqs v-eq
			      indent indent+)))
	       (if (not output) (close-output-port output1)))
		 

	     ;; steady state solver
	     (let ((output1 (or output (open-output-file (s+ prefix "_steadystate.m")))))
	       (with-output-to-port output1
		 (lambda ()
		   (output-steadystate sysname method globals steady-state-index-map pool-ions
				       const-defs init-eq-defs asgn-eq-defs rate-eq-defs reaction-eq-defs 
				       indent indent+)))
	       (if (not output) (close-output-port output1)))

	     ;; initial values function
	     (let ((output1 (or output (open-output-file (s+ prefix "_init.m")))))
	       (with-output-to-port output1
		 (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)))
	       (if (not output) (close-output-port output1)))

	   ;; 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))
			     (let ((output1 (or output (open-output-file (s+ (matlab-name (car fndef)) ".m")))))
			       (with-output-to-port output1
				 (lambda ()
				   (apply define-fn (cons indent fndef))
				   (pp indent ,nl)))
			       (if (not output) (close-output-port output1)))))
		       defuns))
	 
	   (if output (close-output-port output)))

	   )))))


(define (nemo:matlab-translator sys . rest)
  (apply matlab-translator1 (cons sys (cons 'multiple (cons 'odepkg rest)))))

(define (nemo:octave-translator sys . rest)
  (apply matlab-translator1 (cons sys (cons 'single rest))))

#|
  
(define (make-define-dfn fenv)

  (define (D vars body)
    (let loop ((vars vars) (exprs (list)))
      (if (null? vars) (if (pair? (cdr exprs)) 
			   `(+ . ,(reverse exprs))  (car exprs))
	  (let ((de (let ((de (differentiate fenv (car vars) body)))
		       (let loop ((e (simplify de)) (de de))
			 (if (equal? de e) de
			     (loop (simplify e) e))))))
	    (loop (cdr vars) (cons de exprs))))))

  (lambda (indent n proc)
    (let ((lst (procedure-data proc))
	  (indent+ (+ 2 indent)))
      (let ((retval   (matlab-name (gensym "retval")))
	    (rt       (lookup-def 'rt lst))
	    (formals  (lookup-def 'formals lst))
	    (vars     (lookup-def 'vars lst))
	    (body     (lookup-def 'body lst)))
	(pp indent ,nl (function ,retval = ,(s+ "d_" (matlab-name n) ) (,(slp ", " vars)) ))
	(let* ((dbody (D vars body))
	       (dbody (canonicalize-expr/MATLAB (rhsexpr/MATLAB dbody))))
	  (pp indent+ ,(expr->string/MATLAB dbody retval))
	  (pp indent end))
	  )))
)

(define (output-jac sysname method globals fenv state-index-map 
		    rate-eq-defs reaction-eq-defs asgn-eq-defs 
		    mcap i-eqs v-eq
		    defuns indent indent+)


  (define (D var expr)
    (let* ((dexpr  (distribute expr))
	   (de (differentiate fenv var dexpr))
	   (se (let loop ((e (simplify de)) (de de))
		 (if (equal? de e) de
		     (loop (simplify e) e)))))
      se))
    
  (pp indent ,nl (function jac = ,(s+ sysname "_jac") 
			   (,(slp ", " `(t y)))))

  (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
  
  (if (member 'v globals)
      (let ((vi (lookup-def 'v state-index-map)))
	(if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
  
  (for-each (lambda (def)
	      (let ((name (first def)))
		(pp indent+ ,(expr->string/MATLAB 
			      (s+ "y(" (lookup-def name state-index-map) ")") name))))
	    rate-eq-defs)

  (for-each (lambda (def)
	      (let ((n (first def)) (b (second def)))
		(pp indent+ 
		    ,(expr->string/MATLAB b n)))) 
	    reaction-eq-defs)
  
  (pp indent+ ,(expr->string/MATLAB `(zeros (length y) (length y)) 'jac))


  (for-each
   (lambda (def)
     
     (let* ((name  (first def))
	    (rhs   (second def))
	    (fv    (enum-freevars rhs '() '() ))
	    (drvs  (map (lambda (st) 
			  (let ((var (first st))
				(expr `(let ,(or (and v-eq (list v-eq)) '())
					 (let ,rate-eq-defs
					   (let ,reaction-eq-defs
					     (let ,(filter-map (lambda (x)
								 (or (assoc x i-eqs) 
								       (assoc x asgn-eq-defs)))
							       fv))
					     ,rhs)))))
			    (D var expr)))
			state-index-map))
	    (vs (map (lambda (st)
		       (let ((ni (lookup-def name state-index-map))
			     (i (second st)))
			 (s+ "jac(" (slp "," (list ni i)) ")")))
		     state-index-map))
	    (ndrvs  (map (lambda (x) (canonicalize-expr/MATLAB (rhsexpr/MATLAB x))) drvs)))
       
       (pp indent+ ,(s+ "% state " name))
       (for-each (lambda (ndrv v)
		   (pp indent+ ,(expr->string/MATLAB ndrv v)))
		 ndrvs vs)
       ))
	    
     (if (and mcap (member 'v globals))
	 (let ((vi (lookup-def 'v state-index-map)))
	   (if vi (cons
		   (list 'v `(/ (neg ,(sum (map first i-eqs))) ,mcap))
		   rate-eq-defs)
	       rate-eq-defs))
	 rate-eq-defs))
  
  (pp indent ,nl (end))
  
)

|#

)