#!perl -w

#
# Created by Paul Kulchenko (paul@kulchenko.com); November 22, 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 qw(min max);
use Image::Magick;

my @types = (
#       text                                    a              b                c             d    C     k   vr   vt              vp    v     u         x   tau    max  i
  (map ['(RS) regular spiking',              0.03,            -2,             -50,          100, 100,  0.7, -60, -40,             35,   0,    0,        0,    1,  1000, $_, i1(1000, $_)], 200, 100, 70, 55),
  (map ['(IB) intrinsically bursting',       0.01,             5,             -56,          130, 150,  1.2, -75, -45,             50,   0,    0,        0,    1,  1000, $_, i1(1000, $_)], 550, 500, 370, 300),
  (map ['(CH) chattering',                   0.03,            -2,             -40,          150,  50,  1.5, -60, -40,             25,   0,    0,        0,    1,   200, $_, i1(200, $_)], 600, 400, 300, 200),
  (map ['(NS) neostriatal spiny',            0.01,            -20,            -55,          150,  50,    1, -80, -25,             40,   0,    0,        0,    1,   500, $_, i2(500, $_)], 640, 520, 420, 400),
  (map ['(TCt) talamocortical; tone',        0.01, b1(-65, 15, 0),  c1(-60, -0.1),           10, 200,  1.6, -60, -50,   vp1(35, 0.1),   0,    0,        0,    1,  1000, $_, i1(1000, $_)], 200, 100, 50, 30),
  (map ['(TCb) talamocortical; burst',       0.01, b1(-65, 15, 0),  c1(-60, -0.1),           10, 200,  1.6, -60, -50,   vp1(35, 0.1), -85, -220,        0,    1,  1000, $_, i1(1000, $_)], 200, 100, 50, 30),
  (map ['(RTNt) reticular thalamic; tone',  0.015, b1(-65, 10, 2),            -55,           50,  40, 0.25, -65, -45,              0,   0,    0,        0,    1,   600, $_, i1(600, $_)], 110, 75, 50, 40),
  (map ['(RTNb) reticular thalamic; burst', 0.015, b1(-65, 10, 2),            -55,           50,  40, 0.25, -65, -45,              0, -85, -220,        0,    1,   600, $_, i1(600, $_)], 110, 75, 50, 40),
  (map ['(FS) fast spiking',                  0.2, b2(-55, 0.025),            -45,            0,  20,    1, -55, -40,             25,   0,    0,        0, 0.05,   130, $_, i1(130, $_)], 400, 200, 100, 73.2),
  (map ['(LS) late spiking; +noise',         0.15,              6,            -40,          100,  30,  0.4, -56, -37,             30,   0,    0, vd1(-60),  0.1,   800, $_, i2r(800, $_, 1/8, 7/8, 0.1, 2*$_)], 200, 120, 70, 60),
  (map ['(LTS) low threshold spiking; mod',  0.03,              8,            -45,            0, 100,    1, -56, -42,  vp1(40, -0.1),   0,    0,        0, 0.05,   500, $_, i1(500, $_, 0.04)], 275, 200, 125, 100),
  (map ['(TI) thalamic interneuron',         0.05,              7,  c1(-65, 0.08),  d1(50, 530),  20,  0.5, -60, -50, vp1(20, -0.08),   0,    0,        0,  0.2,   400, $_, i2(400, $_, 0.07, 0.8)], 250, 200, 100, 50),
# almost working (except with large currents)...
  (map ['(ST) stellate cell',                0.01,             15,            -50,            0, 200, 0.75, -60, -45,             30,   0,    0,        0,    1,  2000, $_, i2(2000, $_)], -500, 100, 165, 200),
  (map ['(STn) stellate cell; +noise',       0.01,             15,            -50,            0, 200, 0.75, -60, -45,             30,   0,    0,        0,    1,  5000, $_, i2r(5000, $_)], 173, 170, 167, 165),
# not working...
#  (map ['(LTS) low threshold spiking',       0.03,              8,   c1(-53, 0.4),           20, 100,    1, -56, -42,  vp1(40, -0.1),   0,    0,        0, 0.05,   50, $_, i1(50, $_, 0.25)], 275, 200, 125, 100),
# playing with...
#  (map ['(RSp) regular spiking',             0.03,            -2,             -50,          100, 100,  0.7, -60, -40,             35,   0,    0,        0,    5,  1000, $_, i1(1000, $_)], 200, 100, 70, 55),
#  (map ['(RSp) regular spiking',             0.03,            -2,             -50,          100, 100,  0.7, -60, -40,             35,   0,    0,        0,    1,  100,  $_, i1(100, $_)], 1000, 800, 500, 100),
#  (map ['(RSp) regular spiking',             0.03,            -2,             -50,          100, 100,  0.7, -60, -40,             35,   0,    0,        0,    1,  1000, $_, i2(1000, $_, 0.1, 0.12)], 200, 100, 70, 55),
#  (map ['(IBp) intrinsically bursting',      0.01,             5,             -56,          130, 150,  1.2, -75, -45,             50,   0,    0,        0,    1,  1000, $_, i2(1000, $_, 0.1, 0.12)], 550, 500, 370, 300),
);

our $DEBUG = shift if @ARGV && $ARGV[0] eq '-d';

my($show_types) = join '|', map {quotemeta} (@ARGV ? @ARGV : (map {$_->[0]=~/\((.+)\)/ && $1} @types));
my($xsize, $ysize) = (300, 300);
my($tsize, $usize, $isize, $border) = (22, 30, 30, 10);

my $image = Image::Magick->new;

for (grep {$_->[0] =~ /\((?:$show_types)\)/i} @types) {

  my(@parameters) = @$_;
  my($text, $tau, $max, $cur) = @parameters[0, 13, 14, 15];
  my($mv, $ms) = (10, $max/10);

  $max /= 1000; # ms -> s
  $parameters[0] = ($text .= " (${max}s/${tau}ms/${cur}pA)");

  print STDERR "$text\n";

  # remove 'text' and 'I' parameters; they are not used in calculations
  splice(@parameters, 15, 1); shift(@parameters); 
  my($v, $i, $u) = calculate(@parameters);

  my $tile = Image::Magick->new;
  $tile->Set(size=>"${xsize}x${ysize}");
  $tile->ReadImage('xc:white');

  plot($tile, [$border, $tsize, $xsize-2*$border, $ysize-$tsize-$usize-$isize-3*$border], $v, $text, $ms, $mv, $tau);
  plot($tile, [$border, $ysize-$usize-$isize-2*$border, $xsize-2*$border, $usize], $u, (undef) x 4, {stroke => 'green'});
  plot($tile, [$border, $ysize-$isize-$border, $xsize-2*$border, $isize], $i, (undef) x 4, {stroke => 'blue'}); 

  push @$image, $tile;
}

@$image or die "No neuron type selected to work on (@ARGV)\n";

$image->Montage(geometry=>"${xsize}x${ysize}+5+7", @$image > 4 ? (tile=>'4x') : ())
      ->Write('neuron.gif');

sub i1  {my $thr = $_[0] * ($_[2] || 0.1); my $a = $_[1]; sub {$_[0] > $thr ? $a : 0}}
sub i2  {my($max, $thr) = ($_[0] * ($_[3] || 0.8), $_[0] * ($_[2] || 0.2)); my $a = $_[1]; sub {$_[0] > $thr && $_[0] < $max ? $a : 0}}
sub i2r {my($max, $thr) = ($_[0] * ($_[3] || 0.8), $_[0] * ($_[2] || 0.2)); my($a, $tau, $mul, $noise) = (@_[1, 4, 5], 0); sub {$noise += $tau*(2*rand()-1-$noise)/2; $_[0] > $thr && $_[0] < $max ? $a + $mul * $noise : 0}}
sub b1  {my($thr, $a, $b) = @_; sub {$_[1] <= $thr ? $a : $b}}
sub b2  {my($a, $b) = @_; sub{$_[1] < $a ? 0 : $b*($_[1] - $a)**2}}
sub c1  {my($a, $b) = @_; sub{$a+$b*$_[0]}}
sub d1  {my($a, $b) = @_; sub{min($_[0]+$a, $b)-$_[0]}}
sub vp1 {goto &c1}
sub vd1 {my $vd = $_[0]; sub {$vd = 0.01*($_[1]-$vd); 1.2*($vd-$_[1])}}

sub calculate {
  my($a, $b, $c, $d, $C, $k, $vr, $vt, $vp, $v, $u, $x, $tau, $max, $current) = @_;

  $v ||= $vr;
  $u ||= 0;

  $DEBUG and print "# T I u v; spikes are marked as: . (v >= v peak)\n";

  my(@v, @u, @i);
  for (my $t = 0; $t < $max; $t += $tau) {
    my $i = $current->($t);

    $v += $tau * ($k * ($v - $vr) * ($v - $vt) + (ref $x ? $x->($u, $v) : $x) - $u + $i) / $C;
    $u += $tau * $a * ((ref $b ? $b->($u, $v) : $b) * ($v - $vr) - $u);
  
    my $vp = (ref $vp ? $vp->($u, $v) : $vp);
    if ($v >= $vp) {
    $DEBUG and print ". ($v >= $vp)\n";
      push(@v, $vp);
      $v = ref $c ? $c->($u, $v) : $c;
      $u += ref $d ? $d->($u, $v) : $d;
    } else {
      push(@v, $v);
    }
    push(@u, $u);
    push(@i, $i);
    $DEBUG and printf "%.2f %s\n", $t, "$i $v $u";
  }

  return (\@v, \@i, \@u);
}

sub plot {
  my($image, $xy, $f, $text, $xscale, $yscale, $tau, $params) = @_;

  my($xoff, $yoff, $xpx, $ypx) = @$xy;
  my $xrange = [1, scalar @$f];
  my $yrange = [min(@$f), max(@$f)];
  my $scale = '0.8,0.5';
  my($ymax, $ymin, $ybeg, $xbegend);

  my $xratio = (($xrange->[1] - $xrange->[0]) / $xpx) || 1;
  my $yratio = (($yrange->[1] - $yrange->[0]) / $ypx) || 1;

  my $n = 0;
  my($xo, $yo);
  foreach my $y (@$f) {
    my $x = ++$n;

    # if any of variables is out of range, then just reset old variables
    unless ($y >= $yrange->[0] && $y <= $yrange->[1] &&
            $x >= $xrange->[0] && $x <= $xrange->[1]) {
      $xo = $yo = undef;
      next;
    }

    my $xp = int($x / $xratio);
    my $yp = $ypx - int(($y - $yrange->[0])/ $yratio);

    $image->Draw(points => "${\($xo+$xoff)},${\($yo+$yoff)} ${\($xp+$xoff)},${\($yp+$yoff)}", 
                 stroke => 'black', primitive => 'line', antialias => 'false', strokewidth => 0.5, %$params) 
      if defined $xo;

    # ($x_ppem, $y_ppem, $ascender, $descender, $width, $height, $max_advance) = $image->QueryFontMetrics(parameters);

    # mark beginning
    if (!defined $ybeg) {
      my $show = ($ybeg = $y);
      my($width, $height) = ($image->QueryFontMetrics(text => $show, scale => $scale))[4, 5]; # get text width
      $xbegend = $xp+$xoff+$width;
      $image->Annotate(text => $show, y => $yp+$yoff-3, x => $xp+$xoff, scale => $scale, fill => 'red', align => 'Left');
      $image->Draw(points => "${\($xp+$xoff)},${\($yp+$yoff)} ${\($xp+$xoff+1)},${\($yp+$yoff)}", stroke => 'red', primitive => 'circle', strokewidth => 2);
    }

    # mark first max
    if (!defined $ymax && $y != $ybeg && ($y == $yrange->[1])) {
      my $show = sprintf('%.5g', $ymax = $y);
      my($width, $height) = ($image->QueryFontMetrics(text => $show, scale => $scale))[4, 5]; # get text width
      $image->Draw(points => "${\($xp+$xoff)},${\($yp+$yoff)} ${\($xp+$xoff+1)},${\($yp+$yoff)}", stroke => 'red', primitive => 'circle', strokewidth => 2);
      $image->Annotate(text => $show, y => $yp+$yoff+$height-5, x => max($width, $xp+$xoff-3), scale => $scale, fill => 'red', align => 'Right');
    }

    # mark first min
    if (!defined $ymin && $y != $ybeg && ($y == $yrange->[0])) {
      my $show = sprintf('%.5g', $ymin = $y);
      my $width = ($image->QueryFontMetrics(text => $show, scale => $scale))[4]; # get text width
      $image->Draw(points => "${\($xp+$xoff)},${\($yp+$yoff)} ${\($xp+$xoff+1)},${\($yp+$yoff)}", stroke => 'red', primitive => 'circle', strokewidth => 2);
      $image->Annotate(text => $show, y => $yp+$yoff-3, x => max($xbegend+$width, $xp+$xoff-3), scale => $scale, fill => 'red', align => 'Right');
    }

    ($xo, $yo) = ($xp, $yp);
  }

  if ($xscale) {
    $image->Draw(points => "${\($xsize-$xscale/$tau/$xratio-3)},10 ${\($xsize-3)},10",
                 stroke => 'gray', primitive => 'line', antialias => 'false', strokewidth => 1, %$params);
    $image->Annotate(text => "${xscale}ms", y => 8, x => $xsize-30, scale => $scale);
  }
  if ($yscale) {
    $image->Draw(points => "${\($xsize-3)},15 ${\($xsize-3)},${\(15+int($yscale/$yratio))}",
                 stroke => 'gray', primitive => 'line', antialias => 'false', strokewidth => 1, %$params);
    $image->Annotate(text => "${yscale}mV", y => 18, x => $xsize-30, scale => $scale);
  }

  $image->Annotate(text => $text, y => 15, x => 5) if $text;

  return;
}