salt.scm

 
;;
;; Hybrid dynamical systems modeling.
;;
;; Copyright 2015-2018 Ivan Raikov.
;;
;; This implementation follows the work of Tom Short and his Julia
;; Sims.jl library, which is in turn based on David Broman's MKL
;; simulator and the work of George Giorgidze and Henrik Nilsson in
;; functional hybrid modeling.
;;
;; Following Sims.jl, a nodal formulation is used based on David
;; Broman's thesis:
;; 
;;   David Broman. Meta-Languages and Semantics for Equation-Based
;;   Modeling and Simulation. PhD thesis, Thesis No 1333. Department of
;;   Computer and Information Science, Linköping University, Sweden,
;;   2010.
;;   http://www.bromans.com/david/publ/thesis-2010-david-broman.pdf
;;
;; The type elaboration and environment management routines are based
;; on ideas and code from the Scheme dynamic type inferencer by Andrew
;; Wright.
;;
;; 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 salt 

	(

         parse elaborate simcreate codegen-ODE codegen-ODE/ML codegen-ODE/C 
         math-constant-env math-binop-env math-unop-env 
         model-quantities model-units
         verbose add-trace clear-trace

         empty-env env-lookup extend-env-with-binding env->list
         gen-binding binding-key binding-value

         astdecls? make-astdecls astdecls-decls

         equation-set-model
         equation-set-definitions
         equation-set-discrete-definitions
         equation-set-parameters
         equation-set-fields
         equation-set-externals
         equation-set-externalevs
         equation-set-equations
         equation-set-initial
         equation-set-conditions
         equation-set-pos-responses
         equation-set-neg-responses
         equation-set-functions
         equation-set-nodemap
         equation-set-regimemap
         equation-set-dimenv

         constant parameter
         
         codegen-primop codegen-expr codegen-const-expr codegen-ODE
         name/ML binding->ML stmt->ML value->ML codegen-ODE/ML

         model-time
	 )

	(import scheme (chicken base)
                matchable datatype lalr-driver mathh unitconv with-units fmt fmt-c mathh-consts
                (only srfi-1 first last zip fold fold-right filter filter-map list-tabulate concatenate every delete-duplicates drop-right lset-union)
                (only srfi-4 list->s32vector)
                (only srfi-13 string-null? string-concatenate string-map string<)
		(only (chicken string) ->string conc)
                (only (chicken sort) sort topological-sort)
                (only (chicken format) fprintf sprintf)
                (only (chicken pretty-print) pp)
                (only (chicken port) with-output-to-port port-position port-name)
                (only (chicken condition) print-error-message)
                (chicken memory representation)
		)


(include "units.scm")
(include "parser.scm")
(include "env.scm")
(include "codegen.scm")


(define (s+ . rst) (string-concatenate (map ->string rst)))

(define (salt:warning x . rest)
  (let loop ((port (open-output-string)) (objs (cons x rest)))
    (if (null? objs)
	(begin
	  (newline port)
	  (print-error-message (get-output-string port) 
			       (current-error-port) "salt warning"))
	(begin (display (car objs) port)
	       (display " " port)
	       (loop port (cdr objs))))))


(define (salt:error x . rest)
  (let ((port (open-output-string)))
    (if (port? x)
	(begin
	  (display "[" port)
	  (display (port-name x) port)
	  (display "] " port)))
    (let loop ((objs (if (port? x) rest (cons x rest))))
      (if (null? objs)
	  (begin
	    (newline port)
	    (error 'salt (get-output-string port)))
	  (let ((obj (car objs)))
	    (if (procedure? obj) 
		(with-output-to-port port obj)
		(begin
		  (display obj port)
		  (display " " port)))
	    (loop (cdr objs)))))))


(define verbose (make-parameter 0))

(define trace-locations (make-parameter '()))

(define (add-trace loc)
  (if (not (member loc (trace-locations)))
      (trace-locations (cons loc (trace-locations)))))

(define (clear-trace loc)
  (trace-locations '()))


(define (d loc fstr . args)
  (let ([port (current-error-port)])
    (if (positive? (verbose)) 
	(begin 
          (fprintf port "~A: " loc)
          (apply fprintf port fstr args)
          (flush-output port) ) )))

(define (trace loc fstr . args)
  (let ([port (current-error-port)])
    (if (member loc (trace-locations))
	(begin 
          (fprintf port "~A: " loc)
          (apply fprintf port fstr args)
          (flush-output port) ) )))


(define syntactic-keywords
  '(function if cond and or let else))


(define-record-type variable
  (make-variable  name label value history dim)
  variable?
  (name        variable-name)
  (label       variable-label)
  (value       variable-value)
  (history     variable-history)
  (dim        variable-dim)
  )


(define-record-printer (variable x out)
  (fprintf out "#(variable ~S label=~S value=~S history=~S dim=~A)"
	   (variable-name x)
	   (variable-label x)
	   (variable-value x)
	   (variable-history x)
	   (variable-dim x)
           ))


(define (unknown value label dim)
  (make-variable (gensym 'u) label value #t dim))


(define-record-type free-variable
  (make-free-variable name)
  free-variable?
  (name free-variable-name)
  )


(define-record-printer (free-variable x out)
  (fprintf out "#(free ~S)"
	   (free-variable-name x)
           ))


(define-record-type derivative-variable
  (make-derivative-variable parent)
  derivative-variable?
  (parent      derivative-variable-parent)
  )


(define-record-printer (derivative-variable x out)
  (fprintf out "#(deriv ~S)"
	   (derivative-variable-parent x)
           ))



(define-record-type discrete-variable
  (make-discrete-variable name label value dim)
  discrete-variable?
  (name       discrete-variable-name)
  (label      discrete-variable-label)
  (value      discrete-variable-value)
  (dim        discrete-variable-dim)
  )

(define (discrete value label dim)
  (make-discrete-variable (gensym 'd) label value dim))


(define-record-printer (discrete-variable x out)
  (fprintf out "#(discrete ~S (~a) [~A] = ~A)"
	   (discrete-variable-name x)
	   (discrete-variable-label x)
	   (discrete-variable-dim x)
	   (discrete-variable-value x)
           ))


(define-record-type regime-variable
  (make-regime-variable name)
  regime-variable?
  (name       regime-variable-name)
  )

(define-record-printer (regime-variable x out)
  (fprintf out "#(regime-variable ~S)"
	   (regime-variable-name x)
           ))


(define-record-type reduce-variable
  (make-reduce-variable op parent)
  reduce-variable?
  (op          reduce-variable-op)
  (parent      reduce-variable-parent)
  )


(define-record-printer (reduce-variable x out)
  (fprintf out "#(reduce ~S ~S)"
	   (reduce-variable-op x)
	   (reduce-variable-parent x)
           ))


(define-record-type parameter
  (make-system-parameter name label value dim)
  parameter?
  (name parameter-name)
  (label parameter-label)
  (value parameter-value)
  (dim  parameter-dim)
  )


(define-record-printer (parameter x out)
   (fprintf out "#(parameter ~S (~A) [~A] = ~A)"
 	   (parameter-name x)
 	   (parameter-dim x)
 	   (parameter-label x)
 	   (parameter-value x)
           ))


(define-record-type field
  (make-field name label value dim)
  field?
  (name field-name)
  (label field-label)
  (value field-value)
  (dim  field-dim)
  )


(define-record-printer (field x out)
  (fprintf out "#(field ~S (~A) [~A] = ~A)"
	   (field-name x)
	   (field-dim x)
	   (field-label x)
	   (field-value x)
           ))


(define-record-type declared-constant
  (make-declared-constant label content)
  declared-constant?
  (label declared-constant-label)
  (content declared-constant-content)
  )


(define-record-type constant
  (make-constant type value unit)
  constant?
  (type constant-type)
  (value constant-value)
  (unit constant-unit)
  )



(define-record-printer (constant x out)
  (if (constant-unit x)
      (fprintf out "#(constant ~S = ~A ~A)"
               (constant-type x)
               (constant-value x)
               (constant-unit x)
               )
      (fprintf out "#(constant ~S = ~A)"
               (constant-type x)
               (constant-value x)
               )))


(define-record-type external
  (make-external name label initial dim order)
  external?
  (name external-name)
  (label external-label)
  (initial external-initial)
  (dim external-dim)
  (order external-order)
  )


(define-record-printer (external x out)
  (fprintf out "#(external ~S (~A) [~A] = ~A order: ~A)"
	   (external-name x)
	   (external-dim x)
	   (external-label x)
	   (external-initial x)
	   (or (external-order x) "none")
           ))


(define-record-type externalev
  (make-externalev name label initial dim order link)
  externalev?
  (name externalev-name)
  (label externalev-label)
  (initial externalev-initial)
  (dim externalev-dim)
  (order externalev-order)
  (link externalev-link)
  )


(define-record-printer (externalev x out)
  (fprintf out "#(externalev ~S (~A) [~A] = ~A order: ~A link: ~A)"
	   (externalev-name x)
	   (externalev-label x)
	   (externalev-dim x)
	   (externalev-initial x)
	   (externalev-order x)
	   (or (externalev-link x) "none")
           ))


(define-record-type left-var
  (make-left-var u )
  left-var?
  (u left-var-u)
  )


(define-record-printer (left-var x out)
  (fprintf out "#(left-var ~S)"
	   (left-var-u x)
           ))


(define-record-type ref-var
  (make-make-ref-var  u idx )
  ref-var?
  (u        ref-var-u)
  (idx      ref-var-idx)
  )


(define-record-printer (ref-var x out)
  (fprintf out "#(ref-var ~S [~A])"
	   (ref-var-u x)
	   (ref-var-idx x)
           ))


(define-record-type pair-arg
  (make-pair-arg  fst snd )
  pair-arg?
  (fst      pair-arg-fst)
  (snd      pair-arg-snd)
  )

(define-record-type null-arg
  (make-null-arg)
  null-arg?
  )


(define-record-printer (pair-arg x out)
  (fprintf out "(~A . ~A)"
	   (pair-arg-fst x)
	   (pair-arg-snd x)
           ))


(define-record-type pair-formal
  (make-pair-formal  fst snd )
  pair-formal?
  (fst      pair-formal-fst)
  (snd      pair-formal-snd)
  )


(define-record-printer (pair-formal x out)
  (fprintf out "<~A . ~A>"
	   (pair-formal-fst x)
	   (pair-formal-snd x)
           ))


(define-record-type null-formal
  (make-null-formal )
  null-formal?
  )


(define-record-type var-def
  (make-var-def sym )
  var-def?
  (sym      var-def-sym)
  )


(define-record-printer (var-def x out)
  (fprintf out "#(var-def ~S)"
	   (var-def-sym x)
           ))


(define-record-type function
  (make-function name formals body)
  function?
  (name    function-name)
  (formals function-formals)
  (body    function-body))


(define-record-printer (function x out)
  (fprintf out "#(function ~S (~S) = ~S)"
	   (function-name x)
	   (function-formals x)
	   (function-body x)
           ))


(define-record-type equation
  (make-eq pattern rhs)
  equation?
  (pattern equation-pattern)
  (rhs     equation-rhs))


(define-record-printer (equation x out)
  (fprintf out "#(equation ~S = ~A)"
	   (equation-pattern x)
	   (equation-rhs x)
           ))


(define-record-type initial-equation
  (make-init variable expr)
  initial-equation?
  (variable equation-variable)
  (expr   equation-expr))


(define-record-type event
  (make-event name condition pos neg)
  event?
  (name event-name)
  (condition event-condition)
  (pos event-pos)
  (neg event-neg))

(define-record-printer (event x out)
  (fprintf out "#")
  (pp
   `(event
    (name=,(event-name x))
    (condition=,(event-condition x))
    (pos=,(event-pos x))
    (neg=,(event-neg x))
    )
   out))


(define-record-type evcondition
  (make-evcondition name expr)
  evcondition?
  (name evcondition-name)
  (expr evcondition-expr))

(define-record-printer (evcondition x out)
  (fprintf out "#")
  (pp
   `(evcondition
    (name=,(evcondition-name x))
    (expr=,(evcondition-expr x))
    )
   out))

(define-record-type evresponse
  (make-evresponse name expr)
  evresponse?
  (name evresponse-name)
  (expr evresponse-expr))

(define-record-printer (evresponse x out)
  (fprintf out "#")
  (pp
   `(evresponse
    (name=,(evresponse-name x))
    (expr=,(evresponse-expr x))
    )
   out))


(define-record-type extevlink
  (make-extevlink src dst)
  extevlink?
  (src extevlink-src)
  (dst extevlink-dst))


(define-record-printer (extevlink x out)
  (fprintf out "#")
  (pp
   `(extevlink
    (src=,(extevlink-src x))
    (dst=,(extevlink-dst x))
    )
   out))

(define-record-type transition 
  (make-transition event target condition response)
  transition?
  (event transition-event)
  (target transition-target)
  (condition transition-condition)
  (response transition-response)
  )


(define-record-type structural-event
  (make-structural-event name label regime transition )
  structural-event?
  (name structural-event-name)
  (label structural-event-label)
  (regime structural-event-regime)
  (transition structural-event-transition)
  )


(define-record-type astdecls
  (make-astdecls decls)
  astdecls?
  (decls astdecls-decls)
  )



(define (expr-value x)
  (cond ((variable? x) (variable-value x))
        ((ref-var? x) (vector-ref (ref-var-u x) (ref-var-idx x)))
        ((list? x) (eval `(,(car x) . ,(map expr-value x))))
        (else x)))


(define (variable-identity x)
  (cond ((variable? x)  
         (variable-name x))
        ((derivative-variable? x)   
         (variable-identity (derivative-variable-parent x)))
        ((discrete-variable? x)   
         (discrete-variable-name x))
        ((regime-variable? x)   
         (regime-variable-name x))
        ((reduce-variable? x)   
         (variable-identity (reduce-variable-parent x)))
        ((external? x)   
         (external-name x))
        (else (error 'variable-identity "invalid argument" x))
        ))

(define (reinit e x y)
  (cond ((left-var? x) `(signal.reinit ,e ,x ,y))
        ((variable? x)  (reinit e (make-left-var x) y))
        ((ref-var? x)   (reinit e (make-left-var x) y))
        ((derivative-variable? x)   
         (reinit e (make-left-var (derivative-variable-parent x)) y))
        ((discrete-variable? x)   
         (reinit e (make-left-var x) y))
        ((regime-variable? x)   
         (reinit e (make-left-var x) y))
        ((reduce-variable? x)   
         (reinit e (make-left-var (reduce-variable-parent x)) y))
        (else (error 'reinit "invaliusing using using using d argument to reinit" x))
        ))



(define math-constant-env
  (fold (lambda (k v env)
          (let ((binding (gen-binding k v)))
            (extend-env-with-binding env binding)))
          empty-env
        `(UNITZERO 
          E 1/E E^2 E^PI/4 LOG2E LOG10E LN2 LN3 LNPI LN10 1/LN2 1/LN10 PI PI/2
          PI/4 1/PI 2/PI 2/SQRTPI SQRTPI PI^2 DEGREE SQRT2 1/SQRT2 SQRT3 SQRT5
          SQRT10 CUBERT2 CUBERT3 4THRT2 GAMMA1/2 GAMMA1/3 GAMMA2/3 PHI LNPHI
          1/LNPHI EULER E^EULER SIN1 COS1 ZETA3)
        (cons (make-constant 'number 0.0 'unitbottom)
              (map (lambda (x) (make-constant 'number x unitless))
                   (list e 1/e e^2 e^pi/4 log2e log10e ln2 ln3 lnpi ln10 1/ln2 1/ln10 pi pi/2
                         pi/4 1/pi 2/pi 2/sqrtpi sqrtpi pi^2 degree sqrt2 1/sqrt2 sqrt3 sqrt5
                         sqrt10 cubert2 cubert3 4thrt2 gamma1/2 gamma1/3 gamma2/3 phi lnphi
                         1/lnphi euler e^euler sin1 cos1 zeta3)))
        ))
     

(define math-binop-env
  (let ((function-names `(+ - * / pow max min == < <= > >= random.unifrange))
        (op-names `(signal.add
                    signal.sub signal.mul signal.div signal.pow signal.max signal.min signal.eqnum 
                    signal.lt signal.lte signal.gt signal.gte
                    random.unifrange)))
    (fold (lambda (k v env)
            (let ((binding (gen-binding k v)))
              (extend-env-with-binding env binding)))
          empty-env
          function-names
          (map (lambda (f fn)
                 (make-function
                  (make-free-variable f)
                  (make-pair-formal (make-var-def 'a) (make-pair-formal (make-var-def 'b) (make-null-formal)))
                  `(signal.call ,fn ,(make-pair-arg (make-var-def 'a) (make-pair-arg (make-var-def 'b) (make-null-arg))))
                  ))
               function-names 
               op-names))
    ))


(define math-unop-env
  (let ((function-names
         `(neg abs atan asin acos sin cos exp ln
               sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
               round ceiling floor 
               heaviside random.exponential))
        (op-names
         `(signal.neg signal.abs signal.atan signal.asin signal.acos signal.sin signal.cos signal.exp signal.ln
                      signal.sqrt signal.tan signal.cosh signal.sinh signal.tanh signal.hypot signal.gamma signal.lgamma signal.log10
                      signal.log2 signal.log1p signal.ldexp signal.cube
                      signal.round signal.ceiling signal.floor
                      signal.heaviside 
                      random.exponential)))
    (fold (lambda (k v env)
            (let ((binding (gen-binding k v)))
              (extend-env-with-binding env binding)))
          empty-env
          function-names
          (map (lambda (f fn)
                 (make-function (make-free-variable f)
                                (make-pair-formal (make-var-def 'a) (make-null-formal))
                                `(signal.call ,fn ,(make-pair-arg (make-var-def 'a) (make-null-arg))))
                 )
               function-names
               op-names)
          ))
  )
  

(define math-zop-env
  (let ((function-names
         `(random.uniform))
        (op-names
         `(random.uniform)))
    (fold (lambda (k v env)
            (let ((binding (gen-binding k v)))
              (extend-env-with-binding env binding)))
          empty-env
          function-names
          (map (lambda (f fn)
                 (make-function (make-free-variable f)
                                (make-null-formal)
                                `(signal.call ,fn ,(make-null-arg)))
                 )
               function-names
               op-names)
          ))
  )
  


;; A representation of a flattened model, normally created with
;; `elaborate(model)`. The code generation procedures take an
;; elaborated model as input.
;;
;; Contains the hierarchical equations, flattened equations, flattened
;; initial equations, events, event response functions, and a map of
;; nodes.


(define-record-type equation-set
  (make-equation-set model definitions discrete-definitions parameters fields externals externalevs
                     equations initial conditions pos-responses neg-responses 
                     functions nodemap regimemap dimenv)
  equation-set?
  (model equation-set-model)
  (definitions equation-set-definitions)
  (discrete-definitions equation-set-discrete-definitions)
  (parameters equation-set-parameters)
  (fields equation-set-fields)
  (externals equation-set-externals)
  (externalevs equation-set-externalevs)
  (equations equation-set-equations)
  (initial equation-set-initial)
  (conditions equation-set-conditions)
  (pos-responses equation-set-pos-responses)
  (neg-responses equation-set-neg-responses)
  (functions equation-set-functions)
  (nodemap equation-set-nodemap)
  (regimemap equation-set-regimemap)
  (dimenv equation-set-dimenv)
  )


(define-record-printer (equation-set x out)
  (fprintf out "#")
  (pp
   `(equation-set
    (definitions=,(equation-set-definitions x))    
    (discrete-definitions=,(equation-set-discrete-definitions x))
    (parameters=,(equation-set-parameters x))
    (fields=,(equation-set-fields x))
    (externals=,(equation-set-externals x))
    (externalevs=,(equation-set-externalevs x))
    (equations=,(equation-set-equations x))
    (functions=,(equation-set-functions x))
    (conditions=,(equation-set-conditions x))
    (pos-responses=,(equation-set-pos-responses x))
    (neg-responses=,(equation-set-neg-responses x))
    (nodemap=,(equation-set-nodemap x))
    (regimemap=,(equation-set-regimemap x)))
   out))


;; runtime representation of a simulation object
(define-record-type simruntime
  (make-simruntime eqset cindexmap dindexmap evindexmap rindexmap extindexmap extevindexmap
                   parameters fields defs discrete-defs external-defs externalev-defs externalev-links
                   eqblock condblock posresp negresp)
  simruntime?
  (eqset simruntime-eqset)
  (cindexmap simruntime-cindexmap)
  (dindexmap simruntime-dindexmap)
  (evindexmap simruntime-evindexmap)
  (rindexmap simruntime-rindexmap)
  (extindexmap simruntime-extindexmap)
  (extevindexmap simruntime-extevindexmap)
  (parameters simruntime-parameters)
  (fields simruntime-fields)
  (defs simruntime-definitions)
  (discrete-defs simruntime-discrete-definitions)
  (external-defs simruntime-external-definitions)
  (externalev-defs simruntime-externalev-definitions)
  (externalev-links simruntime-externalev-links)
  (eqblock simruntime-eqblock)
  (condblock simruntime-condblock)
  (posresp simruntime-posresp)
  (negresp simruntime-negresp)
  )


(define-record-printer (simruntime x out)
    (fprintf out "#")
    (pp `(simruntime 
          (cindexmap=,(simruntime-cindexmap x))
          (dindexmap=,(simruntime-dindexmap x))
          (evindexmap=,(simruntime-evindexmap x))
          (extindexmap=,(simruntime-extindexmap x))
          (extevindexmap=,(simruntime-extevindexmap x))
          (parameters=,(simruntime-parameters x))
          (fields=,(simruntime-fields x))
          (defs=,(simruntime-definitions x))
          (external-defs=,(simruntime-external-definitions x))
          (externalev-defs=,(simruntime-externalev-definitions x))
          (externalev-links=,(simruntime-externalev-links x))
          (eqblock=,(simruntime-eqblock x))
          (evblock=,(simruntime-condblock x))
          (posresp=,(simruntime-posresp x))
          (negresp=,(simruntime-negresp x))
          )
        out)
    )
    
(define-record-printer (astdecls x out)
    (fprintf out "#")
    (pp `(ast . ,(astdecls-decls x)) out))


(define model-time (unknown (make-constant 'number 0.0 millisecond) 't millisecond))


(define (map-equations f eqs)
  (map 
   (match-lambda*
    [($ salt#initial-equation s expr) 
     (make-init s (f expr))]
    [($ salt#equation s expr) 
     (make-eq s (f expr))]
    [eqn (error 'map-equation "unknown equation type" eqn)])
   eqs))
      


(define (replace-fixed expr)
  (map
   (lambda (x)
     (cond ((variable? x) (or (variable-value x) x))
           ((list? x) `(,(car x) . ,(map replace-fixed (cdr x))))
           (else x)))
   expr))



;; Equation parsing

(define default-parse-env
  (fold (lambda (ext-env env) (extend-env-with-env env ext-env))
        empty-env
        (list math-constant-env
              math-binop-env
              math-unop-env
              math-zop-env)
        ))


(define (parse decls #!key (env default-parse-env))
    (letrec
        ((parse0
          (lambda (decls)
            (d 'parse "decls = ~A~%" decls)
            (map (lambda (x) 
                   (cond
                    ((astdecls? x) x)
                    (else (parse-declaration env x))
                    ))
                 decls))
          ))
      (make-astdecls (parse0 decls)))
    )


;; Expression resolution

(define (resolve expr env-stack)

  (trace 'resolve "expr: ~A env-stack: ~A~%" 
         expr (env-stack-show env-stack))

  (let recur ((e expr))
    (trace 'resolve "e: ~A~%" e)

    (cond
     
     ((free-variable? e)
      (let ((assoc-var-def (env-stack-lookup (free-variable-name e) env-stack)))
        (trace 'resolve "free-var: e = ~A assoc-var-def = ~A~%" (free-variable-name e) assoc-var-def)
        (if assoc-var-def
            (let recur ((resval (binding-value assoc-var-def)))
              (let ((resval1 (resolve resval env-stack)))
                (if (equal? resval resval1) resval
                    (recur resval1))))
            (let ((assoc-unit (assv (free-variable-name e) (model-units))))
              (trace 'resolve "free-var: e = ~A assoc-unit = ~A~%" (free-variable-name e) assoc-unit)
              (if assoc-unit
                  (make-constant 'number 1.0 (cdr assoc-unit))
                  e)))
        ))

     ((derivative-variable? e)
      (make-derivative-variable (resolve (derivative-variable-parent e) env-stack)))

     ((variable? e)
      (let ((name    (variable-name e))
            (label   (variable-label e))
            (value   (let recur ((resval (variable-value e)))
                       (let ((resval1 (resolve resval env-stack)))
                         (if (equal? resval resval1) resval
                             (recur resval1)))))
            (history (variable-history e))
            (dim     (variable-dim e)))
        (make-variable name label value history dim)))

     ((discrete-variable? e)
      (let ((name    (discrete-variable-name e))
            (label   (discrete-variable-label e))
            (value   (let recur ((resval (discrete-variable-value e)))
                       (let ((resval1 (resolve resval env-stack)))
                         (if (equal? resval resval1) resval
                             (recur resval1)))))
            (dim     (discrete-variable-dim e)))
        (make-discrete-variable name label value dim)))

     ((parameter? e)
      (let ((name    (parameter-name e))
            (label   (parameter-label e))
            (value   (let recur ((resval (parameter-value e)))
                       (let ((resval1 (resolve resval env-stack)))
                         (if (equal? resval resval1) resval
                             (recur resval1)))))
            (dim     (parameter-dim e)))
        (make-system-parameter name label value dim)))

     ((field? e)
      (let ((name    (field-name e))
            (label   (field-label e))
            (value   (let recur ((resval (field-value e)))
                       (let ((resval1 (resolve resval env-stack)))
                         (if (equal? resval resval1) resval
                             (recur resval1)))))
            (dim     (field-dim e)))
        (make-field name label value dim)))

     ((reduce-variable? e)
      (make-reduce-variable
       (resolve (reduce-variable-op e) env-stack)
       (resolve (reduce-variable-parent e) env-stack)))

     ((external? e)
      (let ((name    (external-name e))
            (label   (external-label e))
            (initial (let recur ((resval (external-initial e)))
                       (let ((resval1 (resolve resval env-stack)))
                         (if (equal? resval resval1) resval
                             (recur resval1)))))
            (dim     (external-dim e))
            (order   (external-order e)))
        (make-external name label initial dim order)))

     ((externalev? e)
      (let ((name    (externalev-name e))
            (label   (externalev-label e))
            (initial (let recur ((resval (externalev-initial e)))
                       (let ((resval1 (resolve resval env-stack)))
                         (if (equal? resval resval1) resval
                             (recur resval1)))))
            (dim     (externalev-dim e))
            (order   (externalev-order e))
            (link    (externalev-link e))
            )
        (make-externalev name label initial dim order link)))

     ((pair-arg? e)
      (make-pair-arg (recur (pair-arg-fst e))
                     (recur (pair-arg-snd e))))

     ((function? e)
      (let ((name    (function-name e))
            (formals (function-formals e))
            (body    (function-body e)))
        (let* ((env-stack1 (push-env-stack
                            (function-formal-env formals empty-env) 
                            env-stack))
               (body1 (resolve body env-stack1)))
          (d 'resolve "function: name = ~A body = ~A~%" name body)
          (d 'resolve "function: body1 = ~A~%" body1)
          (make-function name formals body1)
          ))
      )


     ((pair? e)
      (let ((op (car e)) (args (cdr e)))
        (trace 'resolve "op = ~A~%" op)
        (trace 'resolve "args = ~A~%" args)
        (case op
          ((signal.if)   `(,op . ,(map recur args)))
          ((signal.cond) `(,op . ,(map recur args)))
          ((signal.and)  `(,op . ,(map recur args)))
          ((signal.or)   `(,op . ,(map recur args)))
          ((signal.let)  `(,op ,(map (lambda (b) (cons (car b) (recur (cdr b)))) (car args))
                               ,(recur (cadr args))))
          ((signal.call) `(,op . ,(map recur args)))
          (else   (map recur e)))))

     (else e)
     ))
  )

;; Resolution of reinit expressions
(define (resolve-reinit evname expr)

  (d 'resolve-reinit "expr: ~A~%" expr)

  (let recur ((e expr))

    (trace 'resolve-reinit "e: ~A~%" e)

    (cond

     ((pair? e)
      (let ((op (car e)) (args (cdr e)))
        (trace 'resolve "op = ~A~%" op)
        (trace 'resolve "args = ~A~%" args)
        (case op
          ((reinit)   (reinit evname (car args) (cadr args)))
          (else       e))
        ))

     (else e)

     ))
  )

;; Merge the equations of the particular regime into the main list of
;; equations with the requisite conditionals.
(define (merge-regime-eq name regime-eq equations env-stack)
  (d 'merge-regime-eq "regime-eq = ~A~%equations = ~A~%" regime-eq equations)
  (match-let ((($ salt#equation s expr) regime-eq))
             (let (
                   (equations1
                    (let ((s1 (resolve s env-stack))
                          (expr1 (resolve expr env-stack)))
                      (trace 'merge-regime-eq "s1 = ~A~%expr1 = ~A~%" s1 expr1) 
                      (if (null? equations)
                          (list (list s1 `(signal.if ,(make-regime-variable name) 
                                                     ,expr1 ,(make-constant 'number 0.0 'unitbottom))))
                          (match-let
                           (((equations1 found?)
                             (fold (match-lambda* 
                                    [((and eq (s2 expr2)) (eqs found?))
                                     (begin
                                       (trace 'merge-regime-eq "s1 = ~A~%s2 = ~A~%" s1 s2)
                                       (trace 'merge-regime-eq "expr1 = ~A expr2 = ~A~%" expr1 expr2)
                                       (if (equal? (variable-identity s2) (variable-identity s1))
                                           (list (cons (list s1 `(signal.if ,(make-regime-variable name) 
                                                                            ,expr1 ,expr2)) eqs)
                                                 #t)
                                           (list (cons eq eqs) found?)))]
                                     [(eq (eqs found?))
                                      (list (cons eq eqs) found?)])
                                   '(() #f)
                                   equations)))
                           (if found? 
                               equations1 
                               (cons (list s1 `(signal.if ,(make-regime-variable name) 
                                                          ,expr1 ,(make-constant 'number 0.0 'unitbottom)))
                                     equations1))
                           ))
                      ))
                   )
               (trace 'merge-regime-eq "regime-eq = ~A equations1 = ~A~%" regime-eq equations1)
               equations1)
  ))
  

;;
;; Equation elaboration (flattening)
;;
;; The main steps in flattening are:
;;
;; Creates a name resolution environment (parameters,fields,externals,constants,variables,functions).
;; Replaces fixed initial values.
;; Flattens models and populates equation, definition, function lists.
;; Pulls out initial equations and populate initial list.
;; Pulls out event definitions and populates event list.
;; Handle structural events.
;;


(define (elaborate model)

  (define (rename-decls decls)

    (let recur ((decls decls) (ax '()))

      (if (null? decls) (reverse ax)

          (let ((decl (car decls)))
      
            (if (pair? decl)
                
                (recur (append decl (cdr decls)) ax)
                
                (match decl  
                       (($ salt#discrete-variable name label value dim)
                        (let ((decl1 (make-discrete-variable (gensym 'dvar) label value dim)))
                          (recur (cdr decls) (cons decl1 ax))
                          ))
                       (($ salt#variable name label value has-history dim)
                        (let ((decl1 (make-variable (gensym 'var) label value has-history dim)))
                          (recur (cdr decls) (cons decl1 ax))))
                       (($ salt#parameter name label value dim)
                        (let ((decl1 (make-system-parameter (gensym 'p) label value dim)))
                          (recur (cdr decls) (cons decl1 ax))))
                       (($ salt#field name label value dim)
                        (let ((decl1 (make-field (gensym 'fld) label value dim)))
                          (recur (cdr decls) (cons decl1 ax))))
                       (($ salt#external name label initial dim order)
                        (let ((decl1 (make-external (gensym 'ext) label initial dim order)))
                          (recur (cdr decls) (cons decl1 ax))))
                       (($ salt#externalev name label initial dim order link)
                        (let ((decl1 (make-externalev (gensym 'extev) label initial dim order link)))
                          (recur (cdr decls) (cons decl1 ax))))
                       (($ salt#function name formals expr)
                        (recur (cdr decls) (cons decl ax)))
                       (($ salt#structural-event name label regime transition)
                        (let ((decl1 (make-structural-event (gensym 'sevn) label regime transition)))
                          (recur (cdr decls) (cons decl1 ax))))
                       (($ salt#event name condition pos neg)
                        (let ((decl1 (make-event (gensym 'evn) condition pos neg)))
                          (recur (cdr decls) (cons decl1 ax))))
                       (else (recur (cdr decls) (cons decl ax)))
                       )
                ))
          ))
    )

  (define (model-env decls)

    (let recur ((decls decls) (env empty-env))

      (if (null? decls) env

          (let ((decl (car decls)))
      
            (if (pair? decl)
                
                (recur (append decl (cdr decls)) env)
                
                (match decl  
                       (($ salt#discrete-variable name label value dim)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label decl))
                               ))
                       (($ salt#variable name label value has-history dim)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label decl))
                               ))
                       (($ salt#parameter name label value dim)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label decl))
                               ))
                       (($ salt#declared-constant label content)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label content))
                               ))
                       (($ salt#field name label value dim)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label decl))
                               ))
                       (($ salt#external name label initial dim order)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label decl))
                               ))
                       (($ salt#externalev name label initial dim order link)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label decl))
                               ))
                       (($ salt#function name formals expr)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding name decl))
                               ))
                       (($ salt#structural-event name label regime transition)
                        (recur (cdr decls)
                               (extend-env-with-binding env (gen-binding label (make-regime-variable name)))
                               ))
                       (($ unitconv#unit name dims factor abbrevs)
                        (let* ((dimsval (eval-dim-expr dims #f (model-quantities)))
                               (factorval (eval-unit-factor factor (model-units)))
                               (factordim (dim-unit-factor factor (model-units)))
                               (unitdecl (make-unit name (make-quantity (gensym 'dim) dimsval) factorval abbrevs)))
                          (if (or (zero? factordim) (= dimsval factordim))
                              (model-units
                               (cons (cons name unitdecl) (model-units)))
                              (error 'elaborate "unit definition dimension mismatch" dims factor))
                          (recur (cdr decls) env)))
                       (else (recur (cdr decls) env))
                       ))
            ))
      ))

  (let ((decls (rename-decls (astdecls-decls model))))

    (let recur (
                (entries        decls)
                (env-stack      (extend-env-stack-with-binding
                                 (push-env-stack (model-env decls) empty-env-stack)
                                 (gen-binding (variable-label model-time) model-time)
                                 ))
                (definitions    '()) 
                (discrete-definitions '()) 
                (parameters     '())
                (fields         '())
                (externals      '())
                (externalevs    '())
                (equations      '())
                (initial        '())
                (conditions     '())
                (pos-responses  '())
                (neg-responses  '())
                (functions      '())
                (nodemap        empty-env)
                (regimemap      empty-env)
                (dimenv         empty-env)
                )
      
    (if (null? entries)
        (begin
          (make-equation-set model
                             (reverse definitions)
                             (reverse discrete-definitions)
                             (reverse parameters)
                             (reverse fields)
                             (reverse externals)
                             (reverse externalevs)
                             (reverse equations)
                             (map-equations replace-fixed initial)
                             (reverse conditions)
                             (reverse pos-responses)
                             (reverse neg-responses)
                             functions
                             nodemap
                             regimemap
                             dimenv
                             )
          )

        (let ((en (car entries)))
          (trace 'elaborate "en = ~A~%" en)

          (if (astdecls? en)

              (let ((decls1 (rename-decls (astdecls-decls en))))

                (trace 'elaborate "decls1 = ~A~%" decls1)
                (trace 'elaborate "model-env decls1 = ~A~%" (env->list (model-env decls1)))
                (trace 'elaborate "env = ~A~%" (env-stack-show (push-env-stack (model-env decls1) env-stack)))


                (recur (append decls1 (cons 'pop-env-stack (cdr entries)))
                       (push-env-stack (model-env decls1) env-stack)
                       definitions discrete-definitions parameters fields externals externalevs
                       equations initial conditions pos-responses neg-responses 
                       functions nodemap regimemap dimenv))

              (match en  

                     ('pop-env-stack
                      (begin
                        (trace 'elaborate "pop-env-stack: env = ~A~%" (env-stack-show env-stack))
                        (trace 'elaborate "pop-env-stack: (pop-env-stack env) = ~A~%" (env-stack-show (pop-env-stack env-stack)))
                        (recur (cdr entries) (pop-env-stack env-stack)
                               definitions discrete-definitions parameters fields externals externalevs
                               equations initial conditions pos-responses neg-responses 
                               functions nodemap regimemap dimenv
                               )
                        ))

                     (($ salt#variable name label value has-history dim)
                      (let* ((resolved-value (resolve value env-stack))
                             (en1 (make-variable name label resolved-value has-history dim)))
                        (trace 'elaborate "variable: name = ~A label = ~A value = ~A~%" name label value)
                        (recur (cdr entries) env-stack
                               (cons (cons name resolved-value) definitions)
                               discrete-definitions
                               parameters fields externals externalevs equations initial conditions 
                               pos-responses neg-responses functions
                               (extend-env-with-binding nodemap (gen-binding name en1))
                               regimemap
                               (extend-env-with-binding dimenv (gen-binding name dim))
                               )
                        ))

                     (($ salt#discrete-variable name label value dim)
                      (let* ((resolved-value (resolve value env-stack))
                             (en1 (make-discrete-variable name label resolved-value dim)))
                        (recur (cdr entries) env-stack
                               definitions
                               (cons (cons name resolved-value) discrete-definitions)
                               parameters fields externals externalevs equations initial conditions 
                               pos-responses neg-responses functions
                               (extend-env-with-binding nodemap (gen-binding name en1))
                               regimemap
                               (extend-env-with-binding dimenv (gen-binding name dim))
                               )
                        ))

                     (($ salt#parameter name label value dim)
                      (trace 'elaborate "parameter: label = ~A value = ~A~%" label value)
                      (let* ((resolved-value (resolve value env-stack))
                             (en1 (make-system-parameter name label resolved-value dim)))
                        (recur (cdr entries) env-stack
                               definitions discrete-definitions
                               (cons (cons name resolved-value) parameters) fields externals externalevs
                               equations initial conditions pos-responses neg-responses 
                               functions
                               (extend-env-with-binding nodemap (gen-binding name en1))
                               regimemap
                               (extend-env-with-binding dimenv (gen-binding name dim))
                               )
                        ))

                     (($ salt#declared-constant label content)
                      (recur (cdr entries) env-stack
                             definitions discrete-definitions
                             parameters fields externals externalevs
                             equations initial conditions pos-responses neg-responses 
                             functions
                             nodemap
                             regimemap
                             dimenv
                             )
                      )

                     (($ salt#field name label value dim)
                      (trace 'elaborate "field: label = ~A value = ~A~%" label value)
                      (let* ((resolved-value (resolve value env-stack))
                             (en1 (make-field name label resolved-value dim)))
                        (recur (cdr entries) env-stack
                               definitions discrete-definitions
                               parameters (cons (cons name resolved-value) fields) externals externalevs
                               equations initial conditions pos-responses neg-responses 
                               functions
                               (extend-env-with-binding nodemap (gen-binding name en1))
                               regimemap
                               (extend-env-with-binding dimenv (gen-binding name dim))
                               )
                        ))

                     (($ salt#external name label initial-value dim order)
                      (trace 'elaborate "external: label = ~A initial-value = ~A~%" label initial-value)
                      (let* ((resolved-initial (resolve initial-value env-stack))
                             (en1 (make-external name label resolved-initial dim order)))
                        (recur (cdr entries) env-stack
                               definitions discrete-definitions
                               parameters fields (cons (cons name resolved-initial) externals) externalevs
                               equations initial conditions pos-responses neg-responses 
                               functions
                               (extend-env-with-binding nodemap (gen-binding name en1))
                               regimemap
                               (extend-env-with-binding dimenv (gen-binding name dim))
                               )
                        ))

                     (($ salt#externalev name label initial-value dim order link)
                      (trace 'elaborate "externalev: label = ~A initial-value = ~A~%" label initial-value)
                      (let* ((resolved-initial (resolve initial-value env-stack))
                             (resolved-link (if (not link) #f
                                                (make-extevlink (resolve (extevlink-src link) env-stack)
                                                                (resolve (extevlink-dst link) env-stack))))
                             (en1 (make-externalev name label resolved-initial dim order resolved-link)))
                        (recur (cdr entries) env-stack
                               definitions discrete-definitions
                               parameters fields externals (cons (cons name resolved-initial) externalevs)
                               equations initial conditions pos-responses neg-responses 
                               functions
                               (extend-env-with-binding nodemap (gen-binding name en1))
                               regimemap
                               (extend-env-with-binding dimenv (gen-binding name dim))
                               )
                        ))

                     (($ salt#equation s expr)
                      (trace 'elaborate "equation: unresolved s = ~A~%" s)
                      (trace 'elaborate "equation: unresolved rhs = ~A~%" expr)
                      (trace 'elaborate "equation: resolved rhs = ~A~%" (resolve expr env-stack))
                      (trace 'elaborate "equation: resolved s = ~A~%" (resolve s env-stack))
                      (recur (cdr entries) env-stack
                             definitions discrete-definitions parameters fields externals externalevs
                             (cons (list (resolve s env-stack) (resolve expr env-stack)) equations)
                             initial conditions pos-responses neg-responses 
                             functions nodemap regimemap dimenv
                             )
                      )
                     (($ salt#initial-equation s expr)
                      (recur (cdr entries) env-stack
                             definitions discrete-definitions parameters fields  externals externalevs
                             equations (cons (cons s (resolve expr env-stack)) initial) 
                             conditions pos-responses neg-responses 
                             functions nodemap regimemap dimenv
                             )
                      )
                     (($ salt#structural-event name label regime ($ salt#transition event target condition pos))
                      (let (
                            (condition-index (length conditions))
                            (regime-index (env-size regimemap))
                            )
                        (trace 'elaborate "structural-event: name = ~A label = ~A event = ~A ~%" name label event)
                        (trace 'elaborate "structural-event: target = ~A~%" target)
                        (recur (cdr entries) env-stack
                               definitions discrete-definitions parameters fields externals externalevs
                               (fold (lambda (eq ax) (merge-regime-eq name eq ax env-stack)) equations regime) initial
                               ;;(cons (make-evcondition event (resolve condition env-stack)) conditions)
                               (let ((cond-eq `(signal.if ,(make-regime-variable name) 
                                                          ,(resolve condition env-stack) ,(make-constant 'number -1.0 'unitbottom))))
                                 (cons (make-evcondition event cond-eq) conditions))
                               (append (cons
                                        (make-evresponse event (resolve-reinit event `(reinit ,(make-regime-variable name) #f)))
                                         (cons
                                          (make-evresponse event (resolve-reinit event `(reinit ,(resolve target env-stack) #t)))
                                          (map
                                           (lambda (x) 
                                             (make-evresponse event (resolve-reinit event (resolve x env-stack))))
                                           pos)))
                                         pos-responses)
                               neg-responses
                               functions 
                               nodemap
                               (extend-env-with-binding regimemap (gen-binding name (list regime-index condition-index)))
                               dimenv
                               ))
                      )
                     (($ salt#event name condition pos neg)
                      (recur (cdr entries) env-stack
                             definitions discrete-definitions parameters fields externals externalevs
                             equations initial 
                             (cons (make-evcondition name (resolve condition env-stack)) conditions) 
                             (append (map (lambda (x) 
                                            (make-evresponse name (resolve-reinit name (resolve x env-stack)))) pos) 
                                     pos-responses)
                             (or (and neg (append (map (lambda (x) 
                                                         (make-evresponse name (resolve-reinit name (resolve x env-stack)))) neg)
                                                  neg-responses))
                                 neg-responses)
                             functions nodemap regimemap dimenv
                             ))
                     (($ salt#function name formals expr)
                      (recur (cdr entries) env-stack
                             definitions discrete-definitions parameters fields externals externalevs
                             equations initial conditions pos-responses neg-responses 
                             (cons (list name formals (resolve expr env-stack)) functions)
                             nodemap regimemap dimenv
                             ))
                     ((or (? unit?) (? quantity?))
                      (recur (cdr entries) env-stack
                             definitions discrete-definitions parameters fields externals externalevs
                             equations initial conditions pos-responses neg-responses 
                             functions nodemap regimemap dimenv
                             ))

                     (decl
                      (error 'elaborate "unknown declaration type" decl))
                     ))
          ))
    ))
  )



(define (function-formal-env formals env)
  (d 'function-formal-env "formals = ~A~%" formals)
  (match formals
         (($ salt#pair-formal x rst)
          (let ((b (gen-binding (var-def-sym x) x)))
            (function-formal-env rst (extend-env-with-binding env b))))
         (($ salt#null-formal)
          env)
         ))


(define (function-call-env formals args env)
  (d 'function-call-env "formals = ~A args = ~A~%" formals args)
  (match formals
         (($ salt#pair-formal x frst)
          (match args
                 ((or ($ salt#pair-arg xv arst) (xv . arst))
                  (let ((b (gen-binding (var-def-sym x) xv)))
                    (function-call-env frst arst (extend-env-with-binding env b))))
                 ((or ($ salt#null-arg) ())
                  (error 'function-call-env
                         "function arity mismatch"))
                 ))
         (($ salt#null-formal)
          (match args
                 ((or ($ salt#null-arg) ()) env)
                 (else 
                  (error 'function-call-env
                         "function arity mismatch"))
                 ))
         ))


(define (subst-symbol sym env-or-env-stack)
  (let ((assoc-var-def   (if (env-stack? env-or-env-stack)
                             (env-stack-lookup sym env-or-env-stack)
                             (env-lookup sym env-or-env-stack))))
    (if assoc-var-def
        (binding-value assoc-var-def) 
        sym)))
        

(define (subst-symbol* sym env-or-env-stack)
  (let ((assoc-var-def   (if (env-stack? env-or-env-stack)
                             (env-stack-lookup sym env-or-env-stack)
                             (env-lookup sym env-or-env-stack))))
    (if assoc-var-def
        (binding-value assoc-var-def) 
        (error 'subst-symbol* "symbol not found" sym))))
        

(define (subst-if args env-stack)
  (map (lambda (x) (subst-expr x env-stack)) args))

(define (subst-reinit args env-stack)
  (cons (car args) (map (lambda (x) (subst-expr x env-stack)) (cdr args))))

(define (subst-cond args env-stack)
  (map (lambda (x) (subst-expr x env-stack)) args))

(define (subst-and args env-stack)
  (map (lambda (x) (subst-expr x env-stack)) args))

(define (subst-or args env-stack)
  (map (lambda (x) (subst-expr x env-stack)) args))

(define (subst-let args env-stack)
  (let ((lbnds (car args))
        (body (cadr args)))
    (list (map (lambda (lb) (cons (car lb) (subst-expr (cdr lb) env-stack))) lbnds)
          (subst-expr body env-stack))))

(define (subst-pair-arg arg env-stack)
  (match arg
         (($ salt#pair-arg fst snd)
          (cons (subst-expr fst env-stack)
                (subst-pair-arg snd env-stack)))
         (($ salt#null-arg)
          (list))))


(define (subst-function-call f arg env-stack)
  (d 'subst-function-call "f = ~A arg=~A~%" f arg)
  (if (symbol? f)
      (let ((args (subst-pair-arg arg env-stack)))
        `(signal.primop ,f . ,args))
      (match f (($ salt#function name formals body)
                (let ((args (subst-pair-arg arg env-stack)))
                  (let ((fenv (function-call-env formals args empty-env)))
                    (trace 'subst-function-call "fenv = ~A~%" (map car (env->list fenv)))
                    (let ((expr1 (subst-expr body (push-env-stack fenv env-stack))))
                      (trace 'subst-function-call "expr1 = ~A~%" expr1)
                      expr1))
                  ))
             (($ salt#free-variable name)
              (error 'subst-function-call "undefined function" name))
             (else
              (error 'subst-function-call "unknown callable object" f))
             )
      ))


(define (subst-expr e env-stack)
  (d 'subst-expr "e = ~A~%" e)
  (cond ((symbol? e)
         (subst-symbol e env-stack))
        ((pair? e)
         (let ((op (car e)) (args (cdr e)))
           (trace 'subst-expr "op = ~A args = ~A~%" op args)
           (case op
             ((signal.reinit) (cons 'signal.reinit (subst-reinit args env-stack)))
             ((signal.if)     (cons 'signal.if (subst-if args env-stack)))
             ((signal.cond)   (cons 'signal.cond (subst-cond args env-stack)))
             ((signal.and)    (cons 'signal.and (subst-and args env-stack)))
             ((signal.or)     (cons 'signal.or (subst-or args env-stack)))
             ((signal.let)    (cons 'signal.let (subst-let args env-stack)))
             ((signal.call)  
              (let ((f (car args))
                    (fargs (cadr args)))
                (trace 'subst-expr "f = ~A fargs = ~A~%" f args)
                (subst-function-call f fargs env-stack)
                ))
             ((signal.primop)  
              (let ((op (car args))
                    (opargs (cdr args)))
                `(signal.primop ,op . ,(map (lambda (x) (subst-expr x env-stack)) opargs))))
             (else e)
             )))
        ((var-def? e)
         (subst-symbol* (var-def-sym e) env-stack))
        ((free-variable? e)
         (subst-symbol (free-variable-name e) env-stack))
        ((left-var? e)
         (subst-expr (left-var-u e) env-stack))
        (else e)
        ))

(define (units-unify u1 u2)
  (cond ((and (unit? u1) (unit? u2))
         (unit-equal? u1 u2))
        ((and (equal? u1 'unitbottom) (equal? u2 'unitbottom)) #t)
        ((and (unit? u2) (equal? u1 'unitbottom)) #t)
        ((and (unit? u1) (equal? u2 'unitbottom)) #t)
        (else #f)))

(define (units-reinit args env)
  (let ((x (cadr args))
        (y (caddr args)))
    (let ((xu (expr-units x env))
          (yu (expr-units y env)))
      (if (units-unify xu yu)
          xu
          (error 'units-reinit "units mismatch in reinit" args))
      ))
  )

        
(define (units-if args env)
  (match args
         ((c x y)
          (let ((cu (expr-units c env))
                (xu (expr-units x env))
                (yu (expr-units y env)))
            (trace 'units-if "cu = ~A xu = ~A yu = ~A~%" cu xu yu)
            ;; TODO: this is not necessary if regime equations take
            ;; into account units of quantities
            (cond 
             ((units-unify xu yu) (if (unit? xu) xu yu))
             ((and (not (unitless? xu)) (unitless? yu)) xu)
             ((and (unitless? xu) (not (unitless? yu))) yu)
             (else (error 'units-reinit "units mismatch in if" args)))
          ))
  (else (error 'units-reinit "invalid arguments to if" args))
  ))


(define (units-let args env)
  (let ((lbnds (car args))
        (body (cadr args)))
    (let ((env1
           (fold
            (lambda (lb env) 
              (extend-env-with-binding 
               env (gen-binding (car lb) (expr-units (cdr lb) env))))
            env
            lbnds)))
          (expr-units body env1))))


(define (units-pair-arg arg env)
  (match arg
         (($ salt#pair-arg fst snd)
          (cons (expr-units fst env)
                (expr-units snd env)))
         (($ salt#null-arg)
          (list))))


(define (units-primop f args units-args)
  (d 'units-primop "f = ~A units-args = ~A args = ~A~%" f units-args args)
  (let ((u (case f 
             ((signal.add signal.sub) 
              (if (units-unify (car units-args) (cadr units-args))
                  (car units-args) 
                  (error 'units-primop "units mismatch" f units-args
                         (quantity-int (unit-dims (car units-args)))
                         (quantity-int (unit-dims (cadr units-args))))))
             ((signal.mul) 
              (unit* (car units-args) (cadr units-args) ))
             ((signal.div) 
              (unit/ (car units-args) (cadr units-args) ))
             ((signal.pow) 
              (unit-expt (car units-args) (constant-value (cadr args) )))
             ((signal.neg) 
              (car units-args))
             (else unitless))))
    (trace 'units-primop "f = ~A args = ~A units-args = ~A units-args-dims = ~A u = ~A dims(u) = ~A~%" 
           f args units-args (map (compose quantity-int unit-dims) units-args) u (quantity-int (unit-dims u)))
    u))


(define (units-function-call f arg env)
  (d 'units-function-call "f = ~A arg = ~A~%" f arg)
  (if (symbol? f)
      (let ((units-args (units-pair-arg arg env)))
        (trace 'units-function-call "units-args = ~A~%" units-args)
        (units-primop f arg units-args))
      (match-let
       ((($ salt#function name formals body) f))
       (let ((fenv (function-call-env formals arg empty-env)))
         (trace 'units-function-call "fenv = ~A~%" fenv)
         (let ((u (expr-units 
                   (subst-expr body (push-env-stack fenv empty-env-stack))
                   env)))
           (trace 'units-function-call "u = ~A~%" u)
           u)
         ))
      ))


(define (expr-units e env)
  (d 'expr-units "e = ~A env = ~A~%" e (env->list env))
  (cond ((symbol? e)
         (let ((u (assv e (model-units))))
           (trace 'expr-units "(symbol? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           (or (and u (cdr u)) (env-lookup e env) unitless)))
        ((pair? e)
         (let ((op (car e)) (args (cdr e)))
           (let ((u (case op
                      ((signal.primop) 
                       (units-primop (car args) (cdr args) 
                                     (map (lambda (e) (expr-units e env)) (cdr args))))
                      ((signal.reinit) 
                       (units-reinit args env))
                      ((signal.if)     
                       (units-if args env))
                      ((signal.let)    
                       (units-let args env))
                      ((signal.call)  
                       (let ((f (car args))
                             (fargs (cadr args)))
                         (units-function-call f fargs env)
                         ))
                      (else e)
                      )))
             (trace 'expr-units "(pair? e) => op = ~A args = ~A u = ~A dims(u) = ~A~%" e args u (and (unit? u) (quantity-int (unit-dims u))))
             u)))
        ((var-def? e)
         (let ((u (expr-units (subst-symbol (var-def-sym e) env) env)))
           (trace 'expr-units "(var-def? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        ((free-variable? e)
         (error 'expr-units "unknown variable" (free-variable-name e)))
        ((variable? e)
         (let ((u (expr-units (variable-value e) env)))
           (trace 'expr-units "(variable? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        ((discrete-variable? e)
         (let ((u (expr-units (discrete-variable-value e) env)))
           (trace 'expr-units "(discrete-variable? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        ((parameter? e)
         (let ((u (expr-units (parameter-value e) env)))
           (trace 'expr-units "(parameter? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        ((field? e)
         (let ((u (expr-units (field-value e) env)))
           (trace 'expr-units "(field? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        ((external? e)
         (let ((u (expr-units (external-initial e) env)))
           (trace 'expr-units "(external? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        ((left-var? e)
         (let ((u (expr-units (left-var-u e) env)))
           (trace 'expr-units "(left-var? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        ((constant? e)
         (let ((u (constant-unit e)))
           (trace 'expr-units "(constant? e) => e = ~A u = ~A dims(u) = ~A~%" e u (and (unit? u) (quantity-int (unit-dims u))))
           u))
        (else e))
  )
         

(define (reduce-expr expr indexmaps )
  (let (
        (pindexmap  (alist-ref 'pindexmap indexmaps))
        (cindexmap  (alist-ref 'cindexmap indexmaps))
        (dindexmap  (alist-ref 'dindexmap indexmaps))
        (rindexmap  (alist-ref 'rindexmap indexmaps))
        (evindexmap (alist-ref 'evindexmap indexmaps))
        (extindexmap (alist-ref 'extindexmap indexmaps))
        (extevindexmap (alist-ref 'extevindexmap indexmaps))
        )
  (d 'reduce-expr "expr = ~A~%" expr)
  (match expr

         (($ salt#pair-arg fst snd)
          (make-pair-arg
           (reduce-expr fst indexmaps)
           (reduce-expr snd indexmaps)))

         (($ salt#derivative-variable y)
          (let ((yindex (env-lookup (variable-name y) cindexmap)))
            (if (not yindex)
                (error 'reduce-expr "variable not in index" y)
                `(getindex dy ,(cdr yindex)))
            ))
         
         (($ salt#left-var y)
          (cond ((variable? y)
                 (let ((yindex (env-lookup (variable-name y) cindexmap)))
                   (if (not yindex)
                       (error 'reduce-expr "variable not in index" y)
                       `(getindex y ,(cdr yindex)))
                   ))
                ((discrete-variable? y)
                 (let ((yindex (env-lookup (discrete-variable-name y) dindexmap)))
                   (if (not yindex)
                       (error 'reduce-expr "variable not in index" y)
                       `(getindex d ,(cdr yindex)))
                   ))
                ((regime-variable? y)
                 (let ((yindex (env-lookup (regime-variable-name y) rindexmap)))
                   (if (not yindex)
                       (error 'reduce-expr "regime variable not in index" y)
                       `(getindex r ,(car (cdr yindex))))
                   ))
                (else 
                 (error 'reduce-expr "unknown left variable type" expr))
                ))
         
         (($ salt#variable name label value has-history dim)
          (let ((yindex (env-lookup name cindexmap)))
            (if (not yindex)
                (if (equal? name (variable-name model-time))
                    label
                    (error 'reduce-expr "variable not in index" name))
                `(getindex y ,(cdr yindex)))
            ))

         (($ salt#reduce-variable op y)
          (cond 
           ((variable? y)
            (let ((yindex (env-lookup (variable-name y) cindexmap)))
              (if (not yindex)
                  (error 'reduce-expr "variable not in index" y)
                  `(getindex y ,(cdr yindex)))
              ))
           ((external? y)
            (let ((yindex (env-lookup (external-name y) extindexmap)))
              (if (not yindex)
                  (error 'reduce-expr "variable not in external index" y)
                  `(getindex yext ,(cdr yindex)))
              ))
           (else
            (error 'reduce-expr "unsupported reduce variable type" y))
           ))
         
         (($ salt#discrete-variable name label value dim)
          (let ((dindex (env-lookup name dindexmap)))
            (if (not dindex)
                (error 'reduce-expr "discrete variable not in index" name)
                `(getindex d ,(cdr dindex)))
            ))

         (($ salt#regime-variable name)
          (let ((rindex (env-lookup name rindexmap)))
            (if (not rindex)
                (error 'reduce-expr "regime variable not in index" name)
                `(getindex r ,(car (cdr rindex))))
            ))

         (($ salt#parameter name label value dim)
          (let ((yindex (env-lookup name pindexmap)))
            (trace 'reduce-expr "pindexmap = ~A~%" pindexmap)
            (if (not yindex)
                (error 'reduce-expr "parameter not in index" name)
                `(getindex p ,(cdr yindex)))
            ))

         (($ salt#field name label value dim)
          (let ((yindex (env-lookup name pindexmap)))
            (trace 'reduce-expr "pindexmap = ~A~%" pindexmap)
            (if (not yindex)
                (error 'reduce-expr "parameter not in index" name)
                `(getindex fld ,(cdr yindex)))
            ))

         (($ salt#external name label initial dim order)
          (let ((yindex (env-lookup name extindexmap)))
            (trace 'reduce-expr "extindexmap = ~A~%" extindexmap)
            (if (not yindex)
                (error 'reduce-expr "external not in index" name)
                `(getindex yext ,(cdr yindex)))
            ))

         (($ salt#externalev name label initial dim order link)
          (let ((yindex (env-lookup name extevindexmap)))
            (trace 'reduce-expr "extevindexmap = ~A~%" (env->list extevindexmap))
            (if (not yindex)
                (error 'reduce-expr "external event not in index" name)
                `(signal.primop signal.sub ,(variable-label model-time) (getindex extev ,(cdr yindex))))
            ))

         (('signal.if ifcond ift iff)
          (let ((expr1 (reduce-expr ifcond indexmaps))
                (expr2 (reduce-expr ift indexmaps))
                (expr3 (reduce-expr iff indexmaps)))
            (if (and (constant? expr2) (constant? expr3)
                     (= (constant-value expr2) (constant-value expr3)))
                ift `(signal.if ,expr1 ,expr2 ,expr3))
            ))

         (('signal.let lbnds body)
          (subst-expr
           `(signal.let 
             ,(map (lambda (lb) (cons (car lb) (reduce-expr (cdr lb) indexmaps))) lbnds)
             ,(reduce-expr body indexmaps))
           empty-env-stack))

         (('signal.call f args)
          (reduce-expr
           (subst-expr 
            `(signal.call 
              ,f ,(reduce-expr args indexmaps))
            empty-env-stack)
           indexmaps))

         (('signal.reinit e x y)
          (let ((eindex (env-lookup e evindexmap)))
            (subst-expr 
             `(signal.reinit (getindex c ,(cdr eindex))
                             ,(reduce-expr x indexmaps) 
                             ,(reduce-expr y indexmaps))
             empty-env-stack)))

         ((op . args)
           (subst-expr
            (cons op (map (lambda (x) (reduce-expr x indexmaps)) args))
            empty-env-stack))

         (else (subst-expr expr empty-env-stack))

         ))
)


;; reduce expr for constant expressions
(define (reduce-constant-expr expr indexmaps)
  (let (
        (pindexmap  (alist-ref 'pindexmap indexmaps))
        (cindexmap  (alist-ref 'cindexmap indexmaps))
        (dindexmap  (alist-ref 'dindexmap indexmaps))
        (evindexmap (alist-ref 'evindexmap indexmaps))
        )

  (d 'reduce-constant-expr "expr = ~A~%" expr)
  (match expr

         (($ salt#pair-arg fst snd)
          (make-pair-arg
           (reduce-constant-expr fst indexmaps)
           (reduce-constant-expr snd indexmaps)))

         (($ salt#variable name label value has-history dim)
          (reduce-constant-expr value indexmaps))

         (($ salt#parameter name label value dim)
          (reduce-constant-expr value indexmaps))

         (($ salt#field name label value dim)
          (reduce-constant-expr value indexmaps))

         (('signal.let lbnds body)
          (subst-expr
           `(signal.let 
             ,(map (lambda (lb) (cons (car lb) (reduce-constant-expr (cdr lb) indexmaps))) lbnds)
             ,(reduce-constant-expr body indexmaps))
           empty-env-stack))

         (('signal.call f args)
          (subst-expr 
           `(signal.call ,f ,(reduce-constant-expr args indexmaps))
           empty-env-stack))

         ((op . args)
          (subst-expr
           (cons op (map (lambda (x) (reduce-constant-expr x indexmaps)) args))
           empty-env-stack))

         (else expr)
         ))
)


(define (reduce-eq eq indexmaps unit-env)

  (define (make-call fn a b)
    `(signal.call ,fn ,(make-pair-arg a (make-pair-arg b (make-null-arg)))))

  (let (
        (pindexmap  (alist-ref 'pindexmap indexmaps))
        (cindexmap  (alist-ref 'cindexmap indexmaps))
        (dindexmap  (alist-ref 'dindexmap indexmaps))
        (evindexmap (alist-ref 'evindexmap indexmaps))
        (rindexmap  (alist-ref 'rindexmap indexmaps))
        (extindexmap (alist-ref 'extindexmap indexmaps))
        (extevindexmap (alist-ref 'extevindexmap indexmaps))
        )
    (d 'reduce-eq "eq = ~A~%" eq)

    (match eq

           ((($ salt#derivative-variable y) rhs)
            (let ((yindex (env-lookup (variable-name y) cindexmap))
                  (dim (variable-dim y)))
              (if (not yindex)
                  (error 'reduce-eq "variable not in index" y)
                  (let ((expr (reduce-expr rhs indexmaps))
                        (ddim (if (equal? (quantity-int dim) (quantity-int Unity))
                                  (quantity-int dim)
                                  (- (quantity-int dim) (quantity-int Time)))))
                    (trace 'reduce-eq "dim = ~A ddim = ~A rhs = ~A ~%" dim ddim rhs)
                    (let ((rhs-units (expr-units rhs unit-env)))
                      (trace 'reduce-eq "rhs-units = ~A dims(rhs-units) = ~A ~%" rhs-units (unit-dims rhs-units))
                      (if (equal? ddim (quantity-int (unit-dims rhs-units)))
                          `(setindex dy_out ,(cdr yindex) ,expr)
                          (error 'reduce-eq "dimension mismatch in rhs" 
                                 (variable-label y) ddim (unit-dims rhs-units)))
                    ))
              ))
           )
           
           ((($ salt#variable name label value has-history dim) rhs)
            (let ((yindex (env-lookup name cindexmap)))
              (if (not yindex)
                  (error 'reduce-eq "variable not in index" name)
                  (let ((expr (reduce-expr rhs indexmaps)))
                    (let ((rhs-units (expr-units rhs unit-env)))
                      (trace 'reduce-eq "rhs = ~A~%" rhs)
                      (trace 'reduce-eq "dim = ~A dims(rhs-units) = ~A ~%" 
                             dim (and (unit? rhs-units) (unit-dims rhs-units)))
                      (if (or (equal? rhs-units 'unitbottom) (equal? dim (unit-dims rhs-units)))
                          `(setindex y_out ,(cdr yindex) ,expr)
                          (error 'reduce-eq "dimension mismatch in rhs" 
                                 name label
                                 dim (unit-dims rhs-units)))
                          ))
                    ))
              )

           ((($ salt#reduce-variable op y) rhs)
            (cond
             ((variable? y)
              (let ((yindex (env-lookup (variable-name y) cindexmap))
                    (dim (variable-dim y)))
                (if (not yindex)
                    (error 'reduce-eq "variable not in index" y)
                    (let ((expr (reduce-expr (make-call op `(getindex y ,(cdr yindex)) rhs) indexmaps)))
                      (let ((rhs-units (expr-units rhs unit-env)))
                        (trace 'reduce-eq "dim = ~A dims(rhs-units) = ~A ~%" dim (unit-dims rhs-units))
                        (if (equal? dim (unit-dims rhs-units))
                            `(reduceindex y_out ,(cdr yindex) ,expr)
                            (error 'reduce-eq "dimension mismatch in rhs" 
                                   (variable-name y)
                                   (variable-label y)
                                   dim (unit-dims rhs-units)))
                        ))
                    ))
              )
             ((external? y)
              (let ((yindex (env-lookup (external-name y) extindexmap))
                    (dim (external-dim y)))
                (if (not yindex)
                    (error 'reduce-eq "variable not in external index" y)
                    (let ((expr (reduce-expr (make-call op `(getindex yext ,(cdr yindex)) rhs) indexmaps)))
                      (let ((rhs-units (expr-units rhs unit-env)))
                        (trace 'reduce-eq "dim = ~A dims(rhs-units) = ~A ~%" dim (unit-dims rhs-units))
                        (if (equal? dim (unit-dims rhs-units))
                            `(reduceindex yext_out ,(cdr yindex) ,expr)
                            (error 'reduce-eq "dimension mismatch in external rhs" 
                                   (external-name y)
                                   (external-label y)
                                   dim (unit-dims rhs-units)))
                        ))
                    ))
              )
             (else
              (error 'reduce-eq "unsupported reduce variable type" y))
             ))

           
           (($ salt#extevlink src dst)
            (let ((src-index (env-lookup (externalev-name src) extevindexmap))
                  (dst-index (env-lookup (externalev-name dst) extevindexmap)))
                (if (not src-index)
                    (error 'reduce-eq "event not in external event index" src))
                (if (not dst-index)
                    (error 'reduce-eq "event not in external event index" dst))
                `(setindex ext_out ,(cdr dst-index)
                           (signal.reinitPriority (getindex extev ,(cdr src-index))
                                                  (getindex extev ,(cdr src-index))
                                                  ))
                ))
           
           (($ salt#evcondition name rhs)
            (let ((evindex (env-lookup name evindexmap))
                  (expr (reduce-expr rhs indexmaps)))
              (if (not evindex)
                  (error 'reduce-eq "event identifier not found in indexmap" name))
              (trace 'reduce-eq "evcondition name = ~A~%" name)
              (trace 'reduce-eq "expr = ~A~%" expr)
              (if (expr-units rhs unit-env)
                  `(setindex c_out ,(cdr evindex) ,expr)
                  (error 'reduce-eq "dimension mismatch in condition" expr))
              ))

           
           (($ salt#evresponse name rhs)
            (trace 'reduce-eq "evresponse name = ~A rhs = ~A~%" name rhs)
            (let ((y (match rhs
                            (('signal.reinit e ($ salt#left-var u) yindex . rest) u)
                            (error 'reduce-eq "unknown event response equation" eq))))
              (cond
               ((variable? y)
                (let (
                      (yindex (env-lookup (variable-name y) cindexmap))
                      (dim (variable-dim y))
                      )
                  (if (not yindex)
                      (error 'reduce-eq "variable not in index" y)
                      (let ((expr (reduce-expr rhs indexmaps)))
                        (if (equal? dim (unit-dims (expr-units rhs unit-env)))
                            `(setindex y_out ,(cdr yindex) ,expr)
                            (error 'reduce-eq "variable dimension mismatch in event response" 
                                   y dim (unit-dims (expr-units rhs unit-env))))
                        ))
                  ))
               ((discrete-variable? y)
                (let (
                      (yindex (env-lookup (discrete-variable-name y) dindexmap))
                      (dim (discrete-variable-dim y))
                      )
                  (if (not yindex)
                      (error 'reduce-eq "variable not in index" y)
                      (let ((expr (reduce-expr rhs indexmaps)))
                        (if (equal? dim (unit-dims (expr-units rhs unit-env)))
                            `(setindex d_out ,(cdr yindex) ,expr)
                            (error 'reduce-eq "variable dimension mismatch in event response" y))
                        ))
                  ))
               ((regime-variable? y)
                (let (
                      (rindex (env-lookup (regime-variable-name y) rindexmap))
                      )
                  (if (not rindex)
                      (error 'reduce-eq "regime variable not in index" y)
                      (let ((expr (reduce-expr rhs indexmaps)))
                        `(setindex r_out ,(car (cdr rindex)) ,expr)
                        ))
                  ))
               (else (error 'reduce-eq "unknown variable type in reinit equation" eq))
               ))
            )

           (else eq)))
  )          


(define (simcreate eqset)

  (define (next-index index iset)
    (if (member index iset)
        (next-index (+ 1 index) iset)
        index))
  
  (let* ((nodemap    (equation-set-nodemap eqset))
         (regimemap  (equation-set-regimemap eqset))
         (dimenv     (equation-set-dimenv eqset))
         (conditions (equation-set-conditions eqset))

         (nodelst   (reverse (env->list nodemap)))
         (unit-env
          (fold (lambda (p env)
                  (let* ((name (car p)) (rhs (cdr p))
                         (rhs-units (expr-units rhs env)))
                    (if (equal? (quantity-int (binding-value (env-lookup name dimenv)))
                                (quantity-int (unit-dims rhs-units)))
                        (extend-env-with-binding
                         env (gen-binding name rhs-units))
                        (error 'simcreate "dimension mismatch in initial value" 
                               name (binding-value (env-lookup name dimenv))
                               rhs rhs-units))
                    ))
                empty-env 
                (append
                 (equation-set-parameters eqset)
                 (equation-set-fields eqset)
                 (equation-set-externals eqset)
                 (equation-set-externalevs eqset))))

         (cindexmap 
           (let recur ((nodelst nodelst)
                       (indexmap  empty-env)
                       (index     0))
             (if (null? nodelst)
                 indexmap
                 (let ((node (car nodelst)))
                   (if (variable? (cdr node))
                       (recur (cdr nodelst)
                              (extend-env-with-binding indexmap (gen-binding (car node) index))
                              (+ 1 index))
                       (recur (cdr nodelst) indexmap index))))
             ))
          
         (dindexmap 
           (let recur ((nodelst nodelst)
                       (indexmap  empty-env)
                       (index     0))
             (if (null? nodelst)
                 indexmap
                 (let ((node (car nodelst)))
                   (if (discrete-variable? (cdr node))
                       (recur (cdr nodelst)
                              (extend-env-with-binding indexmap (gen-binding (car node) index))
                              (+ 1 index))
                       (recur (cdr nodelst) indexmap index))))
             ))
          

         (pindexmap 
          (let recur ((nodelst nodelst)
                      (indexmap  empty-env)
                      (index     0))
             (if (null? nodelst)
                 indexmap
                 (let ((node (car nodelst)))
                   (cond ((parameter? (cdr node))
                          (recur (cdr nodelst)
                                 (extend-env-with-binding indexmap (gen-binding (car node) index))
                                 (+ 1 index)))
                         
                         (else
                          (recur (cdr nodelst) indexmap index)))))
             ))
          
         (fldindexmap 
          (let recur ((nodelst nodelst)
                      (indexmap  empty-env)
                      (index     0))
            (if (null? nodelst)
                indexmap
                (let ((node (car nodelst)))
                  (cond ((field? (cdr node))
                         (recur (cdr nodelst)
                                (extend-env-with-binding indexmap (gen-binding (car node) index))
                                (+ 1 index)))
                         
                         (else
                          (recur (cdr nodelst) indexmap index)))))
             ))
          
         (evindexmap 
           (let recur ((conditions conditions)
                       (indexmap  empty-env)
                       (index     0))
             (if (null? conditions)
                 indexmap
                 (recur (cdr conditions)
                        (extend-env-with-binding indexmap (gen-binding (evcondition-name (car conditions)) index))
                        (+ 1 index)))
             ))

         (extindexmap 
          (let ((initial-env+iset
                 (fold (lambda (node env+iset)
                         (cond ((external? (cdr node))
                                (let ((env (car env+iset))
                                      (iset (cdr env+iset))
                                      (ord (external-order (cdr node))))
                                  (if ord
                                      (if (member ord iset)
                                          (error 'simcreate "the same order of external quantity is specified more than once"
                                                 (car node) ord)
                                          (cons
                                           (extend-env-with-binding env (gen-binding (car node) ord))
                                           (lset-union eq? iset (list ord)))
                                          )
                                      env+iset)))
                               (else env+iset)))
                       (cons empty-env '())
                       nodelst)))
            
            (let ((initial-env (car initial-env+iset))
                  (iset (cdr initial-env+iset)))
              
              (let recur ((nodelst nodelst)
                          (indexmap  initial-env)
                          (index     0))
                (if (null? nodelst)
                    indexmap
                    (let ((node (car nodelst)))
                      (cond ((and (external? (cdr node))
                                  (not (env-lookup (car node) indexmap)))
                             (let ((index1 (next-index index iset)))
                               (recur (cdr nodelst)
                                      (extend-env-with-binding indexmap (gen-binding (car node) index1))
                                      (+ 1 index1))))
                            (else
                             (recur (cdr nodelst) indexmap index)))))
                ))
            ))

          
         (extevindexmap
          (let ((initial-env+iset
                 (fold (lambda (node env+iset)
                         (cond ((externalev? (cdr node))
                                (let ((env (car env+iset))
                                      (iset (cdr env+iset))
                                      (ord (externalev-order (cdr node))))
                                  (if ord
                                      (if (member ord iset)
                                          (error 'simcreate "the same order of external event is specified more than once"
                                                 (car node) ord)
                                          (cons
                                           (extend-env-with-binding env (gen-binding (car node) ord))
                                           (lset-union eq? iset (list ord)))
                                          )
                                      env+iset)))
                               (else env+iset)))
                       (cons empty-env '())
                       nodelst)))
            
            (let ((initial-env (car initial-env+iset))
                  (iset (cdr initial-env+iset)))
              
              (let recur ((nodelst nodelst)
                          (indexmap  initial-env)
                          (index     0))
                (if (null? nodelst)
                    indexmap
                    (let ((node (car nodelst)))
                      (cond ((and (externalev? (cdr node))
                                  (not (env-lookup (car node) indexmap)))
                             (let ((index1 (next-index index iset)))
                               (recur (cdr nodelst)
                                      (extend-env-with-binding indexmap (gen-binding (car node) index1))
                                      (+ 1 index1))))
                            (else
                             (recur (cdr nodelst) indexmap index)))))
                ))
            ))
         
         (indexmaps `((cindexmap  . ,cindexmap)
                      (pindexmap  . ,pindexmap)
                      (fldindexmap  . ,fldindexmap)
                      (dindexmap  . ,dindexmap)
                      (rindexmap  . ,regimemap)
                      (evindexmap . ,evindexmap)
                      (extindexmap . ,extindexmap)
                      (extevindexmap . ,extevindexmap)
                      ))
         
         (param-block 
          (map (lambda (x) (reduce-constant-expr (cdr x) indexmaps))
               (equation-set-parameters eqset)))

         (field-block 
          (map (lambda (x) (reduce-constant-expr (cdr x) indexmaps))
               (equation-set-fields eqset)))

         (init-block 
          (map (lambda (x) (reduce-expr (cdr x) indexmaps))
               (equation-set-definitions eqset)))

         (discrete-init-block
          (map (lambda (x) (reduce-expr (cdr x) indexmaps))
               (equation-set-discrete-definitions eqset)))
         
         (external-init-block
          (map (lambda (x) (reduce-expr (cdr x) indexmaps))
               (equation-set-externals eqset)))
         
         (externalev-init-block
          (map (lambda (x) (reduce-expr (cdr x) indexmaps))
               (equation-set-externalevs eqset)))

         (externalev-link-block
          (filter-map (lambda (node)
                        (and (externalev? (cdr node))
                             (let ((link (externalev-link (cdr node))))
                               (and link
                                    (reduce-eq link indexmaps unit-env)))))
                      nodelst))

         (eq-block
          (map (lambda (x) (reduce-eq x indexmaps unit-env)) 
               (equation-set-equations eqset)))
         
         (cond-block
          (map (lambda (c) (reduce-eq c indexmaps unit-env))
               (equation-set-conditions eqset)))

         (pos-responses 
          (map (lambda (x) (reduce-eq x indexmaps unit-env))
               (equation-set-pos-responses eqset)))

         (neg-responses 
          (map (lambda (x) (reduce-eq x indexmaps unit-env))
               (equation-set-neg-responses eqset)))

         )

    (make-simruntime 
     eqset
     cindexmap dindexmap 
     evindexmap regimemap
     extindexmap extevindexmap
     param-block
     field-block
     init-block
     discrete-init-block
     external-init-block
     externalev-init-block
     externalev-link-block
     eq-block 
     cond-block
     pos-responses
     neg-responses)

    )
  )



        

)