#!perl -w

#
# Created by Paul Kulchenko (paul@kulchenko.com); December 20, 2005
#
# This script is based on the model proposed by Eugene M. Izhikevich and
# described in his article "Simple Model of Spiking Neurons" available
# here: http://www.nsi.edu/users/izhikevich/publications/spikes.htm
# Parameters are taken from his book "Dynamical Systems in Neuroscience:
# The Geometry of Excitability and Bursting" published on his website and
# available here: http://www.nsi.edu/users/izhikevich/publications/dsn.pdf
#

use strict;
use List::Util;
use PDL;
use PDL::IO::Pic;

$| = 1;

my $n = 1000;             # total number of neurons in the model
my $n_ex = int($n * 0.8); # excitatory : inhibitory = 4:1
my $n_in = $n - $n_ex;

my $r_ex = random($n_ex); 

# model parameters for excitatory (regular spiking (RS) and chattering (CH)) 
# and inhibitory (fast spiking (FS)) neurons
#        text                       a               b    c    d    C    k   vr   vt    vp 
#  (map ['(RS) regular spiking', 0.03,             -2, -50, 100, 100, 0.7, -60, -40,   35,
#  (map ['(CH) chattering',      0.03,             -2, -40, 150,  50, 1.5, -60, -40,   25,
#  (map ['(FS) fast spiking',     0.2, b2(-55, 0.025), -45,   0,  20,   1, -55, -40,   25,

my $a =  (0.03 * ones($n_ex)) ->append(0.2 * ones($n_in)); 
my $b =  (-2   * ones($n_ex)) ->append(0 * ones($n_in));     #  RS --  CH
my $c =  (-50  + 10 * $r_ex)  ->append(-45 * ones($n_in));   # -50 -- -40
my $d =  (100  + 50 * $r_ex)  ->append(0 * ones($n_in));     # 100 -- 150
my $k =  (0.7  + 0.8 * $r_ex) ->append(1 * ones($n_in));     # 0.7 -- 1.5
my $C =  (100  - 50 * $r_ex)  ->append(20 * ones($n_in));    # 100 --  50
my $vp = (35   - 10 * $r_ex)  ->append(25 * ones($n_in));    #  35 --  25
my $vr = (-60  + 0 * $r_ex)   ->append(-55 * ones($n_in)); 
my $vt = (-40  * ones($n));

my $t_ex = grandom($n_ex, $n);
my $t_in = grandom($n_in, $n);
my $s = (0.5 * (($t_ex >= 0) * $t_ex + ($t_ex < 0) * (-$t_ex)))->glue(0 
     => -(($t_in >= 0) * $t_in + ($t_in < 0) * (-$t_in)));

my $v = $vr->copy;
my $u = zeroes($n);

my $firings = null;

my $realtime = time;
my $totalfired = 0;

# calculate two (possibly different) time steps for excitatory and inhibitory neurons
my($tau_ex, $tau_in) = (1, 0.2); # assumption: tau_ex > tau_in
my $tau = List::Util::min($tau_ex, $tau_in);
my $tau_reset = List::Util::max($tau_ex, $tau_in) / $tau;
my $tau_diff = 0;

my $i = zeroes($n);

# separate each matrix into two slices: one for excitatory and one for inhibitory neurons
my($v_ex, $k_ex, $vr_ex, $vt_ex, $u_ex, $i_ex, $a_ex, $b_ex, $C_ex) = 
  map {$_->slice(join ':', 0, $n_ex-1)} ($v, $k, $vr, $vt, $u, $i, $a, $b, $C);
my($v_in, $k_in, $vr_in, $vt_in, $u_in, $i_in, $a_in, $b_in, $C_in) = 
  map {$_->slice(join ':', $n_ex, -1)} ($v, $k, $vr, $vt, $u, $i, $a, $b, $C);

my $simultime = 1000; # time in ms
for (my $t = 0; $t < $simultime; $t += $tau) {
  my $f = $v >= $vp;                         # find neurons that fired
  my $fired = $f->which;                     # ...and their indeces
  $firings = $firings->append($f)->copy;     # store fired neurons for diagramming
  $v->nslice($fired) .= $c->nslice($fired);  # membrane potential reset
  $u->nslice($fired) += $d->nslice($fired);  # membrane recovery variable reset

  $i += 10 * sumover($s->dice($fired, 'X')); # calculate synaptic currents
                                             # based on weights from fired neurons

  my($ex, $in) = (!$tau_diff++, 1); 
  $tau_diff = 0 if $tau_diff >= $tau_reset;

  if ($ex) { # time to calculate excitatory neurons
    $v_ex += $tau_ex * ($k_ex * ($v_ex - $vr_ex) * ($v_ex - $vt_ex) - $u_ex + $i_ex) / $C_ex;
    $u_ex += $tau_ex * $a_ex * ($b_ex * ($v_ex - $vr_ex) - $u_ex);
    $i_ex .= 300 * random($n_ex);
  }
  if ($in) { # time to calculate inhibitory neurons
    $v_in += $tau_in * ($k_in * ($v_in - $vr_in) * ($v_in - $vt_in) - $u_in + $i_in) / $C_in;
    $u_in += $tau_in * $a_in * (((($v_in >= -55) * 0.025 * ($v_in - $vr_in)**2)) * ($v_in - $vr_in) - $u_in);
    $i_in .= 300 * random($n_in);
  }

  $totalfired += sum($f);

  print int($t), "ms\r";# unless $t % 10;
}

print "Simulation of $n neurons (${simultime}ms of model time) took ", time-$realtime, " seconds, with ", 
      $totalfired*1000/$simultime/$n, " spikes per neuron per second\n";

# turn firings into black on white instead of white on black
$firings = ! $firings->reshape($n, ($firings->dims)[0]/$n)->transpose;
$firings->wpic('spikes.gif');