ObjSim/src/normalize.cpp

/*Copyright (C) 2005, 2006, 2007 Frank Michler, Philipps-University Marburg, Germany

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 2
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, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include "sys.hpp"   // for libcwd
#include "debug.hpp" // for libcwd


#include "normalize.hpp"
#include "layer.hpp"

AbstractNormalize::AbstractNormalize()
     : SimElement(seNormalize), RewiringOn(false)
{
}

void AbstractNormalize::SetRewiring(float _IncommingConnectivity, float _SynDelThreshold, float _InitialWeights)
{
     RewiringOn=true;
     IncommingConnectivity=_IncommingConnectivity;
     SynDelThreshold=_SynDelThreshold;
     InitialWeights=_InitialWeights;
}

void AbstractNormalize::SetRewiringOff()
{
     RewiringOn=false;
}


////////Normalize//////////////////////////////////

Normalize::Normalize()
     : AbstractNormalize(), Target(0), NTarget(0)
{
}

int Normalize::AddConnection(connection* newcon)
{
     if (Target == 0) 
     {
	  Target = newcon->GetTargetLayer();
	  NTarget = Target->N;
     }
     if (Target == newcon->GetTargetLayer()) ConList.push_back(newcon);
     else cout << "ERROR: Target-Layer not the same\n";
}

int Normalize::proceede(int TotalTime)
{
}

int Normalize::prepare(int Step)
{
}

int Normalize::WriteSimInfo(fstream &fw)
{
     stringstream sstr;
     sstr << "<Target id=\"" << Target->IdNumber << "\"/> \n";
      SimElement::WriteSimInfo(fw, sstr.str());
}

int Normalize::WriteSimInfo(fstream &fw, const string &ChildInfo)
{
     stringstream sstr;
     sstr << "<Target id=\"" << Target->IdNumber << "\"/> \n";
     sstr << ChildInfo;
      SimElement::WriteSimInfo(fw, sstr.str());
}


/////////////////////

FiringRateNormalize::FiringRateNormalize(float _NormThresh, float _NormFactor, float _Tau)
     :  PostSynFirePot(0), PostSynLastFirings(0), Tau(_Tau/dt), NormThreshold(_NormThresh), NormFactor(_NormFactor)
{     
     cout << "FiringRateNormalization\n";
     cout << "Thresh=" << NormThreshold << "  Factor=" << NormFactor << "  Tau=" << Tau*dt << " ms\n";
}


int FiringRateNormalize::AddConnection(connection* newcon)
{
     Normalize::AddConnection(newcon);
     if (PostSynFirePot == 0) 
     {
	  PostSynFirePot = new float [NTarget];
	  PostSynLastFirings = new int [NTarget];
	  int i;
	  for (i=0;i<NTarget;++i) 
	  {
	       PostSynFirePot[i]=0;
	       PostSynLastFirings[i]=0;
	  }
     }
}

int FiringRateNormalize::proceede(int TotalTime)
{
     int t = int(TotalTime % MacroTimeStep);
     int spike = Target->last_N_firings;
     int CurTarget;
     int i,j;
     while (spike < Target->N_firings) {
	  CurTarget = Target->firings[spike][1];
	  PostSynFirePot[CurTarget] *= exp(-(t-PostSynLastFirings[CurTarget])/Tau);
	  if (PostSynFirePot[CurTarget] > NormThreshold) {
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i])  *= NormFactor;
	       }  
	  }
	  PostSynFirePot[CurTarget] += 1;
	  PostSynLastFirings[CurTarget]=t;
	  ++spike;
     }
}


int FiringRateNormalize::prepare(int Step)
{
     int i;
     for (i=0;i<NTarget;++i) PostSynLastFirings[i] -= MacroTimeStep;  
     //FixMe: what if neuron never fires?? prevent negative integer overflow??
}

//////////////////////////////

/*! \brief normalizing synapit weights if firing rates are above a threshold

 NormFrequency: above this spike frequency normalization occurs
 Weights are multiplied with (1-NormFactor)
 if spike frequency is higher then weight reduction is larger
 maximum: (1-NormFactor*MaxNormFactor)

 a look up table is used to determine the current normalization factor, 
depending on the current spike frequency (time difference DeltaT between current spike and last spike):

NormLut(DeltaT)=1-MaxNormFactor*NormFactor*exp(-DeltaT/Tau), 

weight is multiplied with NormLut(DeltaT)

@param  _NormFrequency threshold frequency, above this frequency normalization occurs
@param _NormFactor weight normalization factor, 
if postsynaptic neuron fires with _NormFrequency, synaptic weight is mulitiplied with 1-_NormFactor

@param _MaxNormFactor  maximal normalization: (1-NormFactor*MaxNormFactor)
If postsynaptic neuron fires with frequency higher than _NormFrequency, the normalization is stronger.
For infinite firing rate normalization strength can raise up to  _MaxNormFactor times.
@author (fm)
*/
FiringRateNormalize2::FiringRateNormalize2(
     float _NormFrequency, float _NormFactor, float _MaxNormFactor)
     : PostSynLastFirings(0), 
      MaxNormFactor(_MaxNormFactor),
      NormFactor(_NormFactor),
      NormFrequency(_NormFrequency)
{     
     
     NormDeltaT = 1000./(NormFrequency*dt);    
     Tau =  NormDeltaT/log(MaxNormFactor);

     // NormLut(DeltaT)=1-MaxNormFactor*NormFactor*exp(-DeltaT/Tau)
     NormLut = ExpDecayLut(NormLutN, Tau, -MaxNormFactor*NormFactor, 1, dt, NormDeltaT/Tau);

     cout << "FiringRateNormalization2\n";
     cout << "NormLut=" << "\n";
     for (int i=0;i<NormLutN;++i) cout << NormLut[i] << "\n";
     cout << "  Factor=" << NormFactor << " MaxNormFactor=" << MaxNormFactor << "  Tau=" << Tau*dt << " ms\n";
}

int FiringRateNormalize2::WriteSimInfo(fstream &fw)
{
     stringstream sstr;
     sstr << "<MaxNormFactor value=\"" << MaxNormFactor << "\"/>\n";
     sstr << "<NormFrequency value=\"" << NormFrequency << "\"/>\n";
     sstr << "<NormFactor value=\"" << NormFactor << "\"/>\n";     
     Normalize::WriteSimInfo(fw, sstr.str());
}

int FiringRateNormalize2::AddConnection(connection* newcon)
{
     Normalize::AddConnection(newcon);
     if (PostSynLastFirings == 0) 
     {
	  PostSynLastFirings = new int [NTarget];
	  int i;
	  for (i=0;i<NTarget;++i) 
	  {
	       PostSynLastFirings[i]=0;
	  }
     }
}

int FiringRateNormalize2::proceede(int TotalTime)
{
     int t = int(TotalTime % MacroTimeStep);
     int spike = Target->last_N_firings;
     int CurTarget;
     int i,j, TDiff;
     float NormFactor=1;
     while (spike < Target->N_firings) {
	  CurTarget = Target->firings[spike][1];
	  int TDiff = t-PostSynLastFirings[CurTarget];
	  if (TDiff < NormLutN) {
	       NormFactor = NormLut[TDiff];
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i])  *= NormFactor;
	       }  
	  }
	  PostSynLastFirings[CurTarget]=t;
	  ++spike;
     }
}


int FiringRateNormalize2::prepare(int Step)
{
     int i;
     for (i=0;i<NTarget;++i) PostSynLastFirings[i] -= MacroTimeStep;  
     //FixMe: what if neuron never fires?? prevent negative integer overflow??
}



//////////////////////////////

ConstSFNormalize::ConstSFNormalize(float _DesiredSpikeFreq, float _NormFactor, float _Tau)
     :  PostSynFirePot(0), PostSynLastFirings(0), Tau(_Tau/dt)
{     
     cout << "ConstSFNormalization\n";
     cout << "DesiredSpikeFrequency=" << NormThreshold << "  Factor=" << NormFactor << "  Tau=" << Tau*dt << " ms\n";

     // one spike potential value
     float DeltaT = 1000./(_DesiredSpikeFreq*dt); // inter spike interval in number of time steps
     float DesiredSglSpikePot = exp(-DeltaT/Tau);
     // spike train equilibrium
     DesiredFirePot = 1/(1-exp(-DeltaT/Tau));

     NormFactor = _NormFactor/(DesiredFirePot-1);
}


int ConstSFNormalize::AddConnection(connection* newcon)
{
     Normalize::AddConnection(newcon);
     if (PostSynFirePot == 0) 
     {
	  PostSynFirePot = new float [NTarget];
	  PostSynLastFirings = new int [NTarget];
	  int i;
	  for (i=0;i<NTarget;++i) 
	  {
	       PostSynFirePot[i]=0;
	       PostSynLastFirings[i]=0;
	  }
     }
}

int ConstSFNormalize::proceede(int TotalTime)
{
     int t = int(TotalTime % MacroTimeStep);
     int spike = Target->last_N_firings;
     int CurTarget;
     int i,j;
     while (spike < Target->N_firings) {
	  CurTarget = Target->firings[spike][1];
	  PostSynFirePot[CurTarget] *= exp(-(t-PostSynLastFirings[CurTarget])/Tau);
	  PostSynFirePot[CurTarget] += 1;
	  PostSynLastFirings[CurTarget]=t;
	  for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	  {
	       for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i])  += *((*it)->s_pre[CurTarget][i])*NormFactor*(DesiredFirePot-PostSynFirePot[CurTarget]);
	  }  
	  ++spike;
     }
}


int ConstSFNormalize::prepare(int Step)
{
     int i;
     for (i=0;i<NTarget;++i) PostSynLastFirings[i] -= MacroTimeStep;  
     //FixMe: what if neuron never fires?? prevent negative integer overflow??
}



int ConstSFNormalize::WriteSimInfo(fstream &fw)
{
     fw << "<" << seTypeString << " id=\"" << IdNumber <<  "\" Type=\"" << seType << "\" Name=\"" << Name << "\"> \n";
     fw << "<Target id=\"" << Target->IdNumber << "\"/> \n";
     fw << "<DesiredFirePot Value=\"" << DesiredFirePot << "\"/> \n";     
     fw << "<NormFactor Value=\"" << NormFactor << "\"/> \n";     
     fw << "</" << seTypeString << "> \n";
}


///////////////////////////////////

ConstSumNormalize::ConstSumNormalize(float _WeightSum, bool _quadratic)
     : WeightSum(_WeightSum), quadratic(_quadratic)
{     
     cout << "ConstSumNormalization\n";
     cout << "WeightSum=" << WeightSum << "  quadratic=" << quadratic  << " \n";
}


// int ConstSumNormalize::AddConnection(connection* newcon)
// {
//      Normalize::AddConnection(newcon);
// }

int ConstSumNormalize::proceede(int TotalTime)
{
     int t = int(TotalTime % MacroTimeStep);
     int spike = Target->last_N_firings;
     int CurTarget;
     int i,j;
     float CurWeightSum, NormFactor;
     float tmpweight;
     if (quadratic) 
     {
	  while (spike < Target->N_firings) {
	       CurTarget = Target->firings[spike][1];
	       // calculate WeightSum
	       CurWeightSum=0;
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) {
			 tmpweight= (*((*it)->s_pre[CurTarget][i]));
			 if (tmpweight <0) { // delete this thread
			      cout << "EEEEEEEEEERRRROOORRR, weight deletion didn't work\n";
			      fflush(stdout);
			      exit(2);
			 }
			 CurWeightSum += (*((*it)->s_pre[CurTarget][i]))*(*((*it)->s_pre[CurTarget][i]));
		    }
	       }  
	       // DEBUG
	       if (CurWeightSum > 100) {
		    for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
		    {
			 cout << "N_pre=" << (*it)->N_pre[CurTarget] << "\n";
			 for (i=0;i<(*it)->N_pre[CurTarget];++i) {
			      tmpweight= (*((*it)->s_pre[CurTarget][i]));
			      cout << "w"<< i << "=" << tmpweight << "I_pre=" <<  (*it)->I_pre[CurTarget][i] << "\n";
			 }
		    }  
		    fflush(stdout);
		    exit(2);
	       }
	       // END DEBUG
	       NormFactor = WeightSum/sqrt(CurWeightSum);
	       cout << "NormFactor=" << NormFactor << "WeightSum" << CurWeightSum << "\n";
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i]) *= NormFactor;
	       }  
	       ++spike;
	  }

     }
     else {
	  while (spike < Target->N_firings) {
	       CurTarget = Target->firings[spike][1];
	       // calculate WeightSum
	       CurWeightSum=0;
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) CurWeightSum += *((*it)->s_pre[CurTarget][i]);
	       }  
	       NormFactor = WeightSum/CurWeightSum;

	       // multiplicatively normalize weights
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i]) *= NormFactor;
	       }  
	       ++spike;
	  }
     }
}

int ConstSumNormalize::NormalizeAll()
{
     cout << "Normalizing All ...";
     int CurTarget;
     int i,j;
     float CurWeightSum, NormFactor;
     if (quadratic) 
     {
	  for (CurTarget=0;CurTarget<Target->N;++CurTarget) 
	  {
	       CurWeightSum=0; // wie oben
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) CurWeightSum += (*((*it)->s_pre[CurTarget][i]))*(*((*it)->s_pre[CurTarget][i]));
	       }  
	       NormFactor = WeightSum/sqrt(CurWeightSum);
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i]) *= NormFactor;
	       }  
	  }

     } 
     else {
	  for (CurTarget=0;CurTarget<Target->N;++CurTarget) 
	  {
	       // calculate WeightSum
	       CurWeightSum=0;
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) CurWeightSum += *((*it)->s_pre[CurTarget][i]);
	       }  
	       NormFactor = WeightSum/CurWeightSum;
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i]) *= NormFactor;
	       }  
	  }
     }
     cout << " [done]\n";
}


void ConstSumNormalize::CalcInitWeightSum()
{
     cout << "Calculating Initial Weight Sum ...";
     int CurTarget;
     int i,j;
     float CurWeightSum, NormFactor;
     float TmpWeightSum=0;
     if (quadratic) 
     {
	  for (CurTarget=0;CurTarget<Target->N;++CurTarget) 
	  {
	       CurWeightSum=0; // wie oben
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) CurWeightSum += (*((*it)->s_pre[CurTarget][i]))*(*((*it)->s_pre[CurTarget][i]));
	       }  
	       TmpWeightSum += sqrt(CurWeightSum);
	  }
	  WeightSum = TmpWeightSum/Target->N;

     } 
     else {
	  for (CurTarget=0;CurTarget<Target->N;++CurTarget) 
	  {
	       // calculate WeightSum
	       CurWeightSum=0;
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) CurWeightSum += *((*it)->s_pre[CurTarget][i]);
	       }  
	       TmpWeightSum += CurWeightSum;
	  }
	  WeightSum = TmpWeightSum/Target->N;
     }
     cout << " [done]\n";
}

float ConstSumNormalize::GetWeightSum()
{
     return WeightSum;
}

int ConstSumNormalize::SetWeightSum(float NewWeightSum)
{
     WeightSum=NewWeightSum;
}


int ConstSumNormalize::WriteSimInfo(fstream &fw)
{
     stringstream sstr;     
     sstr << "<WeightSum Value=\"" << WeightSum << "\"/> \n";    
     sstr << "<Quadratic Value=\"" << quadratic << "\"/> \n";
     Normalize::WriteSimInfo(fw, sstr.str());     
}


// int ConstSumNormalize::prepare(int Step)
// {
// }


/////////////////////

NormalizePsp::NormalizePsp(float _NormThresh, float _NormFactor)
     :  NormThreshold(_NormThresh), NormFactor(_NormFactor), PspPot(0)
{     
     cout << "NormalizePsp\n";
     cout << "Thresh=" << NormThreshold << "  Factor=" << NormFactor << "\n";
}


int NormalizePsp::AddConnection(connection* newcon)
{
     Normalize::AddConnection(newcon);
     if (PspPot == 0) 
     {
	  PspPot =  Target->GetPspPointer(csimInputChannel_AMPA); //ToDo: Info aus connectio nverwenden
     }

}


int NormalizePsp::proceede(int TotalTime)
{
     int t = int(TotalTime % MacroTimeStep);
     int spike = Target->last_N_firings;
     int CurTarget;
     int i,j;
     while (spike < Target->N_firings) {
	  CurTarget = Target->firings[spike][1];
	  if (PspPot[CurTarget] > NormThreshold) {
	       for(vector<connection*>::iterator it=ConList.begin(); it !=ConList.end(); ++it)
	       {
		    for (i=0;i<(*it)->N_pre[CurTarget];++i) *((*it)->s_pre[CurTarget][i])  *= NormFactor;
	       }  
	  }
	  ++spike;
     }
}

int NormalizePsp::WriteSimInfo(fstream &fw)
{
     stringstream sstr;     
     sstr << "<NormThreshold Value=\"" << NormThreshold << "\"/> \n";    
     sstr << "<NormFactor Value=\"" << NormFactor << "\"/> \n";    
     Normalize::WriteSimInfo(fw, sstr.str());     
}