/* 
    Version: $Id: constructs.cpp 172 2014-02-12 10:06:07Z gk $
    LAModel main imlementation file

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "constructs.h"
//#include "wxglmodel.h"
#include <iterator>
#include <assert.h>

const int DEBUG_NID = 91;
const int DEBUG_BID = 2270;
const int DEBUG_SID = 3024;
const int CUTOFF = 10.0;

int LANetwork::RSEED = 1980;


static void updminmax(float&  cmin, float& cmax, float val)
{
	if (val < cmin )
		cmin = val;
	if (val > cmax )
		cmax = val;
}





#define VEC_REMOVE(vec, t) (vec).erase(std::remove((vec).begin(), (vec).end(), t), (vec).end())

inline float sigmoid(float x, float x0, float broad)
{
	return (1.0 / (1.0 + exp(-x - x0)/broad));
}




inline double nuclearproteinalpha(float x)
{
	return (x>20.)*((x-20.*60.)/(30.*60)) * exp(1. - (x-20.*60. )/(30.*60.));

	//double d =  (x-60.*60)/(40.0*60.);
	//return (exp(-d*d));

	//double d =  ((x-3.*60)/(23.*60)) * exp(1. - (x-(3.*60.) )/(23.*60.)); // 4. or 6.
	//if (d >0.) return d;
	//else return 0.;
	
}




inline double branchproteinalpha(float x)
{
	return ((x)/(15.*60)) * exp(1. - (x )/(15.*60.));

	
}


inline void adjust_with_bounds(float& val, float direction, float max, float min)
{
	if (direction>=0.0) val += direction*(max-val);
	else val += direction*(val - min);
}


inline float step_with_bounds(float cur_val, float dir, float max, float min)
{
	if (dir >0) return dir*(max - cur_val);
	else  return dir*(cur_val - min);
}




inline float caDP(float x)
{
	//return x >0.2 ? 1.0 : 0.0;
	//float f =  (2.0/(1.+exp(-(x*10.-3.5)*10.))) - (1.0/(1.+exp(-(x*10.0-0.5)*19.)));


	float f = (1.3/(1.+exp(-(x*10.-3.5)*10.))) - (0.3/(1.+exp(-(x*10.0-2.0)*19.)));
	return f;
	//return f;//(2.0/(1.+exp(-(x*10.-0.7)*10.))) - (1.0/(1.+exp(-(x*10.0-0.2)*19.)));
	//return (2.0/(1.+exp(-(x*10.-3.1)*8.))) - (1.0/(1.+exp(-(x*10.0-2.1)*6.)));
}


inline void program_input(nrn_list lst, int tstart, int duration, float freq, float randomness, int limitActive = -1)
{
	int skip = 0;
	if (limitActive == -999)
		skip = 3;
	for (nrn_iter n = lst.begin(); n != lst.end(); ++n)
	{
		if (skip>0)
		{
			skip--;
			continue;
		}

		LAInput* in = (LAInput*)(*n);
		in->Program(tstart, duration, freq, randomness);
		//if (limitActive >0 && ++tot >= limitActive) return;
	}
}


void LANetwork::CreateNeurons(int number, int n_branches_per_neuron, char type, vector<LANeuron*>* appendTo = 0, int inputId =-1, int nBranchesTurnover = 0)
{
	for (int i =0 ;i < number; i++)
	{
		LANeuron* n ;
		if (type == 'S')
		{
			n = new LAInput();
			n->network = this;

			n->input_id = inputId;
			((LAInput*)n)->groupIdx = i;
		}
		else
		{
			n = new LANeuron();
			n->network = this;
			for (int tt =0; tt < n_branches_per_neuron; tt++)
			{
				LABranch* bb  = new LABranch;
				bb->bid = this->branches.size();
				bb->neuron = n;
				if (tt < nBranchesTurnover)
				{
					if (this->turnoverHotspots ==2)
						bb->turnoverRate = 0.5 + rgen()*0.5;
					else
						bb->turnoverRate = 1.0;

				}
				else
					bb->turnoverRate = 0.0;

				this->branches.push_back(bb);
				n->branches.push_back(bb);
			}
		}
		n->type = type;

		n->V = param_V_rest +1.0;
		n->nid = this->neurons.size();

		n->glx = 1.9*(n->nid%50)/50. - 0.9;
		n->gly = 1.9*(n->nid/50)/50. - 0.9;

		this->neurons.push_back(n);
		if (appendTo)
			appendTo->push_back(n);
	}
}



void LANetwork::CalculateDistances()
{
	for (nrn_iter na = neurons.begin(); na != neurons.end(); ++na)
		for (nrn_iter nb = neurons.begin(); nb != neurons.end(); ++nb)
			if (nb != na)
			{
				LANeuron* a = *na, *b=*nb;
				
				float dx = b->pos_x - a->pos_x;
				float dy = b->pos_y - a->pos_y;
				float dz = b->pos_z - a->pos_z;
				distances[pair<int,int>(a->nid,b->nid)] = (double)sqrt(dx*dx + dy*dy + dz*dz);
			}
}


/*
int LANetwork::ConnectToRandomBranches(vector<LANeuron*> fromList, vector<LANeuron*> toList,   int nSynapses, float weight)
{
		LANeuron* from = fromList.at(int(rgen()*(float)toList.size()));
		LANeuron* to = toList.at(int(rgen()*(float)toList.size()));

		LASynapse* syn = new LASynapse();
		syn->sid  = this->synapsesCounter++;
		syn->target_branch = b;

		syn->source_nrn = from;
		syn->target_nrn = to;
		syn->isPlastic = isPlastic;
		syn->weight  = weight;

		syn->target_branch = syn->target_nrn->branches[(int)rval];

		syn->source_nrn->outgoing.push_back(syn);
		syn->target_nrn->incoming.push_back(syn);
		syn->target_branch->synapses.push_back(syn);
		syn->pos = rgen();
		this->synapses.push_back(syn);

}

*/


void LANetwork::AddSynapse(LANeuron* a, LABranch* br, float weight, bool isPlastic)
{
		LASynapse* syn = new LASynapse();
		syn->sid  = this->synapsesCounter++;
		syn->source_nrn = a;
		syn->target_nrn = br->neuron;
		syn->isPlastic = isPlastic;

		syn->weight  = weight;
		syn->target_branch = br; 

		syn->source_nrn->outgoing.push_back(syn);
		syn->target_nrn->incoming.push_back(syn);
		syn->target_branch->synapses.push_back(syn);
		syn->pos = rgen();
		this->synapses.push_back(syn);
}



int LANetwork::ConnectByDistance(vector<LANeuron*> fromList, vector<LANeuron*> toList, float fromDistance, float toDistance, int nNeuronPairs, int nSynapsesPerNeuron, float weight, bool isPlastic= false, bool randomizeweight = false, float overlap =-1.0)
{
	int tpairs =0;
	while(true)
	{
		LANeuron* a = fromList.at(int(rgen()*(float)fromList.size()));
		LANeuron* b = toList.at(int(rgen()*(float)toList.size()));

		for (int i =0; i < nSynapsesPerNeuron; i++)
		{

			float rval;
			rval = rgen()*float(b->branches.size());
			LABranch* br =  b->branches[(int)rval];

			this->AddSynapse(a, br, weight, isPlastic);	
		}

		if (++tpairs >= nNeuronPairs) break;
	}
	return tpairs;
}




inline int randomindex(int max)
{
	return (int)(rgen()*float(max));
}

int LANetwork::PurgeInputSynapses(int totalToRemove,float weightLimit)
{

	int totalFound =0;
	int totalTries = totalToRemove*3;
	while (totalFound < totalToRemove && totalTries-->0)
	{
		nrn_list lst = this->inputs_cs.at(randomindex(this->inputs_cs.size()));
		LANeuron* n = lst.at(randomindex(lst.size()));
		if (n->outgoing.size())
		{
			LASynapse* s = n->outgoing.at(randomindex(n->outgoing.size()));
			if (s->target_nrn->type == 'P'  && s->weight <= weightLimit)
			{
				//candidate for deletion

				//pthread_mutex_lock(&this->synapses_mutex);

				VEC_REMOVE(s->source_nrn->outgoing, s);
				VEC_REMOVE(s->target_nrn->incoming, s);
				VEC_REMOVE(s->target_branch->synapses, s);
				VEC_REMOVE(this->synapses, s);
				//cout << " Removing " << s->sid << endl;
				delete s; 

				//pthread_mutex_unlock(&this->synapses_mutex);

				totalFound++;
			}
		}
	}

	if (totalFound < totalToRemove)
	{
		cout << " Warning: not enougn synapses to remove: " << totalToRemove << endl;  
	}

	return totalFound;
}



int LANetwork::CreateInputSynapses( int totalToAdd)
{
	
	int totalFound = 0;
	while (totalFound < totalToAdd)
	{
		this->ConnectByDistance(this->inputs_cs[randomindex(this->inputs_cs.size())], this->pyr_list, 0.0, 10.0,  1, 1, initWeight, true);
		totalFound++;
	}
	return totalFound;
}



void LANetwork::CreateFearNet(int nneurons, int nbranches, int ninputs, int nneuronsperinput)
{
	this->n_neurons = nneurons;
	this->n_inputs = ninputs;
	this->n_neurons_per_input = nneuronsperinput;
	this->n_branches_per_neuron = nbranches;

	this->CreateNeurons(this->n_neurons*0.8, this->n_branches_per_neuron, 'P', &this->pyr_list, -1, this->nBranchesTurnover);

	this->CreateNeurons(this->n_neurons*0.2, 1 , 'I', &this->in_list);
	
	for (int i=0;i < n_inputs;i++)
	{
		vector<LANeuron*> nrnsCS;
		CreateNeurons(this->n_neurons_per_input, 0, 'S', &nrnsCS, i, true);
		this->inputs_cs.push_back(nrnsCS);
	}


	CreateNeurons(10, 0, 'S', &this->bg_list, -1, true);

	this->CalculateDistances();

	if (1) 
	{
	/*
		int nconn = 1.0*float(this->pyr_list.size()*this->in_list.size());
		nrn_list pl, nl;
		for (int i=0; i < nconn; i++)
		{
			pl.clear();
			nl.clear();
			pl.append(this->pyr_list.at(randomindex(this->pyr_list.size())));
			nl.append(this->in_list.at(randomindex(this->pyr_list.size())));
		}
		*/
		ConnectByDistance(this->pyr_list, this->in_list, 0., 10., this->inhibitionParam*2.*.4*20.*float(this->pyr_list.size()), 1, 2.0, false);

		ConnectByDistance(this->in_list, this->pyr_list, 0.0, 10.,this->inhibitionParam*2.*.6*20.*float(this->pyr_list.size()), 1, 2.0, false);

	}

	// CS inputs  projections
	for (int i =0; i < this->n_inputs ; i++)
	{
		
		this->ConnectByDistance(this->inputs_cs[i], this->pyr_list, 0.0, 10.0, this->connectivityParam*float(this->pyr_list.size()*2), 1, initWeight, true, false, this->branchOverlap);

		/* 
		Quote:
		-Synaptic inputs are not clustered in fast-spiking parvalbumin-positive interneurons
		-Inhibitory interneurons exhibited highly localized alcium activity in their aspiny dendrites, 
			as reported previously (S10, 11). Because they are highly excitable and can fire in response to a 
			single excitatory input (S12), they may not require dendritic integration via synaptic clustering. 
		For another type of interneuron, see ref. S10.
		Takahashi et al. Science 2012 (supplement figure S6) 
		*/

	}


	ConnectByDistance(this->bg_list, this->pyr_list, 0.0, 10., 20*float(this->pyr_list.size()), 1, 0.6, false);

	//this->spikesFile = fopen("./data/0/spikes.dat", "w");
	
	/*&
	for (int i=0; i < 0.20*this->pyr_list.size(); i++)
	{
		//LANeuron* a = this->pyr_list.at(int(rgen()*(float)this->pyr_list.size()));
		LANeuron* a = this->pyr_list.at(i);
		if (!this->disableCreb)
		{
			//a->crebLevel = 1.0;
		}
		//a->prpTransients.push_back( pair<float,float>(0, 1.) );
	}	
	*/

	this->RecordInitialWeightSums();
}



void LANetwork::RunPattern(vector<int>& pattern, float hifreq,  float lowfreq, int duration, int limitActive )
{

	int pad = (4000 - duration)/2;
	//printf("Running pattern: [ ");


	for (size_t j=0; j < pattern.size(); j++)
	{
		if (pattern[j]>0)
		{
			
				/*
			double ts =0;
			for (syn_iter si=  this->synapses.begin() ; 
				 si != this->synapses.end() ;
				si++)
			{
				LASynapse* s = *si;
				if (s->source_nrn->input_id == (int)j)
				{
					ts+= s->weight;
				}
			}
			//printf("Total synaptic drive: %f\n", ts);
			//*/


			program_input(this->inputs_cs[j], pad, duration, hifreq, 0.5, limitActive);
		}
		else
			program_input(this->inputs_cs[j], pad, duration, lowfreq, 0.5, limitActive);

		//printf(" %d", pattern[j] ? 1 : 0);
	}

	program_input(this->bg_list, pad, duration, .5, 0.5, -1);

	//printf("] @%f Hz dur=%d \n", hifreq, duration);

	this->ResetSpikeCounters();
	this->StimDynamics(duration+pad+pad);
	
	int tActive =0;
	int tSpikes =0;
	
	for (nrn_iter ni = this->pyr_list.begin(); ni != this->pyr_list.end(); ni++)
	{
		LANeuron* n = *ni;
		if (float(n->total_spikes) /4000.> CUTOFF )
			tActive++;
		tSpikes += n->total_spikes;
	}
	
	printf("Active: %d (%.2f%%), avgF: %.2f\n", tActive, 100.0*float(tActive)/float(this->pyr_list.size()), tSpikes/(float(this->pyr_list.size())*2));
}



void LANetwork::RunStoreTest(int n_patterns, int activePerPattern, int delayMinutes, int testMode, int patternsOverlapping)
{
	vector<int> pattern;

	int inpfrequency=60;
	int lowfreq = 0;
	int multiruns = 1;


	int n_ones = activePerPattern;
	printf("Running net with %d pyr. neurons, %d branches, %d synapses [%s,%s] [%d per pattern]\n", (int)this->pyr_list.size(),  (int)branches.size(),  (int)synapses.size(), localProteins ? "Local" : "Global", disableCreb ? "-CREB" : "+CREB", n_ones);
	//this->mfile.open("./data/bdata.dat");
	//this->vfile.open("./data/bvoltage.dat");

	//fill(pattern.begin(), pattern.end(), 0);
	for (int i =0; i < this->n_inputs; i++) pattern.push_back( (i < n_ones) ? 1 : 0);

	//std::random_shuffle(pattern.begin(), pattern.end());

	int cp=0;
	int np=0;
	for (int i =0; i < n_patterns; i++)
	{
		//fill( pattern.begin(), pattern.end(), 0);

		if (this->repeatedLearning)
		{
			this->patterns.push_back( pattern );
		}
		else
		{
			fill( pattern.begin(), pattern.end(), 0);

			//cp=0;
			for (np = 0; np < n_ones; np++)
				if (cp < n_inputs)
					pattern[cp++] = 1;

			//std::random_shuffle(pattern.begin(), pattern.end());
			
			this->patterns.push_back( pattern );
		}

		cout<< "Pattern  " << i<< " : [";
		copy(pattern.begin(), pattern.end(), ostream_iterator<int>(cout, ""));
		cout << "]" << endl;
	}

	np =0;

	this->Interstim(1*60);

	PrintSynapsesSnapshot( datadir + "/syn-pre.txt");

	this->ReportSumWeights();

	if (this->pretraining)
	{
		this->runningMode = RUN_PRE;
		this->enablePlasticity = false;
		for (int nr=0; nr < multiruns; nr++)
		{
			np =0;
			for (vector<vector<int> >::iterator it = this->patterns.begin(); it != this->patterns.end(); it++)
			{
				pattern = *it;
				cout<< "Pretraining" << np<< endl;
				this->runningPatternNo = np;

				RunPattern(pattern, inpfrequency, lowfreq, stimDurationParam*3800., n_ones/2);

				cout << "Tiny interstim  ..." << endl;
				this->Interstim(5*60);
				cout << "done" << endl;
				np++;
			}
		}
		
	}

	cout << "Training .. " << endl;
	this->enablePlasticity = true;
	np =0;
	this->runningMode = RUN_TRAINING;
	for (vector<vector<int> >::iterator it = this->patterns.begin(); it != this->patterns.end(); it++)
	{
		pattern = *it;
		cout<< "Training" << np<< endl;
		this->runningPatternNo = np;

		if (std::find(isWeakMem.begin(), isWeakMem.end(), np) != isWeakMem.end())
		{
			printf("Weak: %d\n", np);
			RunPattern(pattern, inpfrequency, lowfreq, 2700, -1);
		}
		else
			RunPattern(pattern, inpfrequency, lowfreq,  3800, -1);

		cout << "Interstim " << delayMinutes << " mins ..." << endl;
		//this->SaveCalcs();
		this->Interstim(delayMinutes*60);
		cout << "done" << endl;
		np++;
		this->ReportSumWeights();
	}


	this->runningMode = RUN_TEST;
	cout << "Large interstim ..." << endl;
	this->enablePlasticity = false;
	for (nrn_iter na = neurons.begin(); na != neurons.end(); ++na)
	{
		LANeuron* nrn = *na;
		nrn->crebLevel = 0.0;
	}
	this->Interstim((int)(this->homeostasisTime*3600.));
	cout << " done" << endl;


	this->ReportSumWeights();

	this->isRecall = true;
	cout << "Recall .. " << endl;
	np =0;
	this->enablePlasticity = false;
	int n =0;
	


	PrintSynapsesSnapshot(datadir + "/syn-post.txt");

	this->spikesPerPattern.resize(n_patterns);
	for (int i=0; i < n_patterns; i++)
		this->spikesPerPattern[i].resize(this->neurons.size());


	this->enablePlasticity = false;

	for ( vector<vector<int> >::iterator it = patterns.begin(); it != patterns.end(); it++ )
	{
		pattern = *it;
		n++ ;
		for (nrn_iter na = neurons.begin(); na != neurons.end(); ++na)
		{
			LANeuron* nrn = *na;
			nrn->crebLevel = 0.0;
		}

		cout<< endl<< "Testing " << np<< endl;
		this->runningPatternNo = np;
		RunPattern(pattern, inpfrequency, lowfreq, 3800, n_ones/2);
		this->Interstim(60);

		for (nrn_iter na = neurons.begin(); na != neurons.end(); ++na)
		{
			LANeuron* nrn = *na;
			this->spikesPerPattern[np][nrn->nid] = nrn->total_spikes;
		}
		np++;
	}

	/*

	np =0;
	n=0;


	for ( vector<vector<int> >::iterator it = patterns.begin(); it != patterns.end(); it++ )
	{
		pattern = *it;
		n++ ;

		for (nrn_iter na = neurons.begin(); na != neurons.end(); ++na)
		{
			LANeuron* nrn = *na;
			nrn->crebLevel = 0.0;
		}


		cout<< endl<< "Testing US " << np<< endl;
		this->runningPatternNo = np;
		RunPattern(pattern, inpfrequency, lowfreq, 3800, -999);

		this->Interstim(60);
		for (nrn_iter na = neurons.begin(); na != neurons.end(); ++na)
		{
			LANeuron* nrn = *na;
			this->spikesPerPattern[np][nrn->nid] = nrn->total_spikes;
		}
		np++;
	}

	*/


	//this->mfile.close();

}





template <typename T> static void PrintVector( vector<T>&  ar, ostream& outfile) 
{
	for (typename vector<T>::iterator it = ar.begin(); it != ar.end(); it++)
	{
		outfile << *it << ' ';
	}
	outfile << std::endl;
}





void LANetwork::StoreDataFiles( bool extras = false )
{
	vector<int> pattern;

	string dirname = this->datadir;
	ofstream paramsdat((dirname + "/parameters.txt").c_str());
	paramsdat <<"total_neurons="<< this->neurons.size() << endl;
	paramsdat <<"total_pyramidals="<< this->pyr_list.size() << endl;
	paramsdat <<"branches_per_neuron="<< this->n_branches_per_neuron << endl;
	paramsdat <<"number_inputs="<< this->inputs_cs.size() << endl;
	paramsdat <<"neurons_per_input="<< this->n_neurons_per_input << endl;
	paramsdat <<"rseed="<< RSEED << endl;


	ofstream patternsdat((dirname + "/patterns.txt").c_str());
	for (vector<vector<int> >::iterator it = this->patterns.begin(); it != this->patterns.end(); it++)
	{
		pattern = *it;
		copy(pattern.begin(), pattern.end(), ostream_iterator<int>(patternsdat, " "));
		patternsdat << endl;
	}

	ofstream synstatedat((dirname + "/synstate.dat").c_str());

	ofstream spikesdat((dirname + "/spikes.dat").c_str());
	ofstream crebdat((dirname + "/creb.dat").c_str());
	ofstream voltagedat((dirname + "/voltages.dat").c_str());
	ofstream branchspikesdat((dirname +"/branchspikes.dat").c_str());
	ofstream branchcalcium((dirname + "/branchcalcium.dat").c_str());
	ofstream weightsdat((dirname + "/weights.dat").c_str());
	ofstream branchproteins((dirname + "/branchproteins.dat").c_str());
	ofstream branchstrengths((dirname + "/branchstrengths.dat").c_str());
	ofstream tagsdat((dirname + "/tags.dat").c_str());
	ofstream nrnproteindat((dirname + "/nrnprotein.dat").c_str());
	ofstream weighthistorydat((dirname + "/weighthistory.dat").c_str());
	ofstream dbgneuron((dirname + "/dbgneuron.dat").c_str());

	PrintVector<float>( dbgNeuron, dbgneuron);

	for (nrn_iter na = neurons.begin(); na != neurons.end(); ++na)
	{
		LANeuron* nrn = *na;

		PrintVector<int>( nrn->spikeTimings, spikesdat);
		//PrintVector<float>( nrn->crebHistory, crebdat);
		PrintVector<float>( nrn->proteinHistory, nrnproteindat);

		for (branch_iter bi = nrn->branches.begin(); bi != nrn->branches.end(); ++bi)
		{
			LABranch* b = *bi;

			PrintVector<float>(b->branchSpikesHistory, branchspikesdat);

			for (syn_iter si = b->synapses.begin(); si != b->synapses.end(); ++si)
			{
				LASynapse* s = *si;
				//PrintVector<float>(s->weightHistory, weightsdat);
				synstatedat << s->sid<<" " 
					<< b->bid<<" "
					<< nrn->nid  << " "
					<< s->source_nrn->nid <<" " 
					<< s->source_nrn->input_id<< " "
					<< b->strength  << " "
					<< s->weight << " " 
					<<endl;

				if (s->tagHistory.size())
				{
					tagsdat << s->source_nrn->input_id<< " ";
					PrintVector<float>(s->tagHistory, tagsdat);
				}

				if (s->weightHistory.size())
				{
					weighthistorydat << s->source_nrn->input_id<< " ";
					PrintVector<float>(s->weightHistory, weighthistorydat);
				}
			}
		}
	}

	ofstream sppdat((dirname + "/spikesperpattern.dat").c_str());
	for (uint i =0; i < this->spikesPerPattern.size(); i++)
	{
		for (uint j =0; j < this->spikesPerPattern[i].size(); j++) sppdat << this->spikesPerPattern[i][j] << " ";
		sppdat << endl;
	}
}



void LANetwork::StimDynamics(int duration) // stimulation dynamics, with dt=1msec 
{
	int t = 0;
	bool spikeState[this->neurons.size()+1];
	int lastSpikeT[this->neurons.size()+1];

	fill_n(spikeState, this->neurons.size()+1, 0);
	fill_n(lastSpikeT, this->neurons.size()+1, 0);

	//FILE* outfile 	= fopen("./data/0/wdata.dat", "w");

	for (syn_iter si = this->synapses.begin(); si != this->synapses.end(); ++si)
	{
		(*si)->calcium = 0.0;
	}


	for (nrn_iter ni = this->neurons.begin(); ni != this->neurons.end(); ++ni)
	{
		LANeuron* n = *ni;
		n->wadapt = 0.0;
		n->vspike =0.;
		n->vreset =0.;
		n->V =0.;
	}



	for (branch_iter bi=this->branches.begin(); bi != this->branches.end(); ++bi)
	{
		(*bi)->totcalc = 0.0;
		(*bi)->depol = 0.0;
		(*bi)->dspike = 0.0;
		(*bi)->dspikeT = -1;
	}

	for (t=0; t < duration; t++)
	{
		for (nrn_iter ni = this->neurons.begin(); ni != this->neurons.end(); ++ni)
		{
			LANeuron* n = *ni;
			double s_inh = 0.0;
			double sv =0.0 ; //, ss =0;
			
			for (branch_iter bi=n->branches.begin(); bi != n->branches.end(); ++bi)
			{
				LABranch* b = *bi;
				float sumepsp =0.;
				float sumipsp =0.;

				//if (this->debugMode &&  b->bid==DEBUG_BID) vfile<< b->depol << " " << n->V << endl;

				for (syn_iter si=b->synapses.begin(); si != b->synapses.end(); ++si)
				{
					LASynapse* s = *si;
					if (spikeState[s->source_nrn->nid])
					{

						if (s->source_nrn->type == 'I') sumipsp += ( s->weight);
						else sumepsp += (  s->weight);
					}
				}


				b->depol +=  (4.0*sumepsp) -  b->depol/20.0;

				s_inh += sumipsp;

				//if (b->bid == 40) printf("BID %d f=%f epsp=%f\n", b->bid, b->depol, sumepsp);
				// http://www.cell.com/neuron/abstract/S0896-6273(12)00761-1: Inhibition blocks only weak dendritic spikes

				/*
				if ( b->dspikeT > 0)
				{
					float tds = t - b->dspikeT;
					if (tds < 40)
						b->depol = 40.;
					else if (tds < 60)
					{
					}
					else
					{
						b->dspikeT = -1; // Allow next dend spike
						b->depol =0.;// Reset dspike
					}
				}
				*/

			

				if (n->type == 'P' && b->dspikeT < t-100 && (n->vspike + b->depol) > 30.*dendSpikeThresh) // Generate a dendritic branch spike
				{
					//b->depol = 60.0;
					b->depol = 50;
					b->dspikeT =t;
					b->branch_spikes++;
					n->dend_spikes++;
					if (this->enablePlasticity && b->strengthTag < 1.0) b->strengthTag +=  (1.0 - b->strengthTag)*0.20;
				}

			//	if (b->dspikeT >0 && t-b->dspikeT < 20) b->depol = 50.;



				for (syn_iter si=b->synapses.begin(); si != b->synapses.end(); ++si)
				{
					LASynapse* s = *si;
					if (spikeState[s->source_nrn->nid])
					{
						if (this->enablePlasticity && s->isPlastic) // && t > 50b)
						{
							float depol =  b->depol + n->vspike;
							//float depol =  n->vspike;
							//NMDA dynamics: Jahr_1993
							//if (depol > 10.0)
							if (depol > 1.0) //n->vspike > 10)
							{
								//float ff = exp(-dt/70.0 ) * (1.0/(1.+exp( (-(depol-30.0)/7.0))));
								float ff =  (1.0/(1.+exp( (-(depol-30.0)/5.0))));

								//s->calcium +=  ( 0.5*ff)*(1.0 - s->calcium) ;
								//s->calcium +=  2.0*(ff);
								//cout << " Calc " << s->calcium<< endl;
								s->calcium +=  ff/10.;
							}
						}
					}
				}




				sv +=  (b->depol)*b->strength ;

				//if (b->dspike > 0.0) b->dspike -= (b->dspike*b->dspike)/40.0; // see Losoczy 2008  for epsp durations
				//b->branchVoltageHistory.push_back(b->depol1 + n->vspike + b->dspike);
				//b->branchCalciumHistory.push_back(bcalc);
			}


			n->inh_cur += 0.18*s_inh; 

			if (n->type == 'S') // Source neuron
			{
				LAInput* in = (LAInput*)n;
				//if (t ==0) printf ("nid=%d , T=%d, next=%d\n", in->nid, t, in->nextSpikeT);
				if (in->spikeTimes && t >= in->nextSpikeT && in->nextSpikeT >0)
				{
					//printf("ID %d spike! t=%d\n", in->nid, t );
					if (in->curSpike < in->totalSpikes)
						in->nextSpikeT = in->spikeTimes[in->curSpike++];
					else 
						in->nextSpikeT = -1;
					spikeState[in->nid] = 1;
					n->isSpiking = true;
					lastSpikeT[in->nid] = t;
					in->total_spikes++;
					in->V = 20.0; // threshold

				}
				else
				{
					spikeState[in->nid] = 0;
					n->isSpiking = false;
					in->V = param_E_L;
				}
			}
			else /// PYR OR IN NEURON
			{

				if (spikeState[n->nid])
				{
					//spikeState[n->nid] =0;
					n->V = 0.0;
					n->wadapt += 0.18;
				}
				
				n->exc_cur = sv*0.12;

				if (n->type == 'I') // Interneuron
	 			{
					n->V +=  (n->exc_cur - n->inh_cur) - (n->V)/20. - n->wadapt*(n->V+10.0);
					n->wadapt -= n->wadapt/70.;
				}
				else
				{
					n->V +=  n->exc_cur - 3.0*s_inh - (n->V)/30. -   n->wadapt*(n->V+10.0) ;

					if (this->disableCreb)
						n->wadapt -= n->wadapt/180.;
					else
						n->wadapt -= n->wadapt/((180. - 70.0*(n->crebLevel>0.2 ? 1. : 0.)));
				}


				if ( lastSpikeT[n->nid] < t-2 && n->V > (20.0 - (n->crebLevel>100.2 ? 2.0 : 0.) ))
				{
					spikeState[n->nid] = 1;
					lastSpikeT[n->nid] = t;
					n->total_spikes++;
					n->isSpiking = true;
					n->vspike = 30.0; 
					n->V = 70.;

				}
				else
				{
					spikeState[n->nid] = 0;
					n->isSpiking = false;
				}

				if (n->vspike > 0) n->vspike -= n->vspike / 17.0; 
				// Legenstein&maas eta_reset = 20msec
				

				if (n->wadapt <-0.) n->wadapt =0.;
				//if (n->wadapt < 0.) n->wadapt =0.;

			}

			if (spikeState[n->nid])
			{
				n->spikeTimings.push_back(t+ this->Tstimulation);
			}
		}

		#ifdef WXGLMODEL_H
		if (this->wx_window && t%10==0)
		{
			this->wx_window->UpdatePanel();
		}
		#endif
		

		//usleep(10000);
	}



	for (branch_iter bi=this->branches.begin(); bi != this->branches.end(); ++bi)
	{
		LABranch* b = *bi;
		b->branchSpikesHistory.push_back(b->branch_spikes);
	}


	//fflush(stdout);
	//fflush(spikesFile);

	this->Tstimulation += duration;

}


void LANetwork::SaveSnapshot(char* filename)
{
	ofstream of(filename);
	for (nrn_iter ni = this->pyr_list.begin(); ni != this->pyr_list.end(); ni++)
	{
		LANeuron* nrn = *ni;
		if (nrn->type == 'P')
		{
			float totTag =0;
			float totProtein =0;
			for (branch_iter bi = nrn->branches.begin(); bi != nrn->branches.end(); ++bi)
			{
				LABranch* b = *bi;
				for (syn_iter si = b->synapses.begin(); si != b->synapses.end(); ++si)
				{
					LASynapse* s = *si;
					if (s->isPlastic)
					{
						totTag += s->stag;
						//of<<s->sid << " "<< s->source_nrn->input_id << " "<< b->bid << ' ' << nrn->nid <<' '<< b->protein  <<' '<< s->stag <<' '<< endl;
					}
				}
				totProtein += b->protein;
			}
			of << nrn->nid << ' ' << totProtein << ' '<<  totTag << endl;
		}
	}

}



void LANetwork::Interstim(int durationSecs)
{
	int tstop = T + durationSecs;
	this->isInterstim = true;

	printf("Interstim %d seconds (T=%d) plast=%d G=%d, L=%d ... \n", durationSecs, T, this->enablePlasticity, this->globalProteins, this->localProteins);

	float tstep = 60.0;
	int totalWeak =0;
	float weightLimit =  initWeight + 0.0;

	// Count the total weak afferent synapses
	for (input_iter ii = this->inputs_cs.begin(); ii != this->inputs_cs.end(); ++ii)
	{
		nrn_list lst = *ii;
		for (nrn_iter n = lst.begin(); n != lst.end(); ++n)
			for (syn_iter si=(*n)->outgoing.begin(); si != (*n)->outgoing.end(); ++si)
				if ((*si)->weight <= weightLimit)
					totalWeak++;
	}

	int trec =0, tactrec=0;
	int totTagged=0, totTaggedMinus=0, totBProd =0;
	float totLTP=0., totLTD=0.;
	int totact =0, totSact =0;
	int totbspikes =0, totSpikes=0;
	float maxSpk =0;
	
		float actmin=9999, actmax=0;
		float nactmin=9999, nactmax=0;

	for (nrn_iter ni = this->pyr_list.begin(); ni != this->pyr_list.end(); ni++)
	{
		LANeuron*  nrn = *ni;
		float nrnCalc =0.0;
		if (nrn->total_spikes > 4)
			totSact++;


		for (branch_iter bi = nrn->branches.begin(); bi != nrn->branches.end(); ++bi)
		{
			LABranch* b = *bi;
			totbspikes += b->branch_spikes;

			if (!this->enablePlasticity)
				continue;

			for (syn_iter si = b->synapses.begin(); si != b->synapses.end(); ++si)
			{
				LASynapse* s =*si;
				float ctag = caDP(s->calcium);

				if (fabs(s->stag) < 0.1) /// Avoid erasing existing tags
				{
					s->stag = ctag;
					
					if (s->stag > 0.1)
					{
						totTagged++;
						totLTP += s->stag;
						//if (s->stag > 0.5 && b->prpTransients.size()>0) cout <<"Candidate is "<<s->sid<< " " << b->bid << " "<< nrn->nid<<endl;
					}
					if (s->stag < -0.1)
					{
						totTaggedMinus++;
						totLTD += s->stag;
					}	
				}

				/* */
				//printf("SID=%d Calc=%f tag=%f %f\n", s->sid, s->calcium, ctag, s->stag);
				//s->stag = ctag ;
				//* (1.0 - s->stag) ; // * ((1.- 1. / (1.+ exp(-(totw-2.5)*1.5))));
				//else if (ctag < -0.05) s->stag =  ctag  ; // * (s->stag - (-1.0)) ;  // * ((1.- 1. / (1.+ exp(-(totw-2.5)*1.5))));
				//
				//if (s->stag > 1.5 && nrn->prpTransients.size()==1) printf("SID=%d tag=%f PR=%f\n", s->sid, s->stag, nrn->proteinRate);

				//if (nrn->nid  == 65) printf("ST=%f %d, pr=%f\n", s->stag, s->sid, b->proteinRate);
				//if (s->stag > 0.2) printf("ID=%d %f\n", s->sid, s->stag);



				//s->eltp = caDP(s->calcium);
				//if (s->calcium >0.0) printf("NID=%d sid=%d calc=%f tag=%f\n", nrn->nid, s->sid, s->calcium, s->stag);
				b->totcalc += s->calcium;
				s->calcium = 0.;
			}


			
			if ( b->totcalc  > this->localPRPThresh*BPROD_CUTOFF) // This branch should produce PRPs now BPROD
			{
				b->prpTransients.push_back( pair<float,float>(T, b->totcalc));
				//printf ("BID %d  NID %d BCALC=%f NCALC=%f S=%d\n", b->bid, nrn->nid, b->totcalc, nrnCalc, (int)b->prpTransients.size());
				totBProd++;
			}

			nrnCalc +=  b->totcalc;

		}


		if (nrn->total_spikes > CUTOFF*4.0)
		{
			totSpikes += nrn->total_spikes;
			totact++;
			//printf("ACT %d\t %.2f\t%d\t%.2f\n",  nrn->nid  , nrn->crebLevel  , (nrn->total_spikes) , nrnCalc);

			updminmax(actmin, actmax, nrnCalc);

		}
		else
			updminmax(nactmin, nactmax, nrnCalc);



		if (maxSpk < nrn->total_spikes)
			maxSpk = nrn->total_spikes;


		if (this->enablePlasticity)
		{
			if (nrnCalc > this->globalPRPThresh*GPROD_CUTOFF) /// GPROD
			{
				nrn->prpTransients.push_back( pair<float,float>(T, nrnCalc) );

				if (!this->disableCreb ) 
					nrn->crebLevel=1.0;

				if (nrn->total_spikes > CUTOFF*4)
					tactrec ++;
				trec ++;
				//printf("id %d nrnCalc=%f spk=%d, dspk=%d CREB=%f trans=%d\n", nrn->nid, nrnCalc, nrn->total_spikes, nrn->dend_spikes, nrn->crebLevel,(int)nrn->prpTransients.size());
			}
		}



		nrn->totcalc  = nrnCalc;

	}



	printf("\n\nAct=[%f,%f] nonact=[%f,%f] \n", actmin, actmax, nactmin, nactmax);

	
	if (this->runningMode == RUN_TRAINING)
	{
		char buf[256];
		sprintf(buf,  "./data/r%d.dat", this->runningPatternNo);
		//this->SaveSnapshot(buf);
	}
	
	printf("TAGGED: %d/%d +%.1f/%.1f GPROD:%d (%.1f%%) BPROD:%d, ACT:%d (%.1f%% ), act+gprod:%d FREQ:%.1f max %.1f DSPK:%d sact:%d\n", totTagged, totTaggedMinus, totLTP, totLTD, trec, 100.*(float)trec/(float(this->pyr_list.size())), totBProd, totact, 100.*(float)totact/(float(this->pyr_list.size())), tactrec, (float)totSpikes/((float)this->pyr_list.size()*4.0), float(maxSpk)/4.0, totbspikes, totSact);

	for (; T < tstop; T+= tstep)
	{
		for (nrn_iter ni = this->pyr_list.begin(); ni != this->pyr_list.end(); ni++)
		{
			LANeuron*  nrn = *ni;

			float totalSynapseWeight =0.0;
			float totalBranchStrength =0.0;
			int totalSynapses =0;
			nrn->protein =0.0;

			// "Whole-neuron" distribution of proteins
			nrn->proteinRate =0;
			for (pair_iter ii = nrn->prpTransients.begin(); ii != nrn->prpTransients.end(); ++ii)
			{
				pair<float, float> p = *ii;
				int td= (T - p.first);
				float al = (nuclearproteinalpha(td));
				if (nrn->proteinRate < al)
					nrn->proteinRate =  al;
			}

			
			
			for (branch_iter bi = nrn->branches.begin(); bi != nrn->branches.end(); ++bi)
			{
				LABranch* b = *bi;
				b->proteinRate =0.;

				for (pair_iter ii = b->prpTransients.begin();ii != b->prpTransients.end(); ++ii)
				{
					pair<float, float> p = *ii;
					float td = float(T - p.first);
					float al = (branchproteinalpha(td));
					if (b->proteinRate < al)
						b->proteinRate = al;
				}


				float  f =0.;

				if (this->localProteins)
					f = 1.0*b->proteinRate; 
				else if (this->globalProteins)
					f =  1.0* nrn->proteinRate;
				else
				{
					f = 1.0*b->proteinRate + 1.0* nrn->proteinRate;
					if (f>1.0) f = 1.0;
				}



				b->protein = f; // += tstep*(f/2000. - b->protein/3600.);


				vector<LASynapse*> candidates;


				for (syn_iter si = b->synapses.begin(); si != b->synapses.end(); ++si)
				{
					LASynapse* s =*si;

					if (s->stag != 0.0) 
					{
						s->stag -= (tstep/3600.)* s->stag;
						//s->stag -= (s->stag>0. ? +1. : -1.)*tstep / (3.4*3600);


						//if (s->sid == 3024) printf("STAG=%f\n", s->stag);

						if (b->protein > 0.1 && (s->stag >0.1 || s->stag < 0.1))
						{
							//candidates.push_back(s);
							float fw = s->stag* b->protein;
							//if (s->stag >0.)
							s->weight += tstep * fw/400.;
						}
					}

					if(1)
					{
						if (s->weight  > maxWeight)
						{
							s->weight = maxWeight;
							//s->stag =0.;
						}
						else if (s->weight  < 0.)
						{
							s->weight = 0.;
							//s->stag =0.;
						}
					}


					totalSynapseWeight += s->weight;
					totalSynapses++;

					//Homeostasis
			
					//if (1) //this->runningMode  != RUN_TRAINING)
					//if (this->runningMode  == RUN_TEST)
					if (1)
					{
						s->weight += s->weight * (1.0 - nrn->synScaling )*tstep/(7.*24.0* 3600.*homeostasisTimeParam); // Synaptic scaling (Ref: http://www.sciencedirect.com/science/article/pii/S0896627308002134)
					}
					//if (s->weight > 1.2) printf("SID =  %d %f\n", s->sid, s->weight);


					if (this->debugMode && b->bid == DEBUG_BID)
					{
						//cout << b->proteinRate << " " << nrn->proteinRate << " " << b->protein << " " << s->stag << " " << s->weight << " " << endl;
						//this->mfile << " "  <<   s->stag << " " << s->weight ;
					}
				}

				/*
				if (candidates.size()>0 && b->protein >0.1)
				{
					LASynapse* s = candidates.at(int(rgen()*(float)candidates.size()));

					// Consolidation
					float q = b->protein* (s->stag>0 ? 1.0 : -1.0) *tstep/600.0;
					float dw;

					if (q>=0)
						dw =  q; //(MAXSTRENGTH - s->weight);
					else
						dw = q; // *(s->weight - 0.0);
					s->weight += dw;
					s->stag -= q*1.0;

					totCons += dw;

					if (s->weight  > MAXSTRENGTH)
					{
						s->weight = MAXSTRENGTH;
						//s->stag =0.;
					}
					else if (s->weight  < 0.)
					{
						s->weight = 0.;
						//s->stag =0.;
					}
				}
				*/



				//if (this->debugMode && b->bid == DEBUG_BID) this->mfile << endl;



				if (0)  // Branch strength potentiation 
				{
					if (b->strengthTag>0.01)
					{
						b->strength += tstep*b->strengthTag*(1.5- b->strength) / (2.0*60.*60.*BSPTimeParam); // Max branch strength=2, time to get there: 40 minutes losonczy 2008 figure 3 http://www.nature.com/nature/journal/v452/n7186/abs/nature06725.html
						b->strengthTag -= tstep*b->strengthTag / (10.*60.);
					}
					totalBranchStrength += b->strength;
					b->strength += tstep*b->strength*(1.0 - nrn->branch_scaling) / (3.0*3600.); // Branch homeostatic scaling 
				}


				if (T%800 ==0)
				{
					b->branchProteinHistory.push_back(b->protein);
					//b->branchStrengthHistory .push_back(b->strength);
				}
			}

			//if (T%500==0) nrn->proteinHistory.push_back(nrn->proteinRate);

			// Synaptic homeostasis / synaptic scaling
			//if (nrn->synapticWeightsInitialSum != 0.0 )
			if (totalSynapses>0)
				nrn->synScaling = totalSynapseWeight / (initWeight*float(totalSynapses));
		
			else
				nrn->synScaling = 1.0;

			// Branch plasticity homeostasis
			nrn->branch_scaling = totalBranchStrength/((float(1.0) * float(nrn->branches.size())));


			//creb Drop
			if (nrn->crebLevel >0.0)
			{
				nrn->crebLevel -= tstep/(3600.*8.*CREBTimeParam );
				//if (nrn->nid == DEBUG_NID) printf("CREBlev=%f\n", nrn->crebLevel);
				//nrn->crebLevel -= nrn->crebLevel * tstep / (3600.0*4.);
			}


		
		}


		#ifdef WXGLMODEL_H
		if (this->wx_window && T%20 ==0)
		{
			this->wx_window->UpdatePanel();
		}
		#endif
	}
	


	//printf("CONSOLIDATED: %f\n", totCons);
	


	// zuo et al. 2005
	// We assume 100% of  filopodia turn over over 24 hours. 
	// We assume that weak synapses (W < 0.3) are filopodia

	/*
	if (this->enablePlasticity)
	{
		
		this->isPruningSynapses = true;
		int synapsesToPrune = 0.005*float(totalWeak) * float(durationSecs)/3600.0;
		cout << "removing "<< synapsesToPrune << "synapses" << endl;
		int found = this->PurgeInputSynapses(synapsesToPrune, weightLimit);
		if (found > 0)
			this->CreateInputSynapses(found);
		cout << "Renewed " << found <<" synapses"  << endl;
		this->isPruningSynapses = false;

	}
	*/

	if (this->enablePlasticity)
		this->DoTurnover(durationSecs);

	this->isInterstim = false;
}



void LANetwork::DoTurnover(float durationSecs )
{

	float pTurnOver = 4.*durationSecs/86400.; // per day

	int synsAdded=0;
	for (branch_iter bi = this->branches.begin(); bi != this->branches.end(); bi++)
	{
		LABranch* b = *bi;
		if (b->turnoverRate>0.0)
		{
			vector<LASynapse*> v;
			for (syn_iter si = b->synapses.begin(); si != b->synapses.end(); si++)
			{
				LASynapse* s = *si;
				if (s->isPlastic && s->weight <= this->initWeight)
				{
					float p;
					if (this->turnoverHotspots==0 || this->turnoverHotspots == 2)
						p = b->turnoverRate*0.1981; // uniform
					else 
					{
						float dist = s->pos - b->turnoverPoint;
						p = exp(-(dist*dist)/.1)*b->turnoverRate;
					}

					if (rgen() < p*pTurnOver)
					{
						VEC_REMOVE(s->source_nrn->outgoing, s);
						VEC_REMOVE(s->target_nrn->incoming, s);
						VEC_REMOVE(this->synapses, s);
						//cout << " Removing syn " << s->sid << endl;

						si = b->synapses.erase(si); //VEC_REMOVE(s->target_branch->synapses, s);

						delete s; 

						/* Connect  a random input back here */

						nrn_list lst = this->inputs_cs.at(randomindex(this->inputs_cs.size()));
						LANeuron* n = lst.at(randomindex(lst.size()));
						this->AddSynapse(n, b, initWeight, true);
						synsAdded++;

						continue;
					}
					
				}
			}

		}
	}

	printf("Added/removed %d synapses\n", synsAdded);
		

}



void LANetwork::Begin()
{
	
}