;;
;; Simple Model of Spiking Neurons, Izhikevich  E
;;
;; IEEE Transactions on Neural Networks (2003) 14:1569- 1572
;; 
;; 

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

(define (Izhikevich03:construct t V U k1 k2 k3 theta a b c d spike tspike Isyn )

  `((f . ,(let* ((subthreshold-eq         (ODE h t 
					       (V  (+ (* k1 V V)  (* k2 V) k3 (neg U) Isyn))
					       (U  (* a (- (* b V) U)))))

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

		 (tspike-set              (ASSIGN (tspike (if spike t tspike))))

		 (subthreshold-regime     (SEQUENCE subthreshold-eq (SEQUENCE threshold-detect tspike-set)))
		 
		 (refractory-eq           (ASSIGN (t t) 
						  (spike false)
						  (tspike tspike)
						  (V c)
						  (U (+ U d))))
		 )

	    (TRANSIENT subthreshold-regime refractory-eq 'spike)

	    ))

    (initial .
	     ((h              1e-2)
	     (t               ,t)
	     (V               ,V)
	     (U               ,U)
	     (k1              ,k1)
	     (k2              ,k2)
	     (k3              ,k3)
	     (theta           ,theta)
	     (a               ,a)
	     (b               ,b)
	     (c               ,c)
	     (d               ,d)
	     (spike           ,spike)
	     (tspike          ,tspike)
	     (Isyn            ,Isyn)
	     ))
    ))


(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) (pre #t) (post #t) (solver 'rk3))
  (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 pre: pre post: post solver: solver)))
    ))


;; Parameters common to all spiking regimes
(define k1    0.04)
(define k2    5.0)
(define k3    140.0)
(define theta 30.0)
(define Isyn  10.0)

;; State initial values
(define t 0.0)
(define V -65.0)


;; Regular spiking (RS) parameters
(define RS 
  (let ((a 0.02)
	(b 0.2)
	(c -65.0)
	(d 8.0))
    (Izhikevich03:construct
     ;; t V U k1 k2 k3 theta a b c d spike tspike Isyn 
     t V (* b V) k1 k2 k3 theta a b c d #f t Isyn)))


;; Intrinsically bursting (IB) parameters
(define IB 
  (let ((a 0.02)
	(b 0.2)
	(c -55.0)
	(d 4.0))
    (Izhikevich03:construct
     ;; t V U k1 k2 k3 theta a b c d spike tspike Isyn 
     t V (* b V) k1 k2 k3 theta a b c d #f t Isyn)))


;; Chattering (CH) parameters
(define CH
  (let ((a 0.02)
	(b 0.2)
	(c -50.0)
	(d 2.0))
    (Izhikevich03:construct
     ;; t V U k1 k2 k3 theta a b c d spike tspike Isyn 
     t V (* b V) k1 k2 k3 theta a b c d #f t Isyn)))


;; Fast-spiking (FS) parameters
(define FS
  (let ((a 0.1)
	(b 0.2)
	(c -65.0)
	(d 2.0))
    (Izhikevich03:construct
     ;; t V U k1 k2 k3 theta a b c d spike tspike Isyn 
     t V (* b V) k1 k2 k3 theta a b c d #f t Isyn)))



;; Thalamo-cortical (TC) parameters
(define TC
  (let ((a 0.02)
	(b 0.25)
	(c -65.0)
	(d 0.05))
    (Izhikevich03:construct
     ;; t V U k1 k2 k3 theta a b c d spike tspike Isyn 
     t V (* b V) k1 k2 k3 theta a b c d #f t Isyn)))


;; Resonator (RZ) parameters
(define RZ
  (let ((a 0.1)
	(b 0.26)
	(c -65.0)
	(d 2.0))
    (Izhikevich03:construct
     ;; t V U k1 k2 k3 theta a b c d spike tspike Isyn 
     t V (* b V) k1 k2 k3 theta a b c d #f t Isyn)))



;; Low-threshold spiking (LTS) parameters
(define LTS
  (let ((a 0.02)
	(b 0.25)
	(c -65.0)
	(d 2.0))
    (Izhikevich03:construct
     ;; t V U k1 k2 k3 theta a b c d spike tspike Isyn 
     t V (* b V) k1 k2 k3 theta a b c d #f t Isyn)))


(define models `((RS  . ,RS) (IB . ,IB) (CH . ,CH) 
		 (FS . ,FS)  (RZ . ,RZ) (LTS . ,LTS) ))

  
(with-output-to-file "Izhikevich03_solver.m"
  (lambda () 
    (let recur ((models models) (pre #t))
      (if (pair? models)
	  (let* ((name.model (car models))
		 (name  (car name.model))
		 (model (cdr name.model)))
	    (codegen name model language: 'octave pre: pre solver: 'lsode)
	    (recur (cdr models) #f)
	    )))
    ))

(with-output-to-file "Izhikevich03_solver.scm"
  (lambda () 
    (let recur ((models models) (pre #t))
      (if (pair? models)
	  (let* ((name.model (car models))
		 (name  (car name.model))
		 (model (cdr name.model)))
	    (codegen name model language: 'scheme pre: pre)
	    (recur (cdr models) #f)
	    )))
    ))

  
(with-output-to-file "Izhikevich03_solver.sml"
  (lambda () 
    (let recur ((models models) (pre #t))
      (if (pair? models)
	  (let* ((name.model (car models))
		 (name  (car name.model))
		 (model (cdr name.model)))
	    (codegen name model language: 'ML pre: pre post: (null? (cdr models)))
	    (recur (cdr models) #f)
	    )))
    ))