#include "define.h"
#include "Neur_Class.h"
#include "Corr_Stat.h"
#include <time.h>
int
main(int ac, char **av)
{
time_t start, end1, end2;
double dG, dT, gmutualoldpre, gmutualnewpre, gmutualoldpost,
gmutualnewpost;
int maxsz, minsz;
int maxold=0, maxnew, minold=0, minnew;
time(&start);
int32 seed = time(NULL);
TRandomMersenne rg(seed);
string s1, s2, s3, s4, s5;
if (ac <= 14) {
puts("NN");
puts("endtime");
puts("Iinput external current I_ext Pre neuron");
puts("Baseline current which is kept fixed, to decide the regime Idc post neuron");
puts("self inhibition strength");
puts("mutual inhibition strength");
puts("decay constant for inhibition");
puts("do you want the time series data");
puts("noise level");
puts("1 if Learning present, O otherwise");
puts("output file name");
puts("Self inhibition strength for pre neuron: Heterogeneity through gs");
puts("Tau Rise time");
puts("Spike time output file name");
puts("Iteration Number for Noise Study");
exit(0);
}
int NN = atoi(av[1]);
double endtime = atof(av[2]);
double I_ext = atof(av[3]);
double Idc = atof(av[4]);
// int Synapse_Type = atoi(av[5]);
double gself = atof(av[5]);
double gmutual = atof(av[6]);
double taus = atof(av[7]);
int scale = atoi(av[8]);
double noise = atof(av[9]);
int Learnpresent = atoi(av[10]);
string out = av[11];
double gselfnew = atof(av[12]);
double tauR = atof(av[13]);
int NoiseNum=atoi(av[14]);
double tauD = taus;
double tim=0.0, timestep = epsilon, Threshold = 0;
double gelec=0.0;
vector < double >x1, x2, w, *T;
vector < double >x1p, x2p, wp, *TP;
T = new vector < double >[3000];
TP = new vector < double >[3000];
int *bol, *bolp;
bol = new int [NN + 10];
bolp = new int [NN + 10];
for (int i = 0; i < NN; i++) {
x1.push_back(0.);
x2.push_back(0.);
w.push_back(0.);
x1p.push_back(0.);
x2p.push_back(0.);
wp.push_back(0.);
}
ofstream Outfile, Outfile1, Outfile2;
Outfile.open(out.c_str(), ios::app);
s1 = "Timeseries";
s2 = DoubleToStdStr(I_ext);
s3 = DoubleToStdStr(noise);
s4 = IntToString(Learnpresent);
s1 = s1 + "_" + s2 + "_" + s3 + "_" + s4 + ".dat";
s5 = "Rates";
s5 = s5 + "_" + s2 + "_" + s3 + "_" + s4 + ".dat";
s5 = "RatesNetwork/" + s5;
//Outfile2.open((char *) s5.c_str(), ios: :out);
HHneuron **neurpre, **neurpost;
neurpre = new HHneuron *[NN];
neurpost = new HHneuron *[NN];
for (int i = 0; i < NN; i++) {
neurpre[i] = new HHneuron(.1, 1);
neurpre[i]->parameter[11] = -65;
neurpre[i]->parameter[7] = 0.;
neurpre[i]->x[0] = -65;
neurpost[i] = new HHneuron(.1, 1);
neurpost[i]->parameter[11] = -65;
neurpost[i]->parameter[7] = 0;
neurpost[i]->x[0] = -55.23;
}
Insynapse **Inputpre, **Inputpost;
Inputpre = new Insynapse *[NN];
Inputpost = new Insynapse *[NN];
for (int i = 0; i < NN; i++) {
Inputpre[i] = new Insynapse(neurpre[i], 0.0);
Inputpost[i] = new Insynapse(neurpost[i], 0.0);
}
Couplingcurrent ***coupAB, ***coupBA;
coupAB = new Couplingcurrent **[NN];
coupBA = new Couplingcurrent **[NN];
for (int i = 0; i < NN; i++) {
coupAB[i] = new Couplingcurrent *[NN];
coupBA[i] = new Couplingcurrent *[NN];
for (int j = 0; j < NN; j++) {
coupAB[i][j] =
new Couplingcurrent(neurpre[i], neurpost[j], gelec / NN);
coupBA[i][j] =
new Couplingcurrent(neurpost[i], neurpre[j], gelec / NN);
}
}
TwoDsynapse ***slowself, ***slowselfA, ***slowmutual, ***slowmutualA;
slowself = new TwoDsynapse **[NN];
slowselfA = new TwoDsynapse **[NN];
slowmutual = new TwoDsynapse **[NN];
slowmutualA = new TwoDsynapse **[NN];
for (int i = 0; i < NN; i++) {
slowself[i] = new TwoDsynapse *[NN];
slowselfA[i] = new TwoDsynapse *[NN];
slowmutual[i] = new TwoDsynapse *[NN];
slowmutualA[i] = new TwoDsynapse *[NN];
slowself[i][i] = new TwoDsynapse(neurpre[i], neurpre[i], gselfnew, 0., 0);
slowselfA[i][i] = new TwoDsynapse(neurpost[i], neurpost[i], gself, 0., 0);
slowself[i][i]->parameter[1] = -82;
slowselfA[i][i]->parameter[1] = -82;
slowself[i][i]->parameter[2] = tauR;
slowself[i][i]->parameter[3] = tauD;
slowselfA[i][i]->parameter[2] = tauR;
slowselfA[i][i]->parameter[3] = tauD;
for (int j = 0; j < NN; j++) {
slowmutual[i][j] = new TwoDsynapse(neurpre[i], neurpost[j], 1.0*gmutual / NN, 0., 0);
slowmutualA[i][j] = new TwoDsynapse(neurpost[i], neurpre[j], gmutual / NN, 0., 0);
slowmutual[i][j]->parameter[1] = -82;
slowmutualA[i][j]->parameter[1] = -82;
slowmutual[i][j]->parameter[2] = tauR;
slowmutual[i][j]->parameter[3] = tauD;
slowmutualA[i][j]->parameter[2] = tauR;
slowmutualA[i][j]->parameter[3] = tauD;
}
}
int count = 0;
double *curr_pre, *curr_post, sumprecur=0, sumpostcur=0, Het;
curr_pre = new double[NN + 10];
curr_post = new double[NN + 10];
for (int i = 0; i < NN; i++) {
curr_pre[i] = I_ext + 0.0 * rg.Random();
curr_post[i] = Idc + 0.0 * rg.Random();
sumprecur += curr_pre[i];
sumpostcur += curr_post[i];
}
Het = 100 * (sumprecur - sumpostcur) / (sumprecur + sumpostcur);
/*******************parameters for the TwoDsynapse to work********************/
double *ref_time_pre, *ref_time_post;
double *vnew_pre, *vnew_post, *vold_pre, *vold_post;
double *diff_pre, *diff_post;
int *bol_pre, *bol_post, *bol_postpre, *bol_prepost;
ref_time_pre = new double[NN];
ref_time_post = new double[NN];
vnew_pre = new double[NN];
vnew_post = new double[NN];
vold_pre = new double[NN];
vold_post = new double[NN];
diff_pre = new double[NN];
diff_post = new double[NN];
bol_pre = new int[NN];
bol_post = new int[NN];
bol_prepost = new int[NN];
bol_postpre = new int[NN];
for (int i = 0; i < NN; i++) {
vnew_pre[i] = 0;
vnew_post[i] = 0;
vold_pre[i] = neurpre[i]->x[0];
vold_post[i] = neurpost[i]->x[0];
bol_pre[i] = bol_post[i] = 0;
bol_prepost[i] = bol_postpre[i] = 0;
}
/*****************************************************************************/
gmutualoldpre = gmutualoldpost = gmutual;
double Hetgself=100*(gselfnew-gself)/(gselfnew+gself);
int szpreold=0,szpostold=0;
double sumcurr=0,sumcurrpost=0;
vector<double> PreCurr,PostCurr;
double RuleNoise=0;
while (tim < endtime) {
for (int i = 0; i < NN; i++) {
Inputpre[i]->set_I(0);
Inputpost[i]->set_I(0);
}
if (tim > 199 && tim < endtime) {
for (int i = 0; i < NN; i++) {
Inputpre[i]->set_I(curr_pre[i]);
Inputpost[i]->set_I(curr_post[i]);
}
}
/*if (tim > 199.25 && tim < 200.) {
for (int i = 0; i < NN; i++) {
Inputpre[i]->set_I(50.0);
Inputpost[i]->set_I(50.0);
}
}*/
if (noise != 0.0) {
for (int i = 0; i < NN; i++) {
do {
x1[i] = 2.0 * rg.Random() - 1.0;
x2[i] = 2.0 * rg.Random() - 1.0;
w[i] = x1[i] * x1[i] + x2[i] * x2[i];
}
while (w[i] >= 1.0);
do {
x1p[i] = 2.0 * rg.Random() - 1.0;
x2p[i] = 2.0 * rg.Random() - 1.0;
wp[i] = x1p[i] * x1p[i] + x2p[i] * x2p[i];
}
while (wp[i] >= 1.0);
neurpre[i]->AddNoise(noise, 0., 1., x1[i], x2[i], w[i],
timestep);
neurpost[i]->AddNoise(noise, 0., 1., x1p[i], x2p[i], wp[i],
timestep);
}
}
for (int i = 0; i < NN; i++) {
HH_Run(&tim, neurpre[i], timestep);
HH_Run(&tim, neurpost[i], timestep);
TwoD_Run(&tim,slowself[i][i],neurpre[i]->x[0], timestep);
TwoD_Run(&tim, slowselfA[i][i],neurpost[i]->x[0], timestep);
for (int j = 0; j < NN; j++) {
TwoD_Run(&tim,slowmutual[i][j],neurpre[i]->x[0], timestep);
TwoD_Run(&tim, slowmutualA[i][j],neurpost[i]->x[0], timestep);
}
}
if (tim > 200) {
for (int i = 0; i < NN; i++) {
spiketimes(&tim, &neurpre[i]->x[0], Threshold, bol[i], T[i]);
spiketimes(&tim, &neurpost[i]->x[0], Threshold, bolp[i],
TP[i]);
}
}
/*Getting the total current per cycle entering the neuron*/
if (szpreold!=T[0].size())
{
//cout<<"Entered loop "<<T[0].size()-1<<endl;
PreCurr.push_back(sumcurr);
sumcurr=0;
}
sumcurr+=I_ext+slowself[0][0]->Isyn()+slowmutualA[0][0]->Isyn();
szpreold=T[0].size();
if (szpostold!=TP[0].size())
{
//cout<<"Entered loop "<<TP[0].size()-1<<endl;
PostCurr.push_back(sumcurrpost);
sumcurrpost=0;
}
sumcurrpost+=Idc+slowselfA[0][0]->Isyn()+slowmutual[0][0]->Isyn();
szpostold=TP[0].size();
/* Learning Rule On Line */
if (Learnpresent == 1 ) {
dG = 0;
for (int u = 0; u < NN; u++)
for (int v = 0; v < NN; v++) {
maxsz = max(T[u].size(), TP[v].size());
minsz = min(T[u].size(), TP[v].size());
maxnew = maxsz;
minnew = minsz;
if (maxnew != maxold || minnew != minold) {
// if(TP[v].size()>1)
if (TP[v].size()>0 && T[u].size()>0)
dT = TP[v][TP[v].size() - 1] - T[u][T[u].size() - 1];
else
dT=-1000;
//RuleNoise=rg.Random();
// dG = Rule(dT, 10.0, 1.0, 10.0, 1.0, .01);
dG=Rule(dT,5,.25,5,.25,.015);
}
gmutualnewpre = gmutualoldpre + dG;
if (gmutualnewpre < 0)
gmutualnewpre = 0;
gmutualnewpost = gmutualoldpost - dG;
if (gmutualnewpost < 0)
gmutualnewpost = 0;
slowmutual[u][v]->parameter[0] = 1.0*gmutualnewpre;
slowmutualA[u][v]->parameter[0] = gmutualnewpost;
}
maxold = maxnew;
minold = minnew;
} else;
gmutualoldpre = 1.0*gmutualnewpre;
gmutualoldpost = gmutualnewpost;
count += 1;
if (scale==2)
{
if (tim>0)
Outfile<<tim<<" "<<neurpre[0]->x[0]<<" "<<neurpost[0]->x[0]<<" "<<gmutualnewpost<<" "<<gmutualnewpre<<endl;
}
tim += epsilon;
}
cout<<tim<<" "<<neurpre[0]->x[0]<<" "<<neurpost[0]->x[0]<<" "<<gmutualnewpre<<" "<<gmutualnewpost<<" "<<slowmutualA[0][0]->x[0]<<endl;
vector < double >Periodpre, Periodpost;
double Avgperiodpre, Avgperiodpost;
double sumpree = 0, sumposte = 0;
double **kura, sumkura = 0, avgkura;
kura = new double *[NN + 10];
for (int j = 0; j < NN + 10; j++)
kura[j] = new double[NN + 10];
for (int i = 0; i < NN; i++) {
Periodpre.push_back(ISI(T[i]));
Periodpost.push_back(ISI(TP[i]));
}
for (int i = 0; i < NN; i++) {
sumpree += Periodpre[i];
sumposte += Periodpost[i];
}
Avgperiodpre = sumpree / NN;
Avgperiodpost = sumposte / NN;
//cout<<I_ext<<" "<<Avgperiodpre<<" "<<Avgperiodpost<<endl;
/*
Outfile1.open((char *)SpikeTimeFile.c_str(), ios::out);
if (scale!=0)
{
for (int i=0;i<PreCurr.size();i++)
Outfile1<<PreCurr[i]<<endl;
Outfile1<<-1000<<endl;
for (int i=0;i<PostCurr.size();i++)
Outfile1<<PostCurr[i]<<endl;
for (int i=0;i<T[0].size();i++)
Outfile1<<T[0][i]<<endl;
Outfile1<<-1000<<endl;
for (int i=0;i<TP[0].size();i++)
{Outfile1<<TP[0][i]<<endl;}
}
Outfile1.close();
*/
for (int i = 0; i < NN; i++) {
int findpre = isnan(Periodpre[i]);
int findpreinf = isinf(1. / Periodpre[i]);
for (int j = 0; j < NN; j++) {
int findpost = isnan(Periodpost[j]);
if (findpre != 1 && findpost != 1 && findpreinf != 1
&& Periodpost[j] != 0.) {
kura[i][j] = Kuramoto(T[i], TP[j], Avgperiodpost);
sumkura += kura[i][j];
}
}
}
avgkura = sumkura / (NN * NN);
int MaxTime, MinTime;
MaxTime = max(T[0].size(), TP[0].size());
MinTime = min(T[0].size(), TP[0].size());
//cout << MaxTime << " " << MinTime << endl;
double diffspike = MaxTime - 2 * MinTime;
//Outfile<<I_ext<<" "<<Avgperiodpre<<endl;
if (diffspike < 100 && MinTime > 100) {
cout<< NoiseNum<<" "<<gmutual<<" "<<I_ext << " " << Het << " " <<gmutualnewpre << " "<<gmutualnewpost<<" " <<avgkura<<" "<< Avgperiodpre << " " <<Avgperiodpost << " " << Avgperiodpost / Avgperiodpre << endl;
}
double *rateA, *rateB;
rateA = new double[NN + 10];
rateB = new double[NN + 10];
Outfile.close();
Outfile1.close();
time(&end1);
double diff = difftime (end1, start);
cout<<"Computation time for "<< I_ext<<" is "<<diff<<" sec"<<endl;
delete *neurpre; delete *neurpost;delete *Inputpre;delete *Inputpost;delete **coupAB;
delete **coupBA; delete **slowself; delete **slowselfA; delete **slowmutual;
delete **slowmutualA;
delete[] neurpre;
delete[] neurpost;
delete[] Inputpre;
delete[] Inputpost;
delete[] coupAB;
delete[] coupBA;
delete[] slowself;
delete[] slowselfA;
delete[] slowmutual;
delete[] slowmutualA;
delete[] ref_time_pre;
delete[] ref_time_post;
delete[] vnew_pre;
delete[] vnew_post;
delete[] vold_pre;
delete[] vold_post;
delete[] diff_pre;
delete[] diff_post;
delete[] bol_pre;
delete[] bol_post;
delete[] bol_postpre;
delete[] bol_prepost;
delete[] T;
delete []TP;
delete []bol;
delete []bolp;
return 0;
}