#!/usr/bin/perl

print "\n\n";
for ($i=1;$i<=80;$i++) {print "#"}; print "\n";
print "This script calculates the synchrony of a given set of neurons following the
definition of synchrony by Hansel et al. 1998. Furthermore, the script performs
this task for all data generated by the corresponding SLI script, i.e. for each
couping weight and neuron model the resulting synchrony is calculated.

At the end a Gnuplot figure is generated.\n

Needs:
- Perl
- Gnuplot

AH 09\n";
for ($i=1;$i<=80;$i++) {print "#"}; print "\n";


$nr = 128;              # number of neurons
$timestart = 5000 ;     # start position
$timesteps = 5000 ;     # number of time steps
$strength_inc = 0.2;    # weight steps

for ($sim = 1; $sim <=2 ; $sim++) {     # 1 = canonical, 2 = grid constrained
 $cnt_strength = 0;
 for ($strength=0;$strength<=5.0;$strength+=$strength_inc) {

  if ($sim == 1)
   { $filename =  "voltmeter-Canon-".$strength."-129-0.dat";}
   else
   { $filename =  "voltmeter-Grid-".$strength."-129-0.dat";}

  print "open ",$filename,"\n";

  open (ifile,$filename);               # store data in memory
  $line_nr = 0;
  for $line (<ifile>) {
    @data_l = split (/\s+/,$line);
    $Data[$line_nr][0] = $data_l[0];
    $Data[$line_nr][1] = $data_l[1];
    $Data[$line_nr][2] = $data_l[2];
    $line_nr++;
  }
  $max_line = $line_nr;                 # remember number of data points
  close ifile;


  ### Calculate A_n(t)
  $ii = 1;
  $Vi = $neuron;
  $data_set=0;

  $sum_Vi = 0.0;
  $sum_Ai = 0.0;
  $sum_sqAi = 0.0;

  for ($cnt_array=0;$cnt_array<$max_line;$cnt_array++) {
   if (($data_set >= $timestart) && ($data_set <= $timestart+$timesteps)) {
    if ($ii == 1) {   # new dataset time t
     #print  $sum_Vi/$nr,"\n";
     $sum_Ai += $sum_Vi/$nr;
     $sum_sqAi += ($sum_Vi/$nr)**2;
     $sum_Vi =0;
    }
    $sum_Vi += $Data[$cnt_array][2];
   }

   $ii++;                               # next neuron
   if ($ii>$nr) {                       # next time step
    $data_set++;
    $ii=1;
   };
  }


  $Delta_N = ($sum_sqAi-($sum_Ai**2/$timesteps))/($timesteps);    ### CAUTION in denominator correctly /($timesteps-1) but to be comparable to Python '1' omitted

  $Delta = 0.0;
  for ($neuron = 1; $neuron <= $nr; $neuron++) {
   $ii = 1;
   $Vi = $neuron;
   $data_set=0;

   $sum_Vi = 0.0;
   $sum_sqVi = 0.0;

   for ($cnt_array=0;$cnt_array<$max_line;$cnt_array++)  {
    if (($data_set >= $timestart) && ($data_set < $timestart+$timesteps)) {
     if ($Vi==$ii) {
      $sum_Vi += $Data[$cnt_array][2];
      $sum_sqVi += $Data[$cnt_array][2]**2;
     }
    }
    $ii++;
    if ($ii>$nr) {
     $data_set++;
     $ii=1;
    };
   }

   $Delta += ($sum_sqVi-($sum_Vi**2/$timesteps))/($timesteps); ### CAUTION in denominator correctly /($timesteps-1) but to be comparable to Python changed

  }

  $Delta = $Delta/$nr;
  ######print "\nDelta = ",$Delta,"\n";
  #####print "\nSigma = ",$Delta_N/$Delta,"\n";

  if ($Delta!=0) {
   print $Delta_N/$Delta,"\n";
   $Plot_Data[$sim][$cnt_strength] = $Delta_N/$Delta;
   }
   else {
   print "Error\n";
   $Plot_Data[$sim][$cnt_strength] = 99;
   }
  $cnt_strength++;
 }
 $max_cnt_strength=$cnt_strength;
}


$filename =  ">ArtificialSynchrony.dat";                # Write results to data file
open (ofile,$filename);
print ofile "#Strength [pA]\t#Sigma_Canon\t#Sigma_Grid \n";
for ($cnt_strength=0;$cnt_strength<$max_cnt_strength;$cnt_strength++) {
 print ofile $cnt_strength*$strength_inc,"\t$Plot_Data[1][$cnt_strength]\t$Plot_Data[2][$cnt_strength]\n";
}
close ofile;


$filename =  ">ArtificialSynchrony.gnu";                # Generate Gnuplot file to print final data
open (ofile,$filename);
print ofile 'set terminal postscript enhanced color solid "Times-Roman" 20
set encoding iso

set title "Network Synchrony of 128 IAF neurons, {/Symbol g} = 0.5, h = 2^{-5} ms" font "Times-Roman, 25"
set xlabel "Coupling strength " font "Times-Roman, 25"
set ylabel "Synchrony {/Symbol S}" font "Times-Roman, 25"

set xtics font "Times-Roman, 25"
set ytics font "Times-Roman, 25"


set yrange [0:1]
set xrange [0:1]
plot \'ArtificialSynchrony.dat\' u ($1/3.90625):3 w l title "" lt -1 lw 2, \
     \'ArtificialSynchrony.dat\' u ($1/3.90625):3 w p title "alpha" lt -1 pt 8 ps 2, \
     \'ArtificialSynchrony.dat\' u ($1/3.90625):2 w l title "" lt -1 lw 2, \
     \'ArtificialSynchrony.dat\' u ($1/3.90625):2 w p title "alpha canon" lt -1 pt 6 ps 2';
close ofile;


system ("gnuplot ArtificialSynchrony.gnu > ArtificialSynchrony.ps && gv ArtificialSynchrony.ps");       # Create plot and show it