;;
;; Leaky integrate-and-fire model.
;; 

(use signal-diagram signal-diagram-dynamics mathh srfi-1 )

(define (LeakyIF:construct t V gL vL Isyn C theta Vreset t_rpend tau_rp )

  `((f . ,(let* ((subthreshold-eq         (UNION 
                                           (ODE (variable h) t 
                                                (V  (/ (+ (* (neg gL) (- V vL)) Isyn) C)))
                                           (ASSIGN (refractoryEnd refractoryEnd)
                                                   (t_rpend t_rpend))))

		 (threshold-detect        (ASSIGN (spike (- V theta))))

		 (refractory-eq           (SEQUENCE
                                           (ON (ASSIGN (V Vreset)
                                                       (t_rpend (+ t tau_rp))) 'spike)
                                           (UNION 
                                            (ODE (variable h) t
                                                 (V 0.0))
                                            (ASSIGN (spike spike))
                                            )))

		 (refractory-detect       (ASSIGN (refractoryEnd (- t t_rpend))))
		 )

	    (RTRANSITION subthreshold-eq
                         refractory-eq
                         'spike threshold-detect
                         'refractoryEnd refractory-detect 'st)

	    ))

    (initial .
	     ((h              1e-2)
              (t               ,t)
              (V               ,V)
              (theta           ,theta)
              (C               ,C)
              (Vreset          ,Vreset)
              (gL              ,gL)
              (vL              ,vL)
              (t_rpend         ,t_rpend)
              (tau_rp          ,tau_rp)
              (Isyn            ,Isyn)
              (spike           -1.0)
              (refractoryEnd   -1.0)
              (st              #f)
              ))
    ))


(define (dataflow model)
  (let* ((f         (alist-ref 'f model)) 
	 (initial   (alist-ref 'initial model))
	 (input     (map car initial)))
    (print "input = " input)
    (print "dataflow f: "  (dataflow (construct f) input))
    ))


(define (codegen name model #!key (language 'scheme) (random #f) (pre #t) (post #t) (solver 'rkdp))
  (let* ((f         (alist-ref 'f model)) 
	 (initial   (alist-ref 'initial model)))
    (case language
      ((octave Octave) (codegen/Octave name (construct f) initial: initial pre: pre solver: solver))
      ((scheme Scheme) (codegen/scheme name (construct f) initial: initial pre: pre solver: solver))
      ((ML)     (codegen/ML name (construct f) initial: initial random: random pre: pre post: post solver: solver)))
    ))


;; Parameters common to all spiking regimes
(define theta 30.0)
(define Isyn  10.0)
(define gL 0.1)
(define vL 0.0)
(define  C 10.0)
(define theta 20.0)
(define Vreset 10.0)
(define tau_rp 5.0)

;; State initial values
(define t 0.0)
(define V 0.0)
(define t_rpend 0.0)

(define model
  (LeakyIF:construct
   ;; t V gL vL Isyn C theta Vreset t_rpend tau_rp
   t V gL vL Isyn C theta Vreset t_rpend tau_rp))


  
(with-output-to-file "LeakyIF_solver.sml"
  (lambda () 
    (codegen 'LeakyIF model language: 'ML solver: 'rkdp random: #f  post: #t)))