#ifndef _define_h_
#define _define_h_
#include<iostream.h>
#include<iomanip.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 "../../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.
*/
double MIN(double a, double b)
{
double Temp;
if (a<=b)
Temp=a;
else Temp=b;
return Temp;
}
double MAX(double a, double b)
{
double Temp;
if (a>=b)
Temp=a;
else Temp=b;
return Temp;
}
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();
}
double
ISI (vector < double >&y)
{
vector < double >s;
double average, sum;
int size;
size=y.size();
sum = 0;
if (size>50)
{
for (unsigned int i = size-50; i < y.size (); i++)
s.push_back (y[i] - y[i - 1]);
for (unsigned int i = 0; i < s.size (); i++)
sum = sum + s[i];
average = sum / (s.size ());
}
else
{
cout<<"insufficient spike increase duration of time or the neuron did not spike "<<endl;
average=0;
}
return (average);
}
double
max_vect (vector < double >&y)
{
double val = 0;
for (unsigned int i = 0; i < y.size (); i++)
{
val = max (y[i], val);
}
return (val);
}
double
min_vect (vector < double >&y)
{
double val = 0;
for (unsigned int i = 0; i < y.size (); i++)
{
if (val != 0)
val = min (y[i], val);
}
return (val);
}
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;
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;
}
double LinearRule(double x,double strength)
{
double out;
if (x>=-10 && x<=10)
out=0.1*x;
if (x>=-25 && x<-10)
out=-(1.0/15)*x-5.0/3.0;
if (x>10 && x<=25)
out=-(1.0/15)*x-5.0/3.0;
if (x>25||x<-25)
out=0.0;
out=out*strength;
return(out);
}
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 << 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 (unsigned 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 (unsigned int j = 0; j < Tpre.size (); j++)
{
diff = Tpost[0] - Tpre[j];
Result.push_back (diff);
}
}
if (sizeTpost==2||sizeTpost>2)
{
for (unsigned 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 (unsigned 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 (unsigned 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 = size-300; i < size; i++)
{
phaseterm = fmod (fabs (Tpost[i] - Tpre[i]), period);
costerm = cos (phaseterm);
sinterm = sin (phaseterm);
cossum += costerm;
sinsum += sinterm;
//cout<<phaseterm<<" "<<sinsum<<" "<<cossum<<endl;
}
kura = sqrt (cossum * cossum + sinsum * sinsum) / 300;
return(kura);
}
void
SlowRise (double &tim,double Threshold, double &volt, int &bul, double tauR,double &ref_time, double &out)
{
if (volt > Threshold && bul==0)
{ ref_time=tim;
// cout<<ref_time<<endl;
}
if (volt>Threshold)
bul = 1;
else bul=0;
if (tim >= ref_time && tim < ref_time + tauR)
{out = 1.0;}
else
out = 0.0;
}
double Phase_Compute(double pre,double post,double period)
{
double phase,inter;
inter=fabs((pre-post));
phase=fmod(inter,period);
//phase=(pre-post)-floor(inter)*period;
return(phase);
}
#endif