#include <iostream>

#include <cmath>

#include <cstdlib>

#include <fstream>

#include <iomanip>

#include <ctime>

#include <string>

#include <sstream>

using namespace std;

#include "RandNum.h"

CRandNum randNum;

#define RAN randNum.GenRandReal_10()

#define RAN00 randNum.GenRandReal_00()

#define RANDINT randNum.GenRandInt32()

#define Norm randNum.Gaussian_Noise()

const double vth=-50;

const double vreset=-70;

const double taum=20;

const double ve=50;

const double vi=-75;

const double vr=-70;

const double gl=1;

const double PNconstexcitconmean = 0.24;

const double PNconstinhibconmean = 0.15;

const double PNconstexcitconvar = 0.02;

const double PNconstinhibconvar = 0.01;

const double inputscalePNconst=0.045;

const double inputscaleLNconst=0.013;

const double inputscaleinterPNconst=0.04;

const double inputscaleinterLNconst=0.004;

const double otherinhibconst=1;

const double trefra=2;

const int noofodourant=160;

const int noofodour=16;

const double adaptcurrentmaxPNmean=4.5;

const double adaptcurrentmaxLNmean=1.8;

const double adaptcurrentmaxPNvar=0.4;

const double adaptcurrentmaxLNvar=0.2;

const double tauadapt=25;

const int noofbisect=9;

const int noofiteration=100;

const double corrconnratio = 0.01;

const double noiseconnratio = 0.001;

const double baseconn = 0.006;

const double baseconnnoise = 0.002;

const double noisePN = 0.01;

const double noiseLN = 0.003;

const double noiseinterPN = 0.01;

const double noiseinterLN = 0.001;

const int noofconc=7;

const int minconc=-6;

const int nooffile = noofconc*2;

int main()

{

int i;

int j;

int m;

int k;

int n;

int fileindexcounter;

double condeff;

double temp;

double tempone;

double temptwo;

double tempthree;

double tempfour;

double tempfive;

double effinput;

double efftaum;

double rantempone;

double rantemptwo;

double adaptcurrentmaxPN[noofodourant];

double adaptcurrentmaxLN[noofodourant];

double ORNmaxresponse[noofodour][noofodourant];

double interconnectivity[noofodourant][noofodourant];

double rawinput[noofodourant];

double outputPN[noofodourant];

double outputLN[noofodourant];

double outputLNold[noofodourant];

double inputscalePN[noofodourant];

double inputscaleLN[noofodourant];

double inputscaleinterPN[noofodourant];

double inputscaleinterLN[noofodourant];

double PNconstexcitcon[noofodourant];

double PNconstinhibcon[noofodourant];

double otherinhibinput;

double v;

double t;

double tupper;

double tlower;

double adaptscale;

string filename;

stringstream filenameindex[nooffile];

srand(time(NULL));

for (i=0;i<noofodourant;i++)

{

do {rantempone = Norm;}

while (abs(rantempone)>2);

inputscalePN[i] = inputscalePNconst+noisePN*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

inputscaleLN[i] = inputscaleLNconst+noiseLN*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

inputscaleinterPN[i] = inputscaleinterPNconst+noiseinterPN*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

inputscaleinterLN[i] = inputscaleinterLNconst+noiseinterLN*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

PNconstexcitcon[i] = PNconstexcitconmean+PNconstexcitconvar*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

PNconstinhibcon[i] = PNconstinhibconmean+PNconstinhibconvar*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

adaptcurrentmaxPN[i] = adaptcurrentmaxPNmean+adaptcurrentmaxPNvar*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

adaptcurrentmaxLN[i] = adaptcurrentmaxLNmean+adaptcurrentmaxLNvar*rantempone;

}

//find connectivity between LN and PN+other LN

for (i=0;i<noofodourant;i++)

interconnectivity[i][i]=0;

ifstream inputone;

inputone.open("ORNfr0.txt");

for (k=0;k<noofodour;k++)

for (j=0;j<noofodourant;j++)

inputone>>ORNmaxresponse[k][j];

inputone.close();

//calculate correlation matrix for ORN-PN connectivity

ofstream outtest;

outtest.open("corrmatrix1.txt");

for (j=1;j<noofodourant;j++)

{

for (i=0;i<j;i++)

{

tempone=0;

temptwo=0;

tempthree=0;

tempfour=0;

tempfive=0;

for (k=0;k<noofodour;k++)

{

tempone=tempone+ORNmaxresponse[k][i];

temptwo=temptwo+ORNmaxresponse[k][j];

tempthree=tempthree+ORNmaxresponse[k][i]*ORNmaxresponse[k][i];

tempfour=tempfour+ORNmaxresponse[k][j]*ORNmaxresponse[k][j];

tempfive=tempfive+ORNmaxresponse[k][i]*ORNmaxresponse[k][j];

}

tempone=tempone/(double)noofodour;

temptwo=temptwo/(double)noofodour;

tempthree=tempthree/(double)noofodour;

tempfour=tempfour/(double)noofodour;

tempfive=tempfive/(double)noofodour;

temp=(tempfive-tempone*temptwo)/(sqrt((tempthree-tempone*tempone)*(tempfour-temptwo*temptwo)));

outtest<<temp<<" ";

if (temp<0)

{

do {rantempone = Norm;}

while (abs(rantempone)>2);

interconnectivity[i][j]=baseconn+baseconnnoise*rantempone;

do {rantempone = Norm;}

while (abs(rantempone)>2);

interconnectivity[j][i]=baseconn+baseconnnoise*rantempone;

}

else

{

do {rantempone = Norm;}

while (abs(rantempone)>2);

do {rantemptwo = Norm;}

while (abs(rantemptwo)>2);

interconnectivity[i][j]=baseconn+baseconnnoise*rantempone+temp*(corrconnratio+noiseconnratio*rantemptwo);

do {rantempone = Norm;}

while (abs(rantempone)>2);

do {rantemptwo = Norm;}

while (abs(rantemptwo)>2);

interconnectivity[j][i]=baseconn+baseconnnoise*rantempone+temp*(corrconnratio+noiseconnratio*rantemptwo);

}

}

outtest<<endl;

}

outtest.close();

fileindexcounter=0;

for (m=0;m<noofconc;m++)

{

ifstream infile;

filenameindex[fileindexcounter]<<"ORNfr"<<minconc+m<<".txt";

filename=filenameindex[fileindexcounter].str();

fileindexcounter=fileindexcounter+1;

infile.open(filename.c_str());

ofstream outfile;

filenameindex[fileindexcounter]<<"PNfr"<<minconc+m<<".txt";

filename=filenameindex[fileindexcounter].str();

fileindexcounter=fileindexcounter+1;

outfile.open(filename.c_str());

for (j=0;j<noofodour;j++)

{

for (i=0;i<noofodourant;i++)

infile>>rawinput[i];

for (i=0;i<noofodourant;i++)

outputLNold[i]=0;

//iterate until steady state

for (n=0;n<noofiteration;n++)

{

for (i=0;i<noofodourant;i++)

{

//find total (raw) inhibitory input received

otherinhibinput=0;

for (k=0;k<noofodourant;k++)

otherinhibinput=otherinhibinput+interconnectivity[i][k]*outputLNold[k];

otherinhibinput=otherinhibinput*otherinhibconst;

//find response of PN

condeff=1/(1+(rawinput[i]*inputscalePN[i]+PNconstexcitcon[i]+PNconstinhibcon[i]+inputscaleinterPN[i]*otherinhibinput));

effinput=(ve*(rawinput[i]*inputscalePN[i]+PNconstexcitcon[i])+vi*(inputscaleinterPN[i]*otherinhibinput+PNconstinhibcon[i])+vr*gl)*condeff;

efftaum=taum*condeff;

v=vreset;

t=0;

if (effinput<=vth)

outputPN[i]=0;

else

{

adaptscale = tauadapt*adaptcurrentmaxPN[i]*(sqrt(outputPN[i]))/(tauadapt-efftaum);

while (v<vth)

{

t=t+1;

v=vreset*exp(-t/efftaum)+effinput*(1-exp(-t/efftaum))-adaptscale*(exp(-t/tauadapt)-exp(-t/efftaum));

}

tupper=t;

tlower=t-1;

t=(tupper+tlower)*0.5;

for(k=0;k<noofbisect;k++)

{

v=vreset*exp(-t/efftaum)+effinput*(1-exp(-t/efftaum))-adaptscale*(exp(-t/tauadapt)-exp(-t/efftaum));

if (v>vth)

tupper=t;

else

tlower=t;

t=(tupper+tlower)*0.5;

}

outputPN[i]=1000/(t+trefra);

}

//find response of LN

condeff=1/(1+(rawinput[i]*inputscaleLN[i]+inputscaleinterLN[i]*otherinhibinput));

effinput=(ve*rawinput[i]*inputscaleLN[i]+vi*inputscaleinterLN[i]*otherinhibinput+vr*gl)*condeff;

efftaum=taum*condeff;

v=vreset;

t=0;

if (effinput<=vth)

outputLN[i]=0;

else

{

adaptscale = tauadapt*adaptcurrentmaxLN[i]*(sqrt(outputLN[i]))/(tauadapt-efftaum);

while (v<vth)

{

t=t+1;

v=vreset*exp(-t/efftaum)+effinput*(1-exp(-t/efftaum))-adaptscale*(exp(-t/tauadapt)-exp(-t/efftaum));

}

tupper=t;

tlower=t-1;

t=(tupper+tlower)*0.5;

for(k=0;k<noofbisect;k++)

{

v=vreset*exp(-t/efftaum)+effinput*(1-exp(-t/efftaum))-adaptscale*(exp(-t/tauadapt)-exp(-t/efftaum));

if (v>vth)

tupper=t;

else

tlower=t;

t=(tupper+tlower)*0.5;

}

outputLN[i]=1000/(t+trefra);

}

}

for (i=0;i<noofodourant;i++)

outputLNold[i]=outputLN[i];

}

for (i=0;i<noofodourant;i++)

outfile<<outputPN[i]<<" ";

outfile<<endl;

}

infile.close();

outfile.close();

}

}