#ifndef _define_h_
#define _define_h_
#include<iostream.h>
#include<fstream.h>
#include<math.h>
#include <vector>
#include<string>
#include<cstring>
#include<sstream>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "../../Neuron.h"
#include "../../Synapse.h"
#include "../../rk4.h"
#include "../../Insynapse/Insynapse.h"
#include "../../Insynapse/Insynapse.cc"
#include "../../HHneuron/HHneuron.h"
//#include "../../QuadraticNeuron/QuadraticNeuron.h"
//#include "../../Minneuron/Minneuron.h"
#include "../../IAcurrent/IA.h"
#include "../../IMcurrent/IM.h"
#include "../../Ionchannel.h"
#include "../../IAHPcurrent/IAHP.h"
#include "../../ILcalcium/IL.h"
#include "../../GABAAsynapse/GABAAsynapse.h"
#include "../../GABABsynapse/GABABsynapse.h"
#include "../../Henrysynapse/Henrysynapse.h"
#include "../../AMPAsynapse/AMPAsynapse.h"
#include "../../TwoDsynapse/TwoDsynapse.h"
#include "../../NMDAsynapse/NMDAsynapse.h"
#include "../../Ihcurrent/Ih.h"
#include "../../matrix.h"
#include "../../Calciumchannel.h"
#include "../../ITcalcium/IT.h"
#include "../../randomc.h"
#include "../../PDrule/PD.h"
#include "../../Couplingcurrent/Couplingcurrent.h"
#include <time.h>
#include <algorithm>
#include <vector>
#define N 4
#ifndef _NO_NAMESPACE
using namespace std;
using namespace math;
#define STD std
#else
#define STD
#endif
#ifndef _NO_TEMPLATE
typedef matrix < double >Matrix;
#else
typedef matrix Matrix;
#endif
/*
Function Spikecount : Stores the times, in the dynamics when v>vpeak..
Basically stores the ISI for Izhikevich type models where u have resetting of the function
when it exceeds certian peak.
Function ISI, computes the average interspike interval from the data generated out of spikecount
Function Max_vect, Min_vect computes the max and min from vector
Function Entropy : Divides the entire range of frequencies, into N bins, and computes the
probability of each bin. From that determines the std informatio entropy. The output is the
Relative entropy which is the Entropy computed divided by the max entropy=log(Nbins)
Function Complexity: Takes in the computed Relative Entropy and spits out, the Complexity
measure: Which characterises the amount of order and disorder in the network.
*/
string
IntToString (int num)
{
ostringstream myStream; //creates an ostringstream object
myStream << num << flush;
/*
* outputs the number into the string stream and then flushes
* the buffer (makes sure the output is put into the stream)
*/
return (myStream.str ()); //returns the string form of the stringstream object
}
string DoubleToStdStr(const double d)
{
std::ostringstream ostr;
ostr << d;
return ostr.str();
}
void
Spikecount (double *v, double *tim, double threshold, int samplesize,
vector < double >&y)
{
if (*tim > 200 && *v > threshold && y.size () < samplesize)
{
y.push_back (*tim);
}
}
double
ISI (vector < double >&y)
{
vector < double >s;
double average, sum;
sum = 0;
for (int i = 1; i < y.size (); i++)
s.push_back (y[i] - y[i - 1]);
for (int i = 0; i < s.size (); i++)
sum = sum + s[i];
average = sum / (s.size ()-1);
return (average);
}
double
max_vect (vector < double >&y)
{
double val = 0;
for (int i = 0; i < y.size (); i++)
{
val = max (y[i], val);
}
return (val);
}
double
min_vect (vector < double >&y)
{
double val = 0;
for (int i = 0; i < y.size (); i++)
{
if (val != 0)
val = min (y[i], val);
}
return (val);
}
long double
Entropy (vector < double >&y, double binsize)
{
int Nbins;
long double En = 0;
sort (y.begin (), y.end ());
double min = y[0];
double max = y[y.size () - 1];
double sum = 0;
if (max > 100)
Nbins = int ((max) / binsize) + 1;
else
Nbins = int (100 / binsize) + 1;
vector < double >hist (Nbins, 0);
vector < double >prob (Nbins, 0);
for (int i = 1; i <= Nbins; i++)
{
for (int j = 0; j < y.size (); j++)
{
if ((i - 1) * binsize <= y[j] && y[j] < i * binsize)
{
hist[i - 1] = hist[i - 1] + 1;
}
}
}
for (int i = 0; i < Nbins; i++)
sum += hist[i];
for (int i = 0; i < hist.size (); i++)
prob[i] = hist[i] / sum;
for (int i = 0; i < Nbins; i++)
{
if (prob[i] != 0)
{
En = En + prob[i] * log10l (1. / prob[i]);
}
}
long double MaxEn = log10l (double (Nbins));
En = En / MaxEn;
return (En);
hist.clear ();
prob.clear ();
}
double
Complexity (double En)
{
double Com;
Com = En * (1 - En);
return (Com);
}
double
Spikepresent (double *tim, double *v, double threshold, double spikewidth,
double vpeak, int &boolean, double *spiketrigger)
{
double spike = 0;
if (*tim > 0 && *v > threshold && boolean == 0)
{
boolean = 1;
spiketrigger[0] = *tim;
}
if (*tim >= spiketrigger[0] && *tim <= spiketrigger[0] + spikewidth)
{
spike = 1;
}
if (*v > vpeak)
boolean = 0;
return (spike);
}
void
Fread (char *y, Matrix & M, int Ntotal)
{
ifstream fread;
vector < double >read (Ntotal, 0);
int row = 0;
fread.open (y, ios::in);
for (int i = 0; i < Ntotal; i++)
fread >> read[i];
if (!fread)
{
cerr << "error in opening" << endl;
exit (0);
}
while (!fread.eof ())
{
for (int j = 0; j < Ntotal; j++)
{
M (row, j) = read[j];
}
for (int j = 0; j < Ntotal; j++)
fread >> read[j];
row += 1;
}
}
/****Function to compute spike times for any given spike train****/
void
spiketimes (double *tim, double *v, double Th, int &bol, vector < double >&y)
{
if (*v > Th && bol == 0)
y.push_back (*tim);
if (*v > Th)
bol = 1;
else
bol = 0;
}
double
Rule (double x, double a, double b, double c, double d,double strength)
{
double out, outmax, outmax1, ap;
outmax = pow (a / b, a) * exp (-a);
outmax1 = pow (c / d, c) * exp (-c);
if (x > 0)
out = strength * (pow (x, 1. * a) * exp (-b * fabs (x))) / outmax;
else
{
if (fmod (c, 2.0) == 0)
{
out = -strength * (pow (x, c) * exp (-d * fabs (x))) / (1 * outmax1);
}
else
{
out = strength * (pow (x, c) * exp (-d * fabs (x))) / (1 * outmax1);
}
}
return out;
}
void
DeltaR (vector < double >In, vector < double >Out, vector < double >&Delta,
vector < double >&result)
{
double diff;
for (int i = 0; i < In.size (); i++)
{
diff = In[i] - Out[i];
//cout<<In[i]<<" "<<Out[i]<<" "<<diff<<endl;
Delta.push_back (diff);
}
for (int i = 0; i < Delta.size (); i++)
result.push_back (Rule (Delta[i], 5.0, 1.0, 5.0, 1.0,.5));
}
double
box_muller (double m, double s, double &x1, double &x2, double &w) /* normal random variate generator */
{
/* mean m, standard deviation s */
double y1;
static float y2;
static int use_last = 0;
if (use_last) /* use value from previous call */
{
y1 = y2;
use_last = 0;
}
else
{
w = sqrt ((-2.0 * log (w)) / w);
y1 = x1 * w;
y2 = x2 * w;
use_last = 1;
}
return (m + y1 * s);
}
double
Average (vector < double >y)
{
double result;
double sum = 0;
int k = y.size ();
if (k != 0)
{
for (int i = 0; i < k; i++)
sum = sum + y[i];
result = sum / k;
}
else
result = 0.0;
return (result);
}
void
Raster (ofstream & Outfile, double &tim, int &bul, double &volt,
double threshold, int &i)
{
//rastercomputation
//takes in the outfile file, tim, bul: a boolean number, volt the neurons spiking, the threshold to detect spike, and spits out the time of neuron "i" firing
bul = 0;
if (volt > threshold && bul == 0)
{
Outfile << tim << " " << i + 1 << endl;
bul = 1;
}
else;
}
void
alpha_syn (double *tim, double g, double alph, vector < double >&Tspk,
vector < double >&sum, int &i)
{
double *x;
x = new double[1000];
for (int j = 0; j < Tspk.size (); j++)
{
x[j] = *tim - Tspk[j];
if (x[j] > 0)
{
sum[i] = sum[i] + g * x[j] * alph * alph * exp (-alph * x[j]);
}
}
delete x;
}
double
Maprule (vector < double >Tpre, vector < double >Tpost, double a, double b,
double c, double d)
{
double out = 0;
double diff;
vector < double >Result;
int sizeTpost = Tpost.size ();
if (sizeTpost == 1)
{
for (int j = 0; j < Tpre.size (); j++)
{
diff = Tpost[0] - Tpre[j];
Result.push_back (diff);
}
}
if (sizeTpost==2||sizeTpost>2)
{
for (int j=0;j<Tpre.size();j++)
{
if (Tpre[j]<Tpost[0])
diff=Tpost[0]-Tpre[j];
else if(Tpre[j]>Tpost[1])
diff=Tpost[1]-Tpre[j];
else if (Tpre[j]>=Tpost[0] && Tpre[j]<=Tpost[1])
{
diff=Tpost[0]-Tpre[j];
Result.push_back(diff);
diff=Tpost[1]-Tpre[j];
}
else;
Result.push_back(diff);
}
}
/* double wold = 1000;
if (sizeTpost >= 2)
{
for (int j = 1; j < Tpre.size (); j++)
{
wold=1000;
for (int i = 0; i < sizeTpost; i++)
{
double u = Tpost[i] - Tpre[j];
// cout << " " << sizeTpost << " " << Tpre.
// size () << " " << u << endl;
if (fabs (u) >= 3 && fabs (u) < fabs (wold))
wold=u;
else;
}
diff = wold;
// cout << diff << endl;
Result.push_back (diff);
}
}
*/
for (int j = 0; j < Result.size (); j++)
out += Rule (Result[j], a, b, c, d,.5);
//out+=Rule(Result[j],3.0,.5,5.0,.25);
return (out);
}
void
Pulsedetect (double *tim, vector < double >&Tspike, int &count, double &pulse)
{
for (int j = 0; j < Tspike.size (); j++)
{
if (*tim > Tspike[count] + .5)
count += 1;
if (*tim > Tspike[count] && *tim < Tspike[count] + .5)
pulse = 1;
else
pulse = 0;
}
}
double Kuramoto(vector<double> &Tpre,vector<double> &Tpost,double period)
{
double phaseterm, costerm, sinterm, cossum = 0, sinsum = 0, kura;
double size=min(Tpre.size(),Tpost.size());
for (int i = 0; i < size; i++)
{
phaseterm = fmod (fabs (Tpost[i] - Tpre[i]), period);
costerm = cos (phaseterm);
sinterm = sin (phaseterm);
cossum += costerm;
sinsum += sinterm;
}
kura = sqrt (cossum * cossum + sinsum * sinsum) / size;
return(kura);
}
void
SlowRise (double &tim, double &volt, double &diff, double tauR,double &ref_time, double &out)
{
double bul_syn;
if (volt > 0 && diff > 0)
bul_syn = 1;
else
bul_syn = 0;
if (bul_syn == 1)
{
ref_time = tim;
}
if (tim >= ref_time && tim < ref_time + tauR)
out = 1;
else
out = 0;
}
#endif