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

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


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

(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 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 (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)
	 (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)
  (list (matlab-name n) (matlab-name (matlab-state-name n open))))


(define (reaction-transition-eqs n initial open transitions power method)
  (match-let (((g  node-subs)  (transitions-graph n open transitions matlab-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  (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 
										  power method) ax))
				      
				      (RATE (name initial rhs)
					    (let ((fbody0  (rhsexpr/MATLAB rhs))
						  (dy      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) 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:lineq? conserve))
						    (cons (state-lineqs (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 table? min-v max-v with)
  (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) (,(sl\ ", " vars)) ))
	(let* ((body1 (canonicalize-expr/MATLAB (rhsexpr/MATLAB body)))
	       (lbs   (enum-bnds body1 (list))))
	  (pp indent+ ,(expr->string/MATLAB body1 retval))
	  (pp indent endfunction))
	  ))))


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

  (pp indent ,nl (function dy = ,sysname (,(sl\ ", " (case method
						       ((lsode)  `(y t))
						       ((odepkg) `(t y))
						       (else     `(y t)))))))
  
  (if (not (null? globals)) (pp indent+ (global ,(sl\ " " 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+ ,(expr->string/MATLAB 
			      (s+ "y(" (lookup-def n state-index-map) ")") n))))
	    rate-eq-defs)
  
  (for-each (lambda (def)
	      (let ((n (matlab-name (first def)) )
		    (b (second def)))
		(pp indent+ ,(expr->string/MATLAB b n)))) 
	    asgn-eq-defs)
  
  (for-each (lambda (def)
	      (let ((n (first def)) (b (second def)))
		(pp indent+ ,(expr->string/MATLAB b n)))) 
	    reaction-eq-defs)
  
  (for-each (lambda (pool-ion)
	      (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
	    pool-ions)
  
  (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 (and mcap (member 'v globals))
      (let ((vi (lookup-def 'v state-index-map)))
	(if vi (pp indent+  ,(expr->string/MATLAB `(/ (neg ,(sum (map first i-eqs))) ,mcap)  
						  (s+ "dy(" vi ")"))))))
  
  (pp indent ,nl (endfunction))
)

  
(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) ) (,(sl\ ", " vars)) ))
	(let* ((dbody (D vars body))
	       (dbody (canonicalize-expr/MATLAB (rhsexpr/MATLAB dbody))))
	  (pp indent+ ,(expr->string/MATLAB dbody retval))
	  (pp indent endfunction))
	  )))
)

(define (output-ddy sysname method globals fenv state-index-map 
		    rate-eq-defs reaction-eq-defs asgn-eq-defs 
		    defuns indent indent+)


  (define (D var expr)
    (let ((de (differentiate fenv var expr)))
      (let loop ((e (simplify de)) (de de))
	(if (equal? de e) de
	    (loop (simplify e) e)))))
    
  (pp indent ,nl (function ddy = ,(s+ sysname "_ddy") 
			   (,(sl\ ", " `(t y)))))
  
  (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)
  
  (pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'ddy))

  (for-each (lambda (def)
	      (let* ((name  (first def))
		     (v     (s+ "ddy(" (lookup-def name state-index-map) ")")) 
		     (b     (second def))
		     (fv    (enum-freevars b '() '() ))
		     (db    (D name 
			       `(let ,(filter-map (lambda (x) (assoc x asgn-eq-defs))
						  fv)
				  ,b)))
		     (cdb    (canonicalize-expr/MATLAB (rhsexpr/MATLAB db))))
		(pp indent+ ,(s+ "# state " name)
		    ,(expr->string/MATLAB cdb v))))
	    rate-eq-defs)
  
  (pp indent ,nl (endfunction))
)


(define (output-steadystate  sysname method globals steady-state-index-map 
			     rate-eq-defs reaction-eq-defs asgn-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 ,(sl\ " " 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))
	
	(for-each (lambda (def)
		    (let ((n (matlab-name (first def)) )
			  (b (second def)))
		      (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs)
	
	(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 (endfunction))))
)


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


(define (output-init sysname globals state-index-map steady-state-index-map 
		     const-defs init-eq-defs rate-eq-defs i-eqs pool-ions
		     pool-ion-i  indent indent+)
  
  (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
  
  (if (not (null? globals)) (pp indent+ (global ,(sl\ " " 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 ")") )))))
  
  (for-each (lambda (def)
	      (let ((n (first def)) (b (second def)))
		(pp indent+ ,(expr->string/MATLAB b n)))) 
	    const-defs)
  
  (for-each (lambda (def)
	      (let* ((n  (first def)) 
		     (ni (lookup-def n state-index-map))
		     (b (second def)))
		(pp indent+ ,(expr->string/MATLAB b n))
		(if ni (pp indent+ ,(expr->string/MATLAB n  (s+ "y0(" ni ")"))))))
	    init-eq-defs)
  
  (for-each (lambda (pool-ion)
	      (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
	    pool-ions)
  
  (if (not (null? steady-state-index-map))
      (begin
	(for-each (lambda (def) 
		    (if (member (first def) pool-ion-i)
			(pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
		  i-eqs)
	
	(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 (pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") (s+ "y0(" ni ")"))))))
	 rate-eq-defs)
	))
  
  
  (for-each (lambda (def) 
	      (if (not (member (first def) pool-ion-i))
		  (pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
	    i-eqs)
  
  (pp indent ,nl (endfunction))
)

(define (nemo:matlab-translator sys . rest)
  (define (cid x)  (second x))
  (define (cn x)   (first x))
  (let-optionals rest ((method 'lsode) (table? #f) (min-v -100) (max-v 100) (step 0.5) )
  (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)))
	     (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))
	     
	     (ionch-info    (nemo:ionch-query sys))
	     (ionchs        (lookup-def 'ion-channels ionch-info))
	     (perm-ions     (map (match-lambda ((comp i e erev) `(,comp ,(matlab-name i) ,(matlab-name e) ,erev)))
				 (lookup-def 'perm-ions ionch-info)))
	     (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map matlab-name (list i in out)))))
				 (lookup-def 'acc-ions ionch-info)))
	     (epools        (lookup-def 'pool-ions ionch-info))
	     (pool-ions     (map (lambda (lst) (map matlab-name lst)) epools))

	     (i-gates       (lookup-def 'i-gates ionch-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-eqs0 (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)))))
			 
			 (cond ((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)))
				  (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)))))
			       
			       ((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)))
				  (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)))
			       (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
			       )))
		     ionchs))
		
	     (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs0))
	     
	     (i-eqs  (fold (lambda (b ax) 
			     (match b 
				    ((and ps ((i e gion) . rst))  
				     (let* ((sum   (if e (sum (map (lambda (b) `(* ,(third b) (- v ,(second b)))) 
								   ps))
						       (sum (map third ps))))
					    (sum0  (rhsexpr/MATLAB sum))
					    (sum1  (canonicalize-expr/MATLAB sum0)))
				       (cons (list i sum1) ax)))
				    
				    ((i e gion)
				     (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)))
	     
	     (pool-ion-i          (map second pool-ions))
	     
	     (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)))))
	     
	     (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))

	       )

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

	   (if (not (null? globals)) (pp indent (global ,(sl\ " " globals))))

	   ;; derivative function
	   (output-dy sysname method globals state-index-map 
		      rate-eq-defs reaction-eq-defs asgn-eq-defs
		      pool-ions mcap i-eqs
		      indent indent+)

	   ;; partial derivative 
	   (output-ddy sysname method globals dfenv state-index-map 
		       rate-eq-defs reaction-eq-defs asgn-eq-defs
		       defuns indent indent+)

	   ;; steady state solver
	   (output-steadystate sysname method globals steady-state-index-map 
			       rate-eq-defs reaction-eq-defs asgn-eq-defs
			       indent indent+)

	   ;; print state vector
	   (output-print-state sysname state-index-map
			       indent indent+)

	   ;; initial values function
	   (output-init sysname globals state-index-map steady-state-index-map 
			const-defs init-eq-defs rate-eq-defs i-eqs pool-ions
			pool-ion-i indent indent+)

	   (pp indent ,nl)
	 
	   ;; user-defined functions
	   (let* ((with       (inexact->exact (round (/ (abs (- max-v min-v)) step))))
		  (define-fn  (make-define-fn table? min-v max-v with)))
	     (for-each (lambda (fndef) 
			 (if (not (member (car fndef) builtin-fns))
			     (apply define-fn (cons indent fndef)))) 
		       defuns))
	 
	   ;; derivatives of user-defined functions
	   (let ((define-dfn  (make-define-dfn dfenv)))
	     (for-each (lambda (fndef) 
			 (if (not (member (car fndef) builtin-fns))
			     (apply define-dfn (cons indent fndef)))) 
		       defuns))


	   )))))
)