;;       
;; An extension for translating NEMO models to NEST code.
;;
;; Copyright 2011-2013 Ivan Raikov and the Okinawa Institute of Science and Technology
;;
;; This program is free software: you can redistribute it and/or
;; modify it under the terms of the GNU General Public License as
;; published by the Free Software Foundation, either version 3 of the
;; License, or (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful, but
;; WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;; General Public License for more details.
;;
;; A full copy of the GPL license can be found at
;; <http://www.gnu.org/licenses/>.
;;

(module nemo-nest

	(nemo:nest-translator)

	(import scheme chicken utils data-structures ports extras srfi-1 srfi-13 srfi-69)
	
	(require-extension lolevel posix matchable strictly-pretty 
			   varsubst datatype nemo-core nemo-utils nemo-units
                           nemo-geometry nemo-defaults nemo-constraints
			   nemo-gate-complex nemo-synapse nemo-currents

                           )


(define C++-ops
  `(+ - * / > < <= >= =))



(define (nest-name s)
  (let ((cs (string->list (->string s))))
    (let loop ((lst (list)) (cs cs))
      (if (null? cs) (string->symbol (list->string (reverse lst)))
	  (let* ((c (car cs))
		 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
			   (else #\_))))
	    (loop (cons c1 lst) (cdr cs)))))))


(define builtin-consts
  (append `(params)
          (map (lambda (x) (nest-name (s+ "M_" (first x)))) nemo:math-constants)))


(define builtin-fns
  `(+ - * / pow neg abs atan asin acos sin cos exp ln
      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
      > < <= >= = and or round ceiling floor max min
      ))

	    
(define (rhsexpr/C++ expr)
  (match expr 
	 (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/C++ x)) es)))
	 (('pow x y)  (if (and (integer? y)  (positive? y))
			  (if (> y 1)  (let ((tmp (gensym "x")))
					 `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
			      x)
			    (if (and (number? y) (zero? y)) 1.0 expr)))
	 ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (nest-name s))
					     (map (lambda (x) (rhsexpr/C++ x)) es)) expr))
	 (id          (if (symbol? id) 
                          (if (assoc id nemo:math-constants)
                              (nest-name (s+ "M_" id))
                              (nest-name id))
                          id))
         ))


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


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


(define group/C++   (doc:block 2 (doc:text "(") (doc:text ")")))
(define block/C++   (doc:block 2 (doc:text "{") (doc:text "}")))
(define (stmt/C++ x) 
  (match x
	 (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
	 (else  (doc:cons x (doc:text ";")))))


(define (ifthen/C++ c e1 e2)
  (doc:nest 2
    (doc:connect (doc:group (doc:connect (doc:text "if") c))
		 (doc:connect (doc:nest 2 e1)
			      (doc:nest 2 (doc:connect 
					   (doc:text "else") 
					   e2))))
    ))


(define (letblk/C++ e1 e2)
  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
	((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
	(else (doc:connect (doc:group (doc:nest 2 (stmt/C++ e1)))
			   (doc:group (doc:nest 2 e2))))))
	

(define (format-op/C++ indent op args)
  (let ((op1 (doc:text (->string op))))
    (if (null? args) op1
	(match args
	       ((x)      (doc:concat (list op1 x)))
	       ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
	       ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
	       (lst      (let* ((n   (length lst))
				(n/2 (inexact->exact (round (/ n 2)))))
			   (doc:concat 
			    (intersperse 
			     (list (format-op/C++ indent op (take lst n/2 )) op1 
				   (format-op/C++ indent op (drop lst n/2 )))
			     (doc:space)))))))))


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


(define (name-normalize expr)
  (match expr 
	 (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
	 (('let bs e)
	  `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
	 ((f . es) 
	  (cons (if (member f builtin-fns) f (nest-name f)) (map name-normalize es)))
	 ((? symbol? ) (nest-name expr))
	 ((? atom? ) expr)))


(define (canonicalize-expr/C++ expr)
  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term)))
    (let* ((expr1 (if-convert expr))
	   (expr2 (subst-convert expr1 subst-empty))
	   (expr3 (let-lift expr2))
	   (expr4 (name-normalize expr3)))
      expr4)))


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

       ((op . rest)  
       (let ((op (case op ((ln) 'log) ((abs) 'fabs) (else op))))
	 (let ((fe
		(if (member op C++-ops)
		    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
			  (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
			  (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
		      (case op
			((/)  
			 (format-op/C++ indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/C++ indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if (or mul? plmin?) (group/C++ fx) fx)))) rest)))
			((*)  
			 (format-op/C++ indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/C++ indent+ x)))
						   (if (or (symbol? x) (number? x)) fx
						       (if plmin? (group/C++ fx) fx)))) rest)))
			
			(else
			 (format-op/C++ indent op 
					  (map (lambda (x) 
						 (let ((fx (format-expr/C++ indent+ x))) fx)) rest)))))
		    
		    (let ((op (case op ((neg) '-) (else op))))
		      (format-fncall/C++ indent op (map (lambda (x) (format-expr/C++ indent+ x)) rest))))))
	   (if rv 
	       (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
	       fe))))
      
      (else  (let ((fe (doc:text (->string expr))))
	       (if rv 
		   (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
		   fe)))))))
	       
	  
(define (expr->string/C++ x . rest)
  (let-optionals rest ((rv #f) (width 72))
    (sdoc->string (doc:format width (format-expr/C++ 2 x rv)))))


(define (make-define-fn sysname )
  (lambda (indent n proc)

    (let ((lst (procedure-data proc))
	  (indent+ (+ 2 indent)))

      (let ((rt       (or (lookup-def 'rt lst) 'double))
	    (formals  (lookup-def 'formals lst))
	    (vars     (lookup-def 'vars lst))
	    (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
	    (body     (lookup-def 'body lst))
	    (rv       (gensym 'rv)))

	(let ((argument-list 
	       (append
		(if (null? vars) '() (map (lambda (x) (sprintf "double ~A" (nest-name x))) vars))
		'("const void* params"))))
	  (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) "{" ))
	  (let* ((body0 (rhsexpr/C++ body))
		 (body1 (canonicalize-expr/C++ (add-params-to-fncall body0 builtin-fns)))
		 (lbs   (enum-bnds body1 (list))))
	    (pp indent+ (double ,rv ";"))
	    (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs) ";")))
	    (if (not (null? consts)) 
		(begin (pp indent+ (double ,(slp ", " (map (compose nest-name car) consts)) ";")
			   (,(sprintf "const ~A::Parameters_" sysname) "& p =  " 
                            ,(sprintf " *(reinterpret_cast< const ~A::Parameters_ *>(params));" sysname)))
		       (for-each (lambda (x) (let ((n (car x))) 
                                               (pp indent+ (,(nest-name n) = ,(sprintf "p.~A;" (nest-name n))))))
                                 consts)
		       ))
	    (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
	    (pp indent+ ,(s+ "return " rv ";"))
	    ))
	(pp indent "}"))
    )))


(define (make-define-fn-header sysname )
  (lambda (indent n proc)

    (let ((lst (procedure-data proc))
	  (indent+ (+ 2 indent)))

      (let ((rt       (or (lookup-def 'rt lst) 'double))
	    (formals  (lookup-def 'formals lst))
	    (vars     (lookup-def 'vars lst))
	    (consts   (filter (lambda (x) (not (procedure? (cdr x)))) (lookup-def 'consts lst)))
	    (body     (lookup-def 'body lst))
	    (rv       (gensym 'rv)))

	(let ((argument-list 
	       (append
		(if (null? vars) '() (map (lambda (x) (sprintf "double ~A" (nest-name x))) vars))
		'("const void* params"))))
	  (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " argument-list)) ";" ))
	  ))
      )))



(define (ith v i) (sprintf "Ith(~A,~A)" v i))


(define (output-prelude sysname steady-state-index-map indent)

  (pp indent ,#<#EOF
##include "#{sysname}.h"
##include "exceptions.h"
##include "network.h"
##include "dict.h"
##include "integerdatum.h"
##include "doubledatum.h"
##include "dictutils.h"
##include "numerics.h"
##include <limits>

##include "universal_data_logger_impl.h"

##include <iomanip>
##include <iostream>
##include <cstdio>
##include <cstring>
##include <cmath>
EOF
))




(define cvode-prelude-header
#<<EOF

  /**
   * Exception to be thrown if a CVode solver does not return CV_SUCCESS
   * @ingroup KernelExceptions
   */
  class CVodeSolverFailure: public KernelException
  {
  public:
    /**
    * @note model should be passed from get_name() to ensure that
    *             names of copied models are reported correctly. 
    * @param model name of model causing problem
    * @param status exit status of the CVode solver
    */
    CVodeSolverFailure(const std::string& model,
                     const int status)
      : KernelException("CVodeSolverFailure"),
        model_(model),
        status_(status)
      {}
    ~CVodeSolverFailure() throw() {}
    
    std::string message()
    {
      std::ostringstream msg;
      msg << "In model " << model_ << ", the CVODE solver "
          << "returned with exit status " << status_ << ".\n";
      return msg.str();
    }

    private:
      const std::string model_;
      const int status_;
  };


EOF
)


(define ida-prelude-header
#<<EOF

  /**
   * Exception to be thrown if an IDA solver does not return IDA_SUCCESS
   * @ingroup KernelExceptions
   */
  class IDASolverFailure: public KernelException
  {
  public:
    /**
    * @note model should be passed from get_name() to ensure that
    *             names of copied models are reported correctly. 
    * @param model name of model causing problem
    * @param status exit status of the CVode solver
    */
    IDASolverFailure(const std::string& model,
                     const int status)
      : KernelException("IDASolverFailure"),
        model_(model),
        status_(status)
      {}
    ~IDASolverFailure() throw() {}
    
    std::string message()
    {
      std::ostringstream msg;
      msg << "In model " << model_ << ", the IDA solver "
          << "returned with exit status " << status_ << ".\n";
      return msg.str();
    }

    private:
      const std::string model_;
      const int status_;
  };


EOF
)


(define sundials-prelude 
#<<EOF

  /* 
   * Check function return value.
   *    opt == 0 means SUNDIALS function allocates memory so check if
   *             returned NULL pointer
   *    opt == 1 means SUNDIALS function returns a flag so check if
   *             flag >= 0
   *    opt == 2 means function allocates memory so check if returned
   *             NULL pointer 
   */

  static int check_flag(void *flagvalue, const char *funcname, int opt)
  {
    int *errflag;

    /* Check if CVode function returned NULL pointer - no memory allocated */
    if (opt == 0 && flagvalue == NULL) 
    {
      fprintf(stderr, 
              "\nCVode error: %s() failed - returned NULL pointer\n\n",
              funcname);
      return(1); 
    }

    /* Check if flag < 0 */
    else if (opt == 1) 
    {
      errflag = (int *) flagvalue;
      if (*errflag < 0) {
        fprintf(stderr, 
                "\nCVode error: %s() failed with flag = %d\n\n",
                funcname, *errflag);
        return(1); }
    }

  /* Check if function returned NULL pointer - no memory allocated */
  else if (opt == 2 && flagvalue == NULL) 
  {
    fprintf(stderr, 
            "\nMemory error: %s() failed - returned NULL pointer\n\n",
            funcname);
    return(1); 
  }

    return(0);
  }


  void adjust_zero_crossings (N_Vector v, double abstol)
  {
    int i;
    for (i = 0; i < NV_LENGTH_S(v); i++)
    {
        if (fabs(NV_Ith_S(v,i)) < abstol) NV_Ith_S(v,i) = 0;
    }
    return;
  }

EOF
)


(define (output-event-handle
	 sysname state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info
	 indent indent+)

  (let ((isyns  (lookup-def 'i-synapses synapse-info))
        (pscs  (lookup-def 'post-synaptic-conductances synapse-info)))
  
  (pp indent  ,nl (,#<#EOF
void #{sysname}::handle(SpikeEvent & e)
  {
    int flag;
    assert(e.get_delay() > 0);
    flag = 0;
EOF
   ))

  (for-each
   (lambda (isyn psc)
     (let ((wscale (fourth isyn))
           (wthreshold (fifth isyn)))

       (pp indent+ (,#<#EOF
    if ((flag==0) && (e.get_rport() == #{(nest-name (first psc))}_SPIKE_RECEPTOR ) && (e.get_weight() > #(or (and wthreshold (sprintf "P_.~A" (nest-name wthreshold))) 0.0)))
      {
	B_.spike_#(nest-name (second psc)).add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
			                             fabs(e.get_weight()) * e.get_multiplicity());
        flag = 1;
      }
EOF
       )) 
   ))
   isyns pscs)

  (pp indent ,nl #\} ,nl)

  (pp indent (,#<#EOF
void #{sysname}::handle(CurrentEvent& e)
  {
    assert(e.get_delay() > 0);

    const double_t c = e.get_current();
    const double_t w = e.get_weight();

    B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), w * c);
  }

void #{sysname}::handle(DataLoggingRequest& e)
  {
    B_.logger_.handle(e);
  }
EOF
  ))

))

(define (output-header 
	 sysname method state-index-map steady-state-index-map 
	 defaults external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info
	 indent indent+)

  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))

    (pp indent
        ("#include \"nest.h\"")
        ("#include \"event.h\"")
        ("#include \"archiving_node.h\"")
        ("#include \"ring_buffer.h\"")
        ("#include \"connection.h\"")
        ("#include \"universal_data_logger.h\"")
        ("#include \"recordables_map.h\"")
        )

    (case method
      ((leapfrog) 
       (pp indent
           ("#define Ith(v,i)    (v[i])   /* Ith component in a vector */")
           ))
      ((cvode) 
       (pp indent
           ("#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */")
           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
           ("#include <cvode/cvode_diag.h>        /* prototype for CVDiag */")
           ("#include <sundials/sundials_types.h> /* definition of type realtype */")
           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
           ))
      ((ida) 
       (pp indent
           ("#include <ida/ida.h>                 /* prototypes for IDA fcts., consts. */")
           ("#include <ida/ida_dense.h>")
           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
           ))
      ((gsl) 
       (pp indent
           ("#include <gsl/gsl_errno.h>")
           ("#include <gsl/gsl_matrix.h>")
           ("#include <gsl/gsl_sf_exp.h>")
           ("#include <gsl/gsl_odeiv2.h>")
           ("#define Ith(v,i)    (v[i])")
           ))
      )

    (pp indent ,nl)

    (pp indent "namespace nest {" ,nl)
  
    (case method
      ((cvode) 
       (pp indent ,cvode-prelude-header))
      ((ida) 
       (pp indent ,ida-prelude-header)))
      
  
    (case method
      ((leapfrog) 
       (pp indent 
           ("typedef int (*rhsfn_t)(double t, const double *y, double *ydot, void *user_data);") ,nl
           (extern "\"C\""
                   int ,(sprintf "~A_dynamics" sysname)
                   "(double, const double *, double *, void*)" #\;) ,nl))
      ((cvode) 
       (pp indent (extern "\"C\""
                          int ,(sprintf "~A_dynamics" sysname)
                          "(double, const N_Vector, N_Vector, void*)" #\;) ,nl
                   (extern "\"C\""
                           int ,(sprintf "~A_event" sysname)
                           "(double, N_Vector, double *, void*)" #\;) ,nl
                           ))
      ((ida) 
       (pp indent (extern "\"C\""
                          int ,(sprintf "~A_residual" sysname)
                          "(double, N_Vector, N_Vector, N_Vector, void*)" #\;) ,nl
                   (extern "\"C\""
                           int ,(sprintf "~A_event" sysname)
                           "(double, N_Vector, N_Vector, double *, void*)" #\;) ,nl
                           ))
      (else
       (pp indent (extern "\"C\""
                          int ,(sprintf  "~A_dynamics" sysname)
                          "(double, const double*, double*, void*)" #\;) ,nl))
      )

    (pp indent 
        (,(sprintf "class ~A:" sysname) "public Archiving_Node { "))

    (pp indent+ ( ,#<#EOF
public: 
 ~#{sysname} ();
 #{sysname} (const #{sysname} &);
 #{sysname} ();
EOF
    ))

    (pp indent ( #<<EOF
    using Node::connect_sender;
    using Node::handle;

    port check_connection(Connection&, port);
    
    void handle(SpikeEvent &);
    void handle(CurrentEvent &);
    void handle(DataLoggingRequest &); 
    
    port connect_sender(SpikeEvent &, port);
    port connect_sender(CurrentEvent &, port);
    port connect_sender(DataLoggingRequest &, port);
    
    void get_status(DictionaryDatum &) const;
    void set_status(const DictionaryDatum &);
    
    void init_node_(const Node& proto);
    void init_state_(const Node& proto);
    void init_buffers_();
    void calibrate();
    
    void update(Time const &, const long_t, const long_t);

EOF
))

  (pp indent (,#<#EOF

    /**
     * Minimal spike receptor type.
     * @note Start with 1 so we can forbid port 0 to avoid accidental
     *       creation of connections with no receptor type set.
     */
    static const port MIN_SPIKE_RECEPTOR = 1;

    /** 
     * Spike receptors.
     */
    enum SpikeSynapseTypes { #{(slp ", " (append
                                         (if (pair? pscs)
                                             (cons (sprintf "~A_SPIKE_RECEPTOR=MIN_SPIKE_RECEPTOR" (first (car pscs)))
                                                   (map (lambda (x) (sprintf "~A_SPIKE_RECEPTOR" (first x))) (cdr pscs))) '())
                                         (list "SUP_SPIKE_RECEPTOR ")))}
                             };


EOF
))
  
  (case method
    ((cvode) 
     (pp indent+ (friend int ,(sprintf "~A_dynamics" sysname)
                         "(double, const N_Vector, N_Vector, void*)" #\; ) ,nl
                 (friend int ,(sprintf "~A_event" sysname)
                         "(double, N_Vector, double *, void*)" #\;) ,nl
                         ))
    ((ida) 
     (pp indent+ (friend int ,(sprintf "~A_residual" sysname)
                         "(double, N_Vector, N_Vector, N_Vector, void*)" #\; ) ,nl
                 (friend int ,(sprintf "~A_event" sysname)
                         "(double, N_Vector, N_Vector, double *, void*)" #\;) ,nl
                         ))
    (else
     (pp indent+ (friend int ,(sprintf "~A_dynamics" sysname)
                         "(double, const double*, double*, void*)" #\;) ,nl))
    )

  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
    (for-each
     (lambda (psc)
       (pp indent+ (,(sprintf "int ~A_transients (long_t lag);" (first psc)))))
     pscs))
  

  (pp indent (,#<#EOF
    // The next two classes need to be friends to access the State_ class/member
    friend class RecordablesMap<#{sysname}>;
    friend class UniversalDataLogger<#{sysname}>;
EOF
   ))

  
  (pp indent+ ,nl "struct Parameters_ { ")

  (for-each 
   (lambda (x) (pp indent+ (,(sprintf "double ~A;" x))))
   (append 
    (map (compose nest-name first) const-defs)
    (map (compose nest-name first) defaults)))

  (pp indent+ 
      ("Parameters_();")
      ("void get(DictionaryDatum&) const;")
      ("void set(const DictionaryDatum&);"))

  (pp indent+ "};")

  (pp indent+ ,nl "struct State_ { ")

  (pp indent+ ,nl 
      "enum StateVecElems {"
      (,(slp ", " (map (lambda (x) 
			 (let ((n (string-upcase (->string (first x))))
			       (ni (second x)))
			   (sprintf "~A = ~A" n ni)))
		       state-index-map)))
      "};")

  (pp indent+ ,nl
      (,(sprintf "double y_[~A];" (length state-index-map)))

      "State_(const Parameters_& p);" 
      "State_(const State_& s);"
      "State_& operator=(const State_& s);"
      "void get(DictionaryDatum&) const;"
      "void set(const DictionaryDatum&, const Parameters_&);")

  (case method
    ((gsl leapfrog) 
     (pp indent+ ,nl "int_t r_; /* refractory counter */")))

  (pp indent+ "};")

  (case method
    ((cvode ida) 
     (pp indent+ ,nl "struct Variables_ {};"))
    (else
     (pp indent+ ,nl "struct Variables_ { int_t RefractoryCounts_; double U_old_; /* for spike-detection */ };"))
    )

  (pp indent+ ,nl "struct Buffers_ { ")

  (pp indent+ ,nl 
      (,(sprintf "Buffers_(~A&);" sysname))
      (,(sprintf "Buffers_(const Buffers_&, ~A&);" sysname))
      (,(sprintf "UniversalDataLogger<~A> logger_;" sysname)))

  (pp indent+ ,nl

      (,(case method

    ((leapfrog)
#<<EOF
  rhsfn_t sys_; // pointer to function that computes dv/dt
  double *u, *v, *dvdt;  // intermediate state vectors used by leapfrog method
  unsigned int  N;  // size of state vector used by leapfrog method
EOF
)

    ((cvode)
#<<EOF
  N_Vector y;    //!< current state vector used by CVode
  void *   sys_;  //!< CVode control structure
EOF
)
    ((ida)
#<<EOF
  N_Vector y;    //!< current state vector used by IDA
  N_Vector yp;    //!< derivatives vector used by IDA
  void *   sys_;  //!< IDA control structure
EOF
)

    ((gsl)
#<<EOF
  gsl_odeiv2_step*    s_;    //!< stepping function
  gsl_odeiv2_control* c_;    //!< adaptive stepsize control function
  gsl_odeiv2_evolve*  e_;    //!< evolution function
  gsl_odeiv2_system   sys_;  //!< struct describing system
  unsigned int N;  // size of state vector used by Jacobian
  double *u, *jac;  // intermediate state vectors used for Jacobian approximation
EOF
)

    )))


  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
    (for-each
     (lambda (psc)
       (pp indent+ (,(sprintf "RingBuffer spike_~A;" (second psc)))))
     pscs))

  (pp indent+ ,nl
#<<EOF
  RingBuffer currents_;

  double_t step_;           //!< step size in ms
  double   IntegrationStep_;//!< current integration time step, updated by solver

  /** 
   * Input current injected by CurrentEvent.
   * This variable is used to transport the current applied into the
   * _dynamics function computing the derivative of the state vector.
   * It must be a part of Buffers_, since it is initialized once before
   * the first simulation, but not modified before later Simulate calls.
   */
  double_t I_stim_;
EOF
)
  (pp indent+ "};")

  (pp indent+
      "template <State_::StateVecElems elem>"
      "double_t get_y_elem_() const { return S_.y_[elem]; }"
      "Parameters_ P_;"
      "State_      S_;"
      "Variables_  V_;"
      "Buffers_    B_;"
      (,(sprintf "static RecordablesMap<~A> recordablesMap_;" sysname))
      "}; ")

  (pp indent+ (,#<#EOF
  inline
  port #{sysname}::check_connection(Connection& c, port receptor_type)
  {
    SpikeEvent e;
    e.set_sender(*this);
    c.check_event(e);
    return c.get_target()->connect_sender(e, receptor_type);
  }

  inline
  port #{sysname}::connect_sender(SpikeEvent&, port receptor_type)
  {
    if ( receptor_type < MIN_SPIKE_RECEPTOR || receptor_type >= SUP_SPIKE_RECEPTOR )
    {
      if ( receptor_type < 0 || receptor_type >= SUP_SPIKE_RECEPTOR )
	throw UnknownReceptorType(receptor_type, get_name());
      else
	throw IncompatibleReceptorType(receptor_type, get_name(), "SpikeEvent");
    }
    return receptor_type;
  }
 
  inline
  port #{sysname}::connect_sender(CurrentEvent&, port receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return 0;
  }

  inline
  port #{sysname}::connect_sender(DataLoggingRequest& dlr, 
				      port receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
  }

  inline
    void #{sysname}::get_status(DictionaryDatum &d) const
  {
    P_.get(d);
    S_.get(d);
    Archiving_Node::get_status(d);

    (*d)[names::recordables] = recordablesMap_.get_list();

    def<double_t>(d, names::t_spike, get_spiketime_ms());

    DictionaryDatum receptor_dict_ = new Dictionary();

    #{(string-concatenate
       (map (lambda (psc) 
              (sprintf "    (*receptor_dict_)[Name(\"~A\")]  = ~A_SPIKE_RECEPTOR;~%"
                       (first psc) (first psc)))
            pscs))}
               
    (*d)[names::receptor_types] = receptor_dict_;
  }

  inline
    void #{sysname}::set_status(const DictionaryDatum &d)
  {
    Parameters_ ptmp = P_;  // temporary copy in case of errors
    ptmp.set(d);                       // throws if BadProperty
    State_      stmp = S_;  // temporary copy in case of errors
    stmp.set(d, ptmp);                 // throws if BadProperty

    // We now know that (ptmp, stmp) are consistent. We do not 
    // write them back to (P_, S_) before we are also sure that 
    // the properties to be set in the parent class are internally 
    // consistent.
    Archiving_Node::set_status(d);

    // if we get here, temporaries contain consistent set of properties
    P_ = ptmp;
    S_ = stmp;

    calibrate();
  }
EOF
  ))

  (pp indent "}" ,nl)
  )
)



(define (output-jac sysname method state-index-map pool-ions i-eqs v-eq indent indent+)

;; Diagonal Jacobian approximation: (f(s+.01) - f(s))/.001 

;; CVODE:
;;  typedef int (*CVDlsDenseJacFn)(long int N, realtype t,
;;                               N_Vector y, N_Vector fy, 
;;                               DlsMat Jac, void *user_data,
;;                               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);

;; GSL:
;; int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);



  (case method

    ((gsl)
     (pp indent ,nl
         (,(sprintf "extern \"C\" int ~A_jacobian" sysname)
          (,(slp ", " `("double t"
                        "const double y[]"
                        "double *dfdy"
                        "double dfdt[]"
                        "void* pnode"
                        ))))
          #\{

          ("// cast the node ptr to an object of the proper type")
          ("assert(pnode);")
          ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
          (,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")

          ("// state is a reference to the model state ")
          (struct ,(sprintf "~A::Buffers_" sysname) "*b;")
          ("b = &(vnode.B_);")
          (,#<#EOF
   for (int i = 0; i < b->N; i++) 
   {
       b->u[i] = y[i] + 0.01;
   }
   #{sysname}_dynamics(t, b->u, b->jac, pnode);
   for (int i = 0; i < b->N; i++) 
   {
       dfdt[i*b->N + i] = (b->jac[i] - dfdy[i]) / .001;
   }
   return 0;
EOF
)
   #\}
   ))
))


(define (output-cvode-event sysname method imports const-defs state-index-map 
	   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
	   defaults i-eqs v-eq indent indent+)

  (let ((vi (lookup-def 'v state-index-map)))
    
    (pp indent ,nl
        (,(sprintf "extern \"C\" int ~A_event" sysname)
         (,(slp ", " `("double t"
                       "N_Vector y"
                       "double *g"
                       "void* pnode"
                       )))
         #\{
         ))

      (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;")

      ,nl

      ("// S is shorthand for the type that describes the model state ")
      (typedef ,(sprintf "~A::State_" sysname) "S;")
           
      ,nl
           
      ("// cast the node ptr to an object of the proper type")
      ("assert(pnode);")
      ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
      
      ,nl
      
      ("// params is a reference to the model parameters ")
      (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
      ("params = &(node.P_);")
      
      ,nl
      
      ("// state is a reference to the model state ")
      (const struct ,(sprintf "~A::State_" sysname) "*state;")
      ("state = &(node.S_);")
      
      )
      
      (let ((vt (lookup-def 'V_t defaults)))
        (if (and vi vt)
            (pp indent+
                (,(expr->string/C++ (ith 'y vi) 'v))
                (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt))
                (,(expr->string/C++ '(- v vt) "g[0]"))
                ,nl
                ("return 0; "))
            (pp indent+
                (,(expr->string/C++ -1.0 "g[0]")))
            ))

      
      (pp indent (#\}))
      
      ))


(define (output-ida-event sysname method imports const-defs state-index-map 
	   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
	   defaults i-eqs v-eq indent indent+)

  (let ((vi (lookup-def 'v state-index-map)))
    
    (pp indent ,nl
        (,(sprintf "extern \"C\" int ~A_event" sysname)
         (,(slp ", " `("double t"
                       "N_Vector y"
                       "N_Vector yp"
                       "double *g"
                       "void* pnode"
                       )))
         #\{
         ))

      (pp indent+ ,nl ("double v, vt; v = -1.0; vt = 0.0;")

      ,nl

      ("// S is shorthand for the type that describes the model state ")
      (typedef ,(sprintf "~A::State_" sysname) "S;")
           
      ,nl
           
      ("// cast the node ptr to an object of the proper type")
      ("assert(pnode);")
      ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
      
      ,nl
      
      ("// params is a reference to the model parameters ")
      (const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
      ("params = &(node.P_);")
      
      ,nl
      
      ("// state is a reference to the model state ")
      (const struct ,(sprintf "~A::State_" sysname) "*state;")
      ("state = &(node.S_);")
      
      )
      
      (let ((vt (lookup-def 'V_t defaults)))
        (if (and vi vt)
            (pp indent+
                (,(expr->string/C++ (ith 'y vi) 'v))
                (,(expr->string/C++ (s+ "params->" (nest-name vt)) 'vt))
                (,(expr->string/C++ '(- v vt) "g[0]"))
                ,nl
                ("return 0; "))
            (pp indent+
                (,(expr->string/C++ -1.0 "g[0]")))
            ))

      
      (pp indent (#\}))
      
      ))




(define (cvode-solve spike-handle)
#<#EOF
        int N = NV_LENGTH_S (B_.y);
        tout = (current_steps+lag)*h;

	// adaptive step integration
	while (tt < tout)
	{
	  const int status = CVode (B_.sys_, tout, B_.y, &tt, CV_NORMAL);

          switch (status)
            {
                case CV_SUCCESS:      continue;
                case CV_ROOT_RETURN:  #{spike-handle};
                default:              throw CVodeSolverFailure (get_name(), 0);
            }
   	}

        for (int i = 0; i < N; i++)
        {
           S_.y_[i] = Ith(B_.y,i);
        }
EOF
)


(define (ida-solve spike-handle)
#<#EOF
        int N = NV_LENGTH_S (B_.y);
        tout = (current_steps+lag)*h;

	// adaptive step integration
	while (tt < tout)
	{
	  const int status = IDASolve (B_.sys_, tout, &tt, B_.y, B_.yp, IDA_NORMAL);

          switch (status)
            {
                case IDA_SUCCESS:      continue;
                case IDA_ROOT_RETURN:  #{spike-handle};
                case IDA_TSTOP_RETURN: break;
                default:               throw IDASolverFailure (get_name(), 0);
            }
   	}

        for (int i = 0; i < N; i++)
        {
           S_.y_[i] = Ith(B_.y,i);
        }
EOF
)


(define (gsl-solve vi spike-handle)
#<#EOF

   V_.U_old_ = S_.y_[#{vi}];

   while (tt < h)
   {

     const int status = gsl_odeiv2_evolve_apply
                         (B_.e_, B_.c_, B_.s_, 
			  &B_.sys_,              // system of ODE
			  &tt,                   // from t...
			  h,             // ...to t=t+h
			 &B_.IntegrationStep_ , // integration window (written on!)
			  S_.y_);	        // neuron state

	  if ( status != GSL_SUCCESS )
	    throw GSLSolverFailure(get_name(), status);

   }

   #{spike-handle}
EOF
)


(define (leapfrog-solve vi spike-handle)
#<#EOF
	V_.U_old_ = S_.y_[#{vi}];

        tt = (current_steps+lag)*h;
        tout = (current_steps+lag+1)*h;

	while (tt < tout)
	{
	   B_.sys_ (tt, S_.y_, B_.dvdt, reinterpret_cast<void*>(this));

	   for (int i = 0; i < B_.N; i++) 
           {
              B_.u[i] = S_.y_[i] + 0.5 * h * B_.dvdt[i];
           }

	   B_.sys_ (tt+0.5*h, B_.u, B_.dvdt, reinterpret_cast<void*>(this));

           for (int i = 0; i < B_.N; i++) 
           {  
              B_.v[i] = S_.y_[i] + h * B_.dvdt[i];
           }

	   tt += h;
	}

        for (int i = 0; i < B_.N; i++)
        {
           S_.y_[i] = B_.v[i];
        }
 
        #{spike-handle}
EOF
)



(define (output-synapse-transients sysname method state-index-map 
                                   const-defs transient-event-defs event-external-eq-defs 
                                   synapse-info imports constraints indent indent+)

  (let (
        (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
                                    (list (nest-name n) v)))
                      (lookup-def 'c-units constraints)))

        (isyns  (lookup-def 'i-synapses synapse-info))
        (pscs   (lookup-def 'post-synaptic-conductances synapse-info))
        (psc-transients (map (lambda (lst) (map nest-name lst)) 
                             (lookup-def 'psc-transients synapse-info)))
        (getstate (case method
                    ((cvode ida) (lambda (i) (sprintf "Ith(B_.y,~A)" i)))
                    (else (lambda (i) (sprintf "Ith(S_.y_,~A)" i)))))
        )


    (if (not (= (length event-external-eq-defs) (length pscs)))
        (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" 
               event-external-eq-defs pscs))
      
    (for-each (lambda (psc transients)

                (let* (
                       (ltransient-event-defs
                        (filter (lambda (x) (member (first x) transients))
                                transient-event-defs))

                       (levent-external-eq-def
                        (car
                         (fold (lambda (def ax)
                                 (let* ((n (nest-name (first def)) )
                                        (b (second def))
                                        (events (let ((fvs (enum-freevars b '() '())))
                                                  (filter (lambda (x) (member (first x) fvs)) event-external-eq-defs)))
                                        )
                                   (append events ax)))
                               '() ltransient-event-defs)))

                       (lconsts
                        (delete-duplicates
                         (fold (lambda (def ax)
                                 (let* ((n (nest-name (first def)) )
                                        (b (second def))
                                        (consts (let ((fvs (enum-freevars b '() '())))
                                                  (filter (lambda (x) (member (first x) fvs)) const-defs)))
                                       )
                                   (append consts ax)))
                               '() ltransient-event-defs)
                         (lambda (x y) (equal? (first x) (first y))))
                        ))
                  
                  (pp indent (,(sprintf "inline int ~A::~A_transients (long_t lag) {" sysname (first psc))))
                  
                  (if (not (null? lconsts)) 
                      (pp indent+ (double ,(slp ", " (map (compose nest-name car) lconsts)) ";")))

                  (let ((n (second levent-external-eq-def)))
                    (pp indent+
                        (,(sprintf "double_t ~A;" n))
                        (,(sprintf "~A = B_.spike_~A.get_value(lag);" 
                                   n (nest-name (second psc))))
                        )
                    )
                  
                  (pp indent+
                      (,(sprintf "if (~A > 0.0)" (second levent-external-eq-def)))
                      (#\{))
                     
                  (let* ((n      (nest-name (first levent-external-eq-def)))
                         (nu     (lookup-def n c-units))
                         (nscale (and nu (nemo:unit-scale nu)))
                         (b      (second levent-external-eq-def))
                         )
                    (pp (+ 2 indent+)
                        (,(sprintf "double_t ~A;" n))
                        (,(expr->string/C++ (if nscale `(* ,nscale ,b) b) n)))
                    )
                  
                  (for-each (lambda (def)
                              (let* ((n  (nest-name (first def)) )
                                     (ni (lookup-def n state-index-map))
                                     (b  (second def))
                                     (consts (let ((fvs (enum-freevars b '() '())))
                                               (filter (lambda (x) (member (first x) fvs)) const-defs)))
                                     )
                                
                              (if (not (null? consts)) 
                                  (for-each (lambda (x) 
                                              (let ((n (first x)))
                                                (pp (+ 2 indent+) (,(nest-name n) = ,(sprintf "P_.~A;" (nest-name n))))))
                                            consts)
                                  )
                              
                              (pp (+ 2 indent+ )
                                  (,(sprintf "double_t ~A;" n))
                                  ,(if ni (expr->string/C++ (getstate ni) n) '())
                                  ,(expr->string/C++ b n)
                                  ,(if ni (expr->string/C++ n  (getstate ni)) '())
                                  )
                              ))

                            ltransient-event-defs)
                  
                  
                  (pp (+ 2 indent+) "return 1;");
                  (pp indent+ (#\}))
                  (pp indent+ "return 0;");
                  (pp indent (#\}))
                  ))

               pscs psc-transients)

    ))


(define (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports indent)

  (let (
        (isyns  (lookup-def 'i-synapses synapse-info))
        (pscs   (lookup-def 'post-synaptic-conductances synapse-info))
        )

    (pp indent (,(sprintf "int events = 0;")))

    (if (not (= (length event-external-eq-defs) (length pscs)))
        (error 'nemo:nest-translator "mismatch between event variables and synaptic conductances" 
               event-external-eq-defs pscs))
      
    (for-each (lambda (psc)
                (pp indent
                    (,(sprintf "events = events + ~A_transients (lag);" (first psc))))
                )
              pscs)

    (case method
      ((cvode)
        (pp indent 
            ( "/* Reinitializes CVode state if any synaptic events have occurred */")
            ("if (events > 0)")
            (#\{)
            ("  int status = CVodeReInit (B_.sys_, tt, B_.y);")
            ("  if (check_flag(&status, \"CVodeReInit\", 1)) throw CVodeSolverFailure (get_name(), status);")
            (#\})
            ))
      ((ida)
        (pp indent 
            ( "/* Reinitializes IDA state if any synaptic events have occurred */")
            ("if (events > 0)")
            (#\{)
            ("  int status = IDAReInit (B_.sys_, tt, B_.y, B_.yp);")
            ("  if (check_flag(&status, \"IDAReInit\", 1)) throw IDASolverFailure (get_name(), status);")
            (#\})
            ))
      ((gsl)
        (pp indent 
            ( "/* Resets the GSL stepping function if any synaptic events have occurred */")
            ("if (events > 0)")
            (#\{)
            ("  gsl_odeiv2_step_reset (B_.s_);")
            (#\})
            ))

       )
    ))



(define (spike-handle method defaults vi)
  (case method
    ((cvode ida)  
#<<EOF
	      {
		set_spiketime(Time::ms(tt));
		SpikeEvent se;
                for (int i = 0; i < N; i++)
                {
                  S_.y_[i] = Ith(B_.y,i);
                }
		network()->send(*this, se, lag);
                adjust_zero_crossings(B_.y, 1e-7);
                continue;
	      }
EOF
        )

    (else
     (if (lookup-def 'V_t defaults)
         (sprintf 
#<<EOF
         if ( S_.r_ > 0 )
	    S_.r_--;
         else
	    if (S_.y_[~A] >= P_.V_t && V_.U_old_ > S_.y_[~A])
	      {
                S_.r_ = V_.RefractoryCounts_;
		set_spiketime(Time::step(origin.get_steps()+lag+1));
		SpikeEvent se;
		network()->send(*this, se, lag);
	      }
EOF
        vi vi)  ))

    ))



(define (output-update
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs transient-event-defs
         i-eqs pool-ions perm-ions synapse-info 
         event-external-eq-defs defaults imports
	 indent indent+)

  (let ((vi (lookup-def 'v state-index-map)))
    
    (pp indent  (,(sprintf "void ~A::update(Time const & origin, const long_t from, const long_t to)" sysname)
                 #\{))

    (pp indent+ 
#<<EOF
    assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
    assert(from < to);

    double tout;
    long_t current_steps = origin.get_steps();

    for ( long_t lag = from ; lag < to ; ++lag )
      {
        double h = B_.step_;    
	double tt = 0.0 ; 
EOF
    )

    (case method

      ((leapfrog)
       (begin
         (pp indent+ (,(leapfrog-solve vi (spike-handle method defaults vi))))
         ))

      ((cvode)
       (begin
         (pp indent+ (,(cvode-solve (spike-handle method defaults vi) )) ,nl)
         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
         ))

      ((ida)
       (begin
         (pp indent+ (,(ida-solve (spike-handle method defaults vi) )) ,nl)
         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
         ))

      ((gsl)
       (begin
         (pp indent+ (,(gsl-solve vi (spike-handle method defaults vi))) ,nl)
         (output-solver-events method transient-event-defs event-external-eq-defs synapse-info imports (+ 6 indent+) ) 
         ))
      )


    (pp (+ 6 indent+)
        ("B_.I_stim_ = B_.currents_.get_value(lag);")
        ("B_.logger_.record_data(current_steps + lag);"))

  (pp indent+  #\})
  (pp indent  #\})
  ))


(define (output-nodes-leapfrog
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info defaults
	 indent indent+)

  (let ((N (length state-index-map))
        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
        )

    (pp indent ,nl  (,#<#EOF
#{sysname}::#{sysname}()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}
EOF
  ))

  (pp indent  ,nl  (,#<#EOF
#{sysname}::#{sysname} (const #{sysname}& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}
EOF
  ))


  (pp indent #<#EOF
#{sysname}::~#{sysname} ()
{
    if ( B_.u != NULL ) free (B_.u);
    if ( B_.v != NULL ) free (B_.v);
    if ( B_.dvdt != NULL ) free (B_.dvdt);
}
EOF
  )

  (pp indent ,nl (,#<#EOF
  void #{sysname}::init_node_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}
EOF
  ))

  (pp indent ,nl (,#<#EOF
void #{sysname}::init_state_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    S_ = State_(pr.P_);
}
EOF
  ))
))


(define (output-nodes-cvode
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info defaults
	 indent indent+)

  (let (
        (N (length state-index-map))
        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
        )

    (pp indent ,nl  (,#<#EOF
#{sysname}::#{sysname}()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}
EOF
  ))


  (pp indent  ,nl  (,#<#EOF
#{sysname}::#{sysname}(const #{sysname}& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}
EOF
  ))

  (pp indent (,#<#EOF
#{sysname}::~#{sysname}()
{

  if ( B_.sys_ )
  {
    /* Free y vector */
    N_VDestroy_Serial(B_.y);

    /* Free integrator memory */
//    CVodeFree(&B_.sys_);
  }
}
EOF
   ))


  (pp indent ,nl (,#<#EOF
void #{sysname}::init_node_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}
EOF
  ))

  (pp indent ,nl (,#<#EOF
void #{sysname}::init_state_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    S_ = State_(pr.P_);
}
EOF
  ))

  (pp indent ,nl (,#<#EOF
void #{sysname}::init_buffers_()
{
EOF
  ))

  (for-each (lambda (psc)
              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
            pscs)

  (pp indent+ (#<<EOF
    B_.currents_.clear();           
    Archiving_Node::clear_history();

    B_.logger_.reset();

    B_.step_ = Time::get_resolution().get_ms();
    B_.IntegrationStep_ = B_.step_;

    B_.I_stim_ = 0.0;
EOF
))

   (pp indent+ (,#<#EOF

    int status, N, rootdir;

    N = #{N};
    // only positive direction (rising edge) of spike events will be detected
    rootdir = 1;

    /* Creates serial vector of length N */
    B_.y = N_VNew_Serial(N);
    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw CVodeSolverFailure (get_name(), 0);

    for (int i = 0; i < N; i++)
    {
       Ith(B_.y,i) = S_.y_[i];
    }
 
    /* Calls CVodeCreate to create the solver memory and specify the 
     * Backward Differentiation Formula and the use of a Newton iteration */
    B_.sys_ = CVodeCreate(CV_BDF, CV_NEWTON);
    if (check_flag((void *)B_.sys_, "CVodeCreate", 0)) throw CVodeSolverFailure (get_name(), 0);

   /* Calls CVodeInit to initialize the integrator memory and specify the
    * right hand side function in y'=f(t,y), the initial time, and
    * the initial values. */
    status = CVodeInit (B_.sys_, #{sysname}_dynamics, 0.0, B_.y);
    if (check_flag(&status, "CVodeInit", 1)) throw CVodeSolverFailure (get_name(), status);

EOF
    )

    ,(if (lookup-def 'V_t defaults) #<#EOF
    /* Spike event handler (detects zero-crossing of V-V_t) */
    status = CVodeRootInit(B_.sys_, 1, (CVRootFn)#{sysname}_event);
    if (check_flag(&status, "CVodeRootInit", 1)) throw CVodeSolverFailure (get_name(), status);

    /* Detect only the rising edge of spikes */
    status = CVodeSetRootDirection(B_.sys_, &rootdir);
    if (check_flag(&status, "CVodeSetRootDirection", 1)) throw CVodeSolverFailure (get_name(), status);
EOF
     '())

     (
#<<EOF
    /* Sets the relative and absolute error tolerances of CVode  */
    status = CVodeSStolerances (B_.sys_, 1e-7, 1e-7);
    if (check_flag(&status, "CVodeSStolerances", 1)) throw CVodeSolverFailure (get_name(), status);

    /* Turn on CVode stability limit detection (only applicable for order 3 and greater) */
    status = CVodeSetStabLimDet (B_.sys_,true);
    if (check_flag(&status, "CVodeSetStabLimDet", 1)) throw CVodeSolverFailure (get_name(), status);

    /* Sets the maximum order of CVode  */
    status = CVodeSetMaxOrd (B_.sys_,4);
    if (check_flag(&status, "CVodeSetMaxOrd", 1)) throw CVodeSolverFailure (get_name(), status);

    /* Sets maximum step size. */
    status = CVodeSetMaxStep (B_.sys_,B_.step_);
    if (check_flag(&status, "CVodeSetMaxStep", 1)) throw CVodeSolverFailure (get_name(), status);

    /* Calls CVodeSetUserData to configure the integrator to pass the 
     * params structure to the right-hand function */
    status = CVodeSetUserData(B_.sys_, reinterpret_cast<void*>(this));
    if (check_flag(&status, "CVodeSetUserData", 1)) throw CVodeSolverFailure (get_name(), status);

    /* Initializes diagonal linear solver. */
    status = CVDiag (B_.sys_);
    if (check_flag(&status, "CVDiag", 1)) throw CVodeSolverFailure (get_name(), status);
    }

EOF
    ))

  (pp indent (,#<#EOF
void #{sysname}::calibrate()
{
    B_.logger_.init();  
}

EOF
   ))
  
))


(define (output-nodes-ida
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info defaults
	 indent indent+)

  (let (
        (N (length state-index-map))
        (pscs (lookup-def 'post-synaptic-conductances synapse-info))
        )

    (pp indent ,nl  (,#<#EOF
#{sysname}::#{sysname}()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}
EOF
  ))


  (pp indent  ,nl  (,#<#EOF
#{sysname}::#{sysname}(const #{sysname}& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}
EOF
  ))

  (pp indent (,#<#EOF
#{sysname}::~#{sysname}()
{

  if ( B_.sys_ )
  {
    /* Free y vector */
    N_VDestroy_Serial(B_.y);
    N_VDestroy_Serial(B_.yp);

    /* Free integrator memory */
//    IDAFree(&B_.sys_);
  }
}
EOF
   ))


  (pp indent ,nl (,#<#EOF
void #{sysname}::init_node_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}
EOF
  ))

  (pp indent ,nl (,#<#EOF
void #{sysname}::init_state_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    S_ = State_(pr.P_);
}
EOF
  ))

  (pp indent ,nl (,#<#EOF
void #{sysname}::init_buffers_()
{
EOF
  ))

  (for-each (lambda (psc)
              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
            pscs)

  (pp indent+ (#<<EOF
    B_.currents_.clear();           
    Archiving_Node::clear_history();

    B_.logger_.reset();

    B_.step_ = Time::get_resolution().get_ms();
    B_.IntegrationStep_ = B_.step_;

    B_.I_stim_ = 0.0;
EOF
))

   (pp indent+ (,#<#EOF

    int status, N, rootdir;

    N = #{N};
    // only positive direction (rising edge) of spike events will be detected
    rootdir = 1;

    /* Creates serial vector of length N */
    B_.y = N_VNew_Serial(N);
    B_.yp = N_VNew_Serial(N);
    if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw IDASolverFailure (get_name(), 0);

    for (int i = 0; i < N; i++)
    {
       Ith(B_.y,i) = S_.y_[i];
    }
 
    /* Calls IDACreate to create the solver memory */ 
    B_.sys_ = IDACreate();
    if (check_flag((void *)B_.sys_, "IDACreate", 0)) throw IDASolverFailure (get_name(), 0);

   /* Calls IDAInit to initialize the integrator memory and specify the
    * resdual function, the initial time, and the initial values. */
    status = CVodeInit (B_.sys_, #{sysname}_dynamics, 0.0, B_.y);
    status = IDAInit (B_.sys, #{sysname}_residual, 0.0, B_.y, B_.yp);

    if (check_flag(&status, "IDAInit", 1)) throw IDASolverFailure (get_name(), status);

EOF
    )

    ,(if (lookup-def 'V_t defaults) #<#EOF
    /* Spike event handler (detects zero-crossing of V-V_t) */
    status = IDARootInit(B_.sys_, 1, (IDARootFn)#{sysname}_event);
    if (check_flag(&status, "CVodeRootInit", 1)) throw IDASolverFailure (get_name(), status);

    /* Detect only the rising edge of spikes */
    status = IDASetRootDirection(B_.sys_, &rootdir);
    if (check_flag(&status, "CVodeSetRootDirection", 1)) throw IDASolverFailure (get_name(), status);
EOF
     '())

     (
#<<EOF
    /* Sets the relative and absolute error tolerances of IDA  */
    status = IDASStolerances (B_.sys_, 1e-7, 1e-7);
    if (check_flag(&status, "IDASStolerances", 1)) throw IDASolverFailure (get_name(), status);

    /* Sets the maximum order of IDA  */
    status = IDASetMaxOrd (B_.sys_,4);
    if (check_flag(&status, "IDASetMaxOrd", 1)) throw IDASolverFailure (get_name(), status);

    /* Sets maximum step size. */
    status = IDASetMaxStep (B_.sys_,B_.step_);
    if (check_flag(&status, "IDASetMaxStep", 1)) throw IDASolverFailure (get_name(), status);

    /* Calls IDASetUserData to configure the integrator to pass the 
     * params structure to the right-hand function */
    status = IDASetUserData(B_.sys_, reinterpret_cast<void*>(this));
    if (check_flag(&status, "IDASetUserData", 1)) throw IDASolverFailure (get_name(), status);

    /* Initializes dense linear solver. */
    status = IDADense (B_.sys_);
    if (check_flag(&status, "IDADense", 1)) throw IDASolverFailure (get_name(), status);
    }

EOF
    ))

  (pp indent (,#<#EOF
void #{sysname}::calibrate()
{
    B_.logger_.init();  
}

EOF
   ))
  
))


(define (output-nodes-gsl
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info defaults
	 indent indent+)

  (let ((N    (length state-index-map))
        (iv   (lookup-def 'v state-index-map))
        (pscs (lookup-def 'post-synaptic-conductances synapse-info)))

    (pp indent ,nl  (,#<#EOF
#{sysname}::#{sysname}()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}
EOF
  ))

  (pp indent  ,nl  (,#<#EOF
#{sysname}::#{sysname}(const #{sysname}& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}
EOF
  ))


  (pp indent ,#<#EOF
#{sysname}::~#{sysname} ()
{
    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
    if ( B_.s_ != NULL) gsl_odeiv2_step_free (B_.s_);
    if ( B_.c_ != NULL) gsl_odeiv2_control_free (B_.c_);
    if ( B_.e_ != NULL) gsl_odeiv2_evolve_free (B_.e_);
    if ( B_.u != NULL) free (B_.u);
    if ( B_.jac != NULL) free (B_.jac);
}
EOF
  )

  (pp indent ,nl (,#<#EOF
  void #{sysname}::init_node_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}
EOF
  ))

  (pp indent ,nl (,#<#EOF
void #{sysname}::init_state_(const Node& proto)
{
    const #{sysname}& pr = downcast<#{sysname}>(proto);
    S_ = State_(pr.P_);
}
EOF
  ))

  (pp indent ,nl (,#<#EOF
void #{sysname}::init_buffers_()
{
EOF
  ))

  (for-each (lambda (psc)
              (pp indent+ (,(sprintf "B_.spike_~A.clear();" (second psc)))))
            pscs)

  (pp indent+ (#<<EOF
    B_.currents_.clear();           
    Archiving_Node::clear_history();

    B_.logger_.reset();

    B_.step_ = Time::get_resolution().get_ms();
    B_.IntegrationStep_ = B_.step_;

    B_.I_stim_ = 0.0;
EOF
))

    (pp indent+ ( 
,#<#EOF
    static const gsl_odeiv2_step_type* T1 = gsl_odeiv2_step_rk2;
    B_.N = #{N};
  
    if ( B_.s_ == 0 )
      B_.s_ = gsl_odeiv2_step_alloc (T1, B_.N);
    else 
      gsl_odeiv2_step_reset(B_.s_);
    
    if ( B_.c_ == 0 )  
      B_.c_ = gsl_odeiv2_control_standard_new (1e-7, 1e-7, 1.0, 0.0);
    else
      gsl_odeiv2_control_init(B_.c_, 1e-7, 1e-7, 1.0, 0.0);
    
    if ( B_.e_ == 0 )  
      B_.e_ = gsl_odeiv2_evolve_alloc(B_.N);
    else 
      gsl_odeiv2_evolve_reset(B_.e_);
  
    B_.sys_.function  = #{sysname}_dynamics; 
    B_.sys_.jacobian  = #{sysname}_jacobian;
    B_.sys_.dimension = B_.N;
    B_.sys_.params    = reinterpret_cast<void*>(this);

    B_.u = (double *)malloc(sizeof(double) * B_.N);
    assert (B_.u);
    B_.jac = (double *)malloc(sizeof(double) * B_.N);
    assert (B_.jac);

EOF
   ))

  (pp indent ,nl #\})


  (pp indent (,#<#EOF
void #{sysname}::calibrate()
{
    B_.logger_.init();  
    #{(if iv (sprintf "V_.U_old_ = S_.y_[~A];" iv) "")}
    #{(if (lookup-def 't_ref defaults) 
          "V_.RefractoryCounts_ = Time(Time::ms(P_.t_ref)).get_steps();" 
          "V_.RefractoryCounts_ = 0;")}
}
EOF
))
))



(define (output-buffers
	 sysname method state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions 
	 synapse-info defaults
	 indent indent+)

  (let ((N (length state-index-map)))

    (case method

      ((leapfrog)

       (pp indent ,nl
	   (,#<#EOF
#{sysname}::Buffers_::Buffers_(#{sysname}& n)
    : logger_(n),
      sys_(0),
      N(0),
      u(0),
      v(0),
      dvdt(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
)))

   ((cvode ida)

    (pp indent ,nl
        (,#<#EOF
#{sysname}::Buffers_::Buffers_(#{sysname}& n)
    : logger_(n),
      sys_(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
        )
       (,#<#EOF
#{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n)
    : logger_(n),
      sys_(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
    )))

    ((gsl)

     (pp indent ,nl
	 ( ,#<#EOF
#{sysname}::Buffers_::Buffers_(#{sysname}& n)
    : logger_(n),
      s_(0),
      c_(0),
      e_(0),
      N(0),
      u(0),
      jac(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}

EOF
          )
        ( ,#<#EOF
#{sysname}::Buffers_::Buffers_(const Buffers_&, #{sysname}& n)
    : logger_(n),
      s_(0),
      c_(0),
      e_(0),
      N(0),
      u(0),
      jac(0)
{
    // Initialization of the remaining members is deferred to
    // init_buffers_().
}
EOF
    )))
   ))
)


(define (output-accessors+modifiers
	 sysname imports defaults state-index-map const-defs asgn-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions constraints
         indent indent+)

  (let ((c-eqs (lookup-def 'c-eqs constraints))
        
        (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
                                     (list (nest-name n) v)))
                      (lookup-def 'c-units constraints)))
        )

    (pp indent ,nl (,(sprintf "void ~A::Parameters_::get (DictionaryDatum &d) const" sysname) ))
    (pp indent  #\{)
    
    (for-each
     (lambda (def)
       (let ((n (first def)))
         (pp indent+ (,(sprintf "def<double_t>(d, ~S, ~A);" (->string n) n)))))
     (append const-defs
             defaults))
    
    (pp indent  #\})
    
    (pp indent ,nl (,(sprintf "void ~A::Parameters_::set (const DictionaryDatum &d)" sysname) ))
    (pp indent  #\{)
    
    (for-each
     (lambda (def)
       (let* ((n (first def))
              (nu (lookup-def n c-units))
              (scale (and nu (nemo:unit-scale nu)))
              )
         (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, ~A);" (->string n) n)))
         (if scale (pp indent+ (,(sprintf "~A = ~A * ~A;" n scale n ))))
         ))
     (append const-defs
             defaults))
    
    (pp indent  #\})
    
    (pp indent ,nl (,(sprintf "void ~A::State_::get (DictionaryDatum &d) const" sysname) ))
    (pp indent  #\{)
    
    (for-each
     (lambda (def)
       (let* (
              (n (first def))
              (i (second def))
              )
         (pp indent+ (,(sprintf "def<double_t>(d, ~S, y_[~A]);" (->string n) i)))
         ))
     state-index-map)
    
    (let ((vi (lookup-def 'v state-index-map)))
      (if vi
          (pp indent+ (,(sprintf "def<double_t>(d, names::V_m, y_[~A]);"  vi) ))
          ))
    
    (pp indent  #\})
    
    (pp indent ,nl (,(sprintf "void ~A::State_::set (const DictionaryDatum &d, const Parameters_&)" sysname) ))
    (pp indent  #\{)
    
    (for-each
     (lambda (def)
       (let* ((n     (first def)) 
              (i     (second def))
              (nu    (lookup-def n c-units))
              (scale (and nu (nemo:unit-scale nu)))
              )
         (pp indent+ (,(sprintf "updateValue<double_t>(d, ~S, y_[~A]);" (->string n) i)))
         (if scale (pp indent+ (,(sprintf "y_[~A] = ~A * y_[~A];" i scale i ))))
         ))
     state-index-map)
    
    (let ((vi (lookup-def 'v state-index-map)))
      (if vi
          (pp indent+ (,(sprintf "updateValue<double_t>(d, names::V_m, y_[~A]);"  vi) ))
          ))
    
    (pp indent  #\})
    
    ))


(define (output-init sysname method state-index-map steady-state-index-map 
		     imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
		     reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints 
                     indent indent+)

  
  (let* ((N (length state-index-map))

         (c-eqs (lookup-def 'c-eqs constraints))

         (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
                                     (list (nest-name n) v)))
                       (lookup-def 'c-units constraints)))

	 (i-eqs 
	  (map 
	   (lambda (def) (list (first def) 
			       (add-params-to-fncall 
				(canonicalize-expr/C++ (second def)) builtin-fns)))
	   i-eqs))

	 (init-eqs 
	  (append 
	     
	   (map (lambda (def)
		  (let ((n (first def))
			(b (second def)))
		    (list (nest-name n) (nest-name b))))
		external-eq-defs)
	   
	   asgn-eq-defs
	   init-eq-defs

	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-in pool-ion))
			(b (pool-ion-inq pool-ion)))
		    (list n b)))
		pool-ions)

	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-out pool-ion))
			(b (pool-ion-outq pool-ion)))
		    (list n b)))
		pool-ions)
	   ))

	 (init-dag 
	  (map (lambda (def)
		 (cons (first def) (enum-freevars (second def) '() '())))
	       init-eqs))

	 (init-order
	  (reverse
	   (topological-sort init-dag 
			     (lambda (x y) (string=? (->string x) (->string y))))))

	 (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))

	 )
    
    (pp indent ,nl
        (,(sprintf "~A::State_::State_" sysname) (const Parameters_& p))
        ( #\{))

    (pp indent+
	(double ,(slp ", " (delete-duplicates
			    (map ->string 
				 (filter (lambda (x) (not (member x builtin-consts)))
					 (append 
					  init-locals
					  init-order
					  (map first external-eq-defs)
					  (map pool-ion-in pool-ions)
					  (map pool-ion-out pool-ions)
					  (map first i-eqs)
					  (map first state-index-map) 
					  (map first const-defs)
					  (map first reaction-eq-defs)
					  )))
				 string=?))
		      ";")
	("const Parameters_ *params = &p;")
	)

    (case method
      ((gsl leapfrog) 
       (pp indent+ "r_ = 0;" ,nl)))

    (pp indent+ ,(sprintf "memset(y_,0,~A*sizeof(double));" N))

    (for-each
     (lambda (x) 
       (let* ((n  (first x)))
         (pp indent+ ,(expr->string/C++ (sprintf "p.~A" n) n))))
     const-defs)

    (for-each
     (lambda (eq) 
       (match-let (((op left right)  eq))
                  (pp indent+ ,(sprintf "if (!((~A) ~A (~A))) { throw BadProperty (~S); }; " 
                                        (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ left))) op
                                        (expr->string/C++ (canonicalize-expr/C++ (rhsexpr/C++ right)))
                                        (sprintf "Constraint ~A is not satisfied." eq)
                                        ))
                  ))
     c-eqs)

    (let ((vi (lookup-def 'v state-index-map))
	  (vrest (or (and (lookup-def 'Vrest const-defs) 'Vrest) -65.0)))
      (if (and vi vrest) 
	  (pp indent+ ,(expr->string/C++  vrest 'v ))))

    (for-each (lambda (n)
		(let ((b  (lookup-def n init-eqs)))
		  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
	      init-order)
    
    (for-each (lambda (def)
		(let* ((n  (first def)) 
		       (ni (lookup-def n state-index-map)))
		  (if ni (pp indent+ ,(expr->string/C++ n  (sprintf "y_[~A]" ni))))))
	      init-eq-defs)

    
#|
    (if (not (null? steady-state-index-map))
	(begin
	  (gsl-fminimizer sysname N indent+)
	  (for-each
	   (lambda (def)
	     (let* ((n (first def)) 
		    (si (lookup-def n steady-state-index-map))
		    (ni (lookup-def n state-index-map)))
	       (if si (begin 
			(pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n))
			(pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")") 
						       (s+ "y_[" ni "]")))))))
	   rate-eq-defs)
      (pp indent+ "gsl_vector_free (ys);")
      ))
|#

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

    (let ((vi (lookup-def 'v state-index-map)))
      (if vi (pp indent+ ,(expr->string/C++ 'v (sprintf "y_[~A]" vi)))))
      
    (pp indent  #\})

    (pp indent ,nl
	(,(sprintf "~A::State_::State_" sysname) (const State_& s)  )
        ( #\{)
        
        ,(case method
          ((gsl leapfrog) 
           `("r_ = s.r_;"))
          (else ""))
        
        (,(sprintf "for ( int i = 0 ; i < ~A ; ++i ) y_[i] = s.y_[i];" N))
        ( #\})
        ,nl
        )

    (pp indent ,nl
        (,(sprintf "~A::State_& ~A::State_::operator=(const State_& s)" sysname sysname))
        ( #\{)
        )

    (case method
      ((gsl leapfrog) 
       (pp indent+ "r_ = s.r_;" ,nl)))

    (pp indent+ ( ,#<#EOF
     assert(this != &s);  
     for ( size_t i = 0 ; i < #{N} ; ++i )
       y_[i] = s.y_[i];
EOF
))

        (pp indent+ "return *this;")

        (pp indent  #\})
))


(define (output-parameters sysname imports const-defs external-eq-defs defaults indent indent+)

  (let* ((parameter-defs
	  (map 
	   (lambda (def) (list (first def) 
			       (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
	   const-defs))
	 (parameter-locals  (find-locals (map second parameter-defs)))

         (const-exprs
	    (map (lambda (def)
		   (let* ((n  (first def)) (b (second def)))
		     (s+ (nest-name n) "  (" (expr->string/C++ b) ")")))
		 const-defs) )
	     
	  (default-exprs 
	    (map (lambda (def)
		   (let ((n (first def))
			 (b (second def)))
		     (expr->string/C++ (nest-name b) (nest-name n))))
		 defaults))
	   )

      (pp indent ,nl (,(sprintf "~A::Parameters_::Parameters_" sysname) () ":"))
      
      (if (not (null? parameter-locals)) 
	  (pp indent+ (double ,(slp ", " parameter-locals) ";")))

      (pp indent+ ,(slp ",\n" const-exprs))
      
      (pp indent  #\{ ,(slp "\n" default-exprs) #\})
      
      ))


(define (output-recordables
	 sysname state-index-map steady-state-index-map 
	 const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	 reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
  
  (pp indent ,nl (,(sprintf "RecordablesMap<~A> ~A::recordablesMap_;" sysname sysname)))
  (pp indent (,(sprintf "template <> void RecordablesMap<~A>::create()" sysname)))
  (pp indent #\{)

  (for-each 
   (lambda (def)
     (let ((n (first def)) (i (second def)))
       (pp indent+ (,(sprintf "insert_(~S, &~A::get_y_elem_<~A::State_::~A>);" 
			      (->string n) sysname sysname (string-upcase (->string n)))))
       ))
   state-index-map)

  (if (lookup-def 'v state-index-map)
      (pp indent+ (,(sprintf "insert_(names::V_m, &~A::get_y_elem_<~A::State_::~A>);" 
			     sysname sysname (string-upcase (->string 'v))))))

  (pp indent #\})
  )



(define (output-dy sysname method imports const-defs state-index-map 
		   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
		   pool-ions i-eqs v-eq constraints indent indent+)

  (let* (
         (c-units (map (lambda (x) 
                         (let ((n (first x)) (v (second x)))
                           (list (nest-name n) v)))
                       (lookup-def 'c-units constraints)))

	 (i-eqs 
	  (map 
	   (lambda (def) (list (first def) 
			       (add-params-to-fncall (canonicalize-expr/C++ (second def)) builtin-fns)))
	   i-eqs))

	 (v-eq
	  (and v-eq
	       (list (first v-eq) 
		     (add-params-to-fncall (canonicalize-expr/C++ (second v-eq)) builtin-fns))))

	 (eqs 
	  (append 

	   external-eq-defs
	   asgn-eq-defs
	   reaction-eq-defs
		 
	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-in pool-ion))
			(b (pool-ion-inq pool-ion)))
		    (list n b)))
		pool-ions)

	   (map (lambda (pool-ion)
		  (let ((n (pool-ion-out pool-ion))
			(b (pool-ion-outq pool-ion)))
		    (list n b)))
		pool-ions)

           i-eqs
	  ))
	 
	 (eq-dag 
	  (map (lambda (def)
		 (cons (first def) (enum-freevars (second def) '() '())))
	       eqs))

	 (eq-order
	  (reverse
	   (topological-sort eq-dag 
			     (lambda (x y) (string=? (->string x) (->string y))))))

	 (eq-locals  (find-locals 
		      (map second
			   (if v-eq (cons v-eq (append i-eqs rate-eq-defs eqs))
			       (append i-eqs rate-eq-defs eqs)))))
	 )


    (case method
      ((leapfrog)
       (pp indent ,nl (
                       ,(sprintf "extern \"C\" int ~A_dynamics" sysname)
                       (,(slp ", " `("double t"
                                     "const double y[]"
                                     "double f[]"
                                     "void* pnode"
                                     )))
                       #\{
                       )))
      ((cvode)
       (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname)
                       (,(slp ", " `("double t"
                                     "N_Vector y"
                                     "N_Vector f"
                                     "void* pnode"
                                     )))
                       #\{
                       )))
      ((ida)
       (pp indent ,nl (,(sprintf "extern \"C\" int ~A_residual" sysname)
                       (,(slp ", " `("double t"
                                     "N_Vector y"
                                     "N_Vector yp"
                                     "N_Vector f"
                                     "void* pnode"
                                     )))
                       #\{
                       )))
      ((gsl)
       (pp indent ,nl (,(sprintf "extern \"C\" int ~A_dynamics" sysname)
                       (,(slp ", " `("double t"
                                     "const double y[]"
                                     "double f[]"
                                     "void* pnode"
                                     )))
                       #\{
                       )))
      )



    (pp indent+ (double ,(slp ", " (delete-duplicates 
				    (map (compose ->string nest-name)
					 (filter (lambda (x) 
						   (not (member x builtin-consts)))
						 (append 
						  eq-locals
						  eq-order
						  (map first i-eqs)
						  (map first external-eq-defs)
						  (map first state-index-map)
						  (map first const-defs)
						  )))
				    string=?))
			";"))

    (pp indent+ ,nl

        ("// S is shorthand for the type that describes the model state ")
	(typedef ,(sprintf "~A::State_" sysname) "S;")
	
	,nl
	
	("// cast the node ptr to an object of the proper type")
	("assert(pnode);")
	("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
	(,sysname "& vnode =  *(reinterpret_cast<" ,sysname "*>(pnode));")
	
	,nl

	("// params is a reference to the model parameters ")
	(const struct ,(sprintf "~A::Parameters_" sysname) "*params;")
	("params = &(node.P_);")
	
	,nl

	("// state is a reference to the model state ")
	(struct ,(sprintf  "~A::State_" sysname) "*state;")
	("state = &(vnode.S_);")
	
	,nl

	("// y[] must be the state vector supplied by the integrator, ")
	("// not the state vector in the node, node.S_.y[]. ")

	,nl

	)

    (for-each (lambda (def)
		(let ((n (first def)) )
		  (pp indent+
		      ,(expr->string/C++ (sprintf "params->~A" n) n))))
	      const-defs)

    (let ((vi (lookup-def 'v state-index-map)))
      (if vi (pp indent+ ,(expr->string/C++ (ith 'y vi) 'v))))

    (for-each (lambda (def)
		(let* ((n (first def))
		       (ni (lookup-def n state-index-map)))
		  (pp indent+ ,(expr->string/C++ (ith 'y ni) (nest-name n)))))
	      rate-eq-defs)

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

    (for-each (lambda (def)
		(let* ((n (first def))
		       (b (second def))
		       (fv (ith 'f (lookup-def n state-index-map)))
		       )
		  (pp indent+ ,(s+ "// rate equation " n)
		      ,(expr->string/C++ b fv))))
	      rate-eq-defs)
    
    (if v-eq
	(let ((vi (lookup-def 'v state-index-map)))
	  (if vi
	      (pp indent+ ,(expr->string/C++ (second v-eq) (ith 'f vi)))
	      )))

    (case method
      ((leapfrog)
	(pp indent+ ,nl ("return 0;")))
      ((cvode)
	(pp indent+ ,nl ("return 0;")))
      ((gsl)
	(pp indent+ ,nl ("return GSL_SUCCESS;")))
      )

    (pp indent  #\})

    ))


(define (nemo:nest-translator sys . rest)

  (define (cid x)  (second x))
  (define (cn x)   (first x))

  (let-optionals rest ((dirname ".") (method #f))

    (let ((method (or method 'gsl)))

      (if (not (member method '(gsl cvode leapfrog)))
	  (nemo:error 'nemo:nest-translator ": unknown method " method))

  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
    (let ((imports  ((dis 'imports)  sys))
	  (exports  ((dis 'exports)  sys)))
      (let* ((indent      0)
	     (indent+     (+ 2 indent ))

	     (sysname     (nest-name ((dis 'sysname) sys)))
	     (prefix      (->string sysname))
	     (deps*       ((dis 'depgraph*)  sys))
	     (consts      ((dis 'consts)     sys))
	     (asgns       ((dis 'asgns)      sys))
	     (states      ((dis 'states)     sys))
	     (reactions   ((dis 'reactions)  sys))
	     (defuns      ((dis 'defuns)     sys))
	     (components  ((dis 'components) sys))
	     
	     (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
	     (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))

	     (const-defs       (filter-map
				(lambda (nv)
				  (and (not (member (first nv) builtin-consts))
				       (let ((v1 (canonicalize-expr/C++ (second nv))))
					 (list (nest-name (first nv)) v1))))
				consts))
	     
	     (defaults             (nemo:defaults-query sys))

             (geometry             (nemo:geometry-query sys))

	     (gate-complex-info    (nemo:gate-complex-query sys))
	     (perm-ions       (map (match-lambda ((comp i e erev val) `(,comp ,(nest-name i) ,(nest-name e) ,erev)))
				   (lookup-def 'perm-ions gate-complex-info)))
	     (acc-ions        (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out)))))
				   (lookup-def 'acc-ions gate-complex-info)))
	     (epools          (lookup-def 'pool-ions gate-complex-info))
	     (pool-ions       (pool-ion-name-map nest-name  epools))

	     (comprc         (any (match-lambda ((name 'membrane-tau id) (list name id)) (else #f)) components))
	     (compcap        (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
	     (mrc            (or (and comprc (car ((dis 'component-exports) sys (cid comprc))))
                                 (lookup-def 'membrane-tau defaults)
                                 (lookup-def 'tau_m defaults)
				 (and compcap (car ((dis 'component-exports) sys (cid compcap))))
                                 (lookup-def 'membrane-capacitance defaults)
                                 (lookup-def 'C_m defaults)
                                 ))

	     (soma-geometry  (lookup-def 'soma geometry))
	     (marea          (and soma-geometry (third soma-geometry)))

	     (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
	     (synapse-info         (nemo:post-synaptic-conductance-query sys))


             (pscs           (lookup-def 'post-synaptic-conductances synapse-info))

	     (i-syns         (lookup-def 'i-synapses synapse-info))
		
	     (i-gates        (lookup-def 'i-gates gate-complex-info))

	     (i-defs         (nemo:ionic-current-definitions
                              gate-complexes i-gates i-syns pscs marea
                              (lambda (x) (state-power sys x))
                              (lambda (x) ((dis 'component-exports) sys x))
                              (lambda (x) ((dis 'component-subcomps) sys x))
                              nest-name rhsexpr/C++ canonicalize-expr/C++
                              builtin-fns))

             (i-eqs          (lookup-def 'i-eqs i-defs))
             (i-names        (lookup-def 'i-names i-defs))

             (constraints    (nemo:constraint-definitions 
                              gate-complexes i-gates i-syns pscs marea imports
                              (lambda (x) (state-power sys x))
                              (lambda (x) (quantity-unit sys x))
                              (lambda (x) ((dis 'component-exports) sys x))
                              (lambda (x) ((dis 'component-subcomps) sys x))
                              nest-name))

	     (external-eq-defs   (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++
                                                        namespace-filter: (lambda (x) (not (equal? x 'event)))))

	     (event-external-eq-defs   (sys->external-eq-defs sys nest-name rhsexpr/C++ canonicalize-expr/C++
                                                              namespace-filter: (lambda (x) (equal? x 'event))))

	     (asgn-eq-defs       (poset->asgn-eq-defs* poset sys nest-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
	     
	     (rate-eq-defs       (reverse (poset->rate-eq-defs* poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)))
	     
	     (reaction-eq-defs   (poset->reaction-eq-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++))

	     (transient-event-defs  (poset->transient-event-defs poset sys method nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns)) 
	     
	     (init-eq-defs       (poset->init-defs* poset sys nest-name nest-state-name rhsexpr/C++ canonicalize-expr/C++ builtin-fns))
	     
	     (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
				      (poset->state-conserve-eq-defs poset sys nest-name nest-state-name)))
	     
	     (imports-sans-v (filter (lambda (x) (not (equal? 'v (first x)))) imports))

	     (v-eq    (and (not (null? i-names))
			   (let ((istim "(node.B_.I_stim_)" )) 
			     (cond

                              ((and mrc marea)
                               (list 'v (rhsexpr/C++ 
                                         `(/ (+ (* ,istim (/ 100. ,marea)) 
                                                (* -1e3 ,(sum i-names))) ,mrc))))
                              (marea
                               (list 'v (rhsexpr/C++ 
                                         `(+ (* ,istim (/ 100. ,marea))
                                             (* -1e-3 ,(sum i-names))))))
                              (mrc
                               (list 'v (rhsexpr/C++ `(/ (+ ,istim (* -1e-3 ,(sum i-names))) ,mrc))))
                              
                              (else
                               (list 'v (rhsexpr/C++ `(+ ,istim (* -1e-3 ,(sum i-names))))))
                              ))
                           ))

	     
	     (state-index-map  (let ((acc (fold (lambda (def ax)
						  (let ((st-name (first def)))
						    (list (+ 1 (first ax)) 
							  (cons `(,st-name ,(first ax)) (second ax)))))
						(list 0 (list)) 
						(if (lookup-def 'v imports)
						    (cons (list 'v) rate-eq-defs)
						    rate-eq-defs)
						)))
				 (second acc)))
	     
	     (steady-state-index-map  (let ((acc (fold (lambda (def ax)
							 (let ((st-name (first def)))
							   (if (not (alist-ref st-name init-eq-defs))
							       (list (+ 1 (first ax)) 
								     (cons `(,st-name ,(first ax)) (second ax)))
							       ax)))
						       (list 0 (list)) 
						       rate-eq-defs)))
					(second acc)))
	     
	     (dfenv 
	      (map (lambda (x) (let ((n (first x))) (list n (nest-name (sprintf "d_~A" n )))))
		   defuns))

	     )
	
	
	
	(for-each 
	 (lambda (a)
	   (let ((acc-ion   (car a)))
	     (if (assoc acc-ion perm-ions)
		 (nemo:error 'nemo:nest-translator 
			     ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
	 acc-ions)


	(let ((cpp-output  (open-output-file (make-output-fname dirname prefix ".cpp")))
	      (hpp-output  (open-output-file (make-output-fname dirname prefix ".h"))))
	  
	  ;; include files and other prelude definitions
	  (with-output-to-port cpp-output
	    (lambda ()
              (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
	      (output-prelude sysname steady-state-index-map indent)
	      ))
	  
	  (with-output-to-port cpp-output
	    (lambda () (pp indent ,nl "namespace nest {" ,nl)))
	  
          (case method
            ((cvode) 
             (with-output-to-port cpp-output
               (lambda ()
                 (pp indent ,sundials-prelude)))))

	  ;; user-defined functions
	  (let ((define-fn  (make-define-fn sysname))
		(define-fn-header (make-define-fn-header sysname)))
	    
	    (for-each (lambda (fndef) 
			(if (not (member (car fndef) builtin-fns))
			    (with-output-to-port cpp-output
			      (lambda ()
				(apply define-fn-header (cons indent fndef))
				(pp indent ,nl)))
			    ))
		      defuns)

	    (for-each (lambda (fndef) 
			(if (not (member (car fndef) builtin-fns))
			    (with-output-to-port cpp-output
			      (lambda ()
				(apply define-fn (cons indent fndef))
				(pp indent ,nl)))
			    ))
		      defuns)
	    )
	  
	  
	  ;; derivative function
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-dy sysname method imports-sans-v const-defs state-index-map 
			 external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
			 pool-ions i-eqs v-eq constraints
			 indent indent+)
	      ))

	  ;; Jacobian function
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-jac  sysname method state-index-map pool-ions i-eqs v-eq indent indent+)
	      ))

	  ;;  system recordables
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-recordables
	       sysname state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; system parameters
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-parameters sysname imports-sans-v const-defs 
				 external-eq-defs defaults indent indent+)
	      ))
	  
	  ;; initial values function
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-init sysname method state-index-map steady-state-index-map 
                           imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
                           reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; accessors and modifiers for states and parameters
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-accessors+modifiers
	       sysname imports-sans-v defaults state-index-map 
	       const-defs asgn-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions 
               constraints
	       indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  
	  ;; buffer and node initialization
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-buffers
	       sysname method state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions 
	       synapse-info defaults
	       indent indent+)
	      (pp indent ,nl)
	      ))
	  

	  (with-output-to-port cpp-output
	    (lambda ()
              ((case method
                 ((leapfrog) 
                  output-nodes-leapfrog)
                 ((gsl)      
                  output-nodes-gsl)
                 ((cvode)
                  output-nodes-cvode)
                 ((ida)
                  output-nodes-ida)
                 )
	       sysname method state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions 
	       synapse-info defaults
	       indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; spike threshold detect (cvode method only)
          (case method
            ((cvode)
              (with-output-to-port cpp-output
                (lambda ()
                  (output-cvode-event
                   sysname method imports const-defs state-index-map 
                   external-eq-defs rate-eq-defs reaction-eq-defs asgn-eq-defs
                   defaults i-eqs v-eq indent indent+)
                  (pp indent ,nl)
                  )))
            )

          ;; synaptic event transients

	  (with-output-to-port cpp-output
	    (lambda ()
              (output-synapse-transients
               sysname method state-index-map
               const-defs transient-event-defs event-external-eq-defs 
               synapse-info imports constraints
               indent indent+)
	      (pp indent ,nl)
              ))

	  ;; state update
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-update
	       sysname method state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs transient-event-defs i-eqs
               pool-ions perm-ions synapse-info 
               event-external-eq-defs
               defaults imports indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  ;; spike handler
	  (with-output-to-port cpp-output
	    (lambda ()
	      (output-event-handle
	       sysname state-index-map steady-state-index-map 
	       const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
	       reaction-eq-defs i-eqs pool-ions perm-ions 
	       synapse-info
	       indent indent+)
	      (pp indent ,nl)
	      ))
	  
	  (with-output-to-port cpp-output
	    (lambda () (pp indent "}" ,nl)))

	     (with-output-to-port hpp-output
	       (lambda ()
                 (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
		 (output-header
		  sysname method state-index-map steady-state-index-map defaults
		  external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
		  reaction-eq-defs i-eqs pool-ions perm-ions 
		  synapse-info
		  indent indent+)
		 (pp indent ,nl)
		 ))
	     
	     (close-output-port cpp-output)
	     (close-output-port hpp-output)
             
             ))
      ))
  ))

))