ObjSim/src/libcsim.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.
*/

// c++ simulation library 
// layers and connection matrices as objects
// goal: structured and expandable but still fast (avoid too much function calls)
// Created by Frank Michler, September 2005, Marburg, Germany
//
// chunks of code were taken from SPNET (Eugene Izhikevich) 
// http://vesicle.nsi.edu/users/izhikevich/publications/spnet.cpp


// 5.10.2005: Idee: ein Objekt (oder structure) mit globalen Einstellungen
//  (dt, MacroTimeStep, ...), und jedes SimElement hat eine Referenz
// oder einen Pointer darauf [done]
// 13.10.2005: gsl_rng * r;  /* global generator */ sollte auch darin enthalten sein
// --> k�nnte alles in SimLoop integriert werden [done]


// 02.11.2005
// ToDo: check existence of directory before saving files
// use assertions [done]

// 03.11.2005
// ToDo: MySimLoop: save sim info with all layer and connection names
// an idl program could then automatically show all layers [done]


// 15.11.2005
// ToDo: ask and create directory if it doesn't exist [done]

// 16.11.2005
// ToDo: bug: bei Delay-Differenz == 0 tritt Gleitkomma-Fehler auf (wahrscheinlich Division durch Null) --> abfangen!! (wahrscheinlich modulo 0 bei Zufalls-Generierung der Delays: getrandom(DelayDiff))
// 18.11.2005 [done]


// 17.11.2005
// ToDo: Idee: Zeiger auf MainSimLoop global machen?? 
// so dass er nicht mehr mit jedem Konstruktor �bergeben werden muss [done]

// 29.11.2005
// Fehlermeldungen nach cerr ausgeben!!

// 09.12.2005
// ToDo: input-Klassen: TestMode abschaffen und durch enum mode ersetzen

// 21.12.2005
// ToDo: Gewichts-Matritzen: im Moment ns*M-Matritzen, Problem: wenn wenige Source-Neuronen auf alle Target-Neuronen konvergieren: M==nt und die Matrix wird sehr gro�
// L�sung: Die Gewichts-Arrays nur so gro� machen wie n�tig, und zus�tzlich die Gr��e speichern


// 22.12.2005
// ToDo: Handling der Delays: MaxDelay muss Connection-spezifisch gesetzt werden (und evtl. minDelay/maxDelay ber�cksichtigen)
// 08.01.2006
// ToDo: Ber�cksichtigung des individuellen maximumDelay [done]
// ToDo: Ber�cksichtigung des individuellen minimumDelay

// 22.02.2006
// ToDo: Format der Gewichtsstrukturen:
// delays[ns][maximumDelay][M] --> zu gro�, muss nicht bis M in der dritten Dim gehen
// Problem: Format muss an vielen Stellen ber�cksichtigt werden (Einlesen, Speichern, IDL, Tcl/Tk)

// 24.02.2006
// ToDo: AnyOptionWrapper::SaveDefaultConfigFile
// aus den gesetzten Parametern ein Default-Config-File generieren [done]

// 19.06.2006
// ToDo: automatically save configuration file to data directory
#include "sys.hpp"
#include "debug.hpp"

#include "objsimlibrary.hpp"

// #include "anyoptionwrapper.cpp"
// #include "simmodules.cpp"

// float gauss(float x, float sigma=1, float x0=0)
// {
//      return exp(-(x-x0)*(x-x0)/(2*sigma*sigma));
// }

double gauss(double x, double sigma, double x0)
{
     return exp(-(x-x0)*(x-x0)/(2*sigma*sigma));
}




/////////////////
LookUpTable::LookUpTable(): value(0), tmpvalue(0), N(0), xRange0(0), xRange1(0), RDiff(0), IndexFactor(0)
{
}

LookUpTable::~LookUpTable()
{
     if (tmpvalue != 0) delete [] tmpvalue;
}


void LookUpTable::Reset()
{
     if (tmpvalue != 0) delete [] tmpvalue;
     tmpvalue=0;
     value=0;
     N=0;
     xRange0=0;
     xRange1=0;
     RDiff=0;
}

void LookUpTable::InitExp(int _N, float Tau, float v0, float Range0, float Range1)
{
     Reset();
     xRange0=Range0;
     xRange1=Range1;
     RDiff = xRange1-xRange0;
     N=_N;
     IndexFactor = N/RDiff;
     LBound = int(xRange0*IndexFactor);
     UBound = int(xRange1*IndexFactor);
     float coef=RDiff/N;

     if (N > 0) {
	  if (RDiff >0) {
	       tmpvalue = new float [N+1];
	       value = tmpvalue-LBound;
	       for (int t=LBound;t<=UBound;++t) {
		    value[t] = v0*exp(-coef*t/Tau);
	       }
	  } else {
	       tmpvalue = new float [N];
	       value = tmpvalue;
	       for (int t=0;t<N;++t) value[t] = v0*exp(-t/Tau);
	  }
     }
}

float LookUpTable::GetValue(int index)
{
     if (index > 0) {
	  if (index <N) {
	       return value[index];
	  } else return value[N-1];
     } else return value[0];
}

double LookUpTable::GetValue(double X)
{
     int index = int(IndexFactor*X);
//     cout << "(X-xRange0)/RDiff=" << N*(X-xRange0)/RDiff << " index=" << index << "\n";
     if (index > LBound) {
	  if (index <UBound) {
	       return value[index];
	  } else return value[UBound];
     } else return value[LBound];
}


/////////////////////
DistanceProfile::DistanceProfile(int _N): N(_N)
{
}


DistanceProfile::DistanceProfile(const vector<float> &Profile)
{
     cout << "DistanceProfile::DistanceProfile, inizialize with vector" << "\n";
     N=Profile.size();
     profile = new float [N];   
     for (int i=0;i<N;++i)
     {
	  profile[i]=Profile[i];		       
     }
}


DistanceProfile::DistanceProfile(int _N, float Sigma1, float Sigma2): N(_N)
{
     cout << "DistanceProfile::DistanceProfile, N=" << N << "\n";
     float MaxX=1;
     profile = new float [N];   
     int i;
     double x;
     double s1, s2;
     for (i=0;i<N;++i)
     {
	  x = i*MaxX/N;
	  s2= gauss(x, Sigma2);
	  s1=gauss(x, Sigma1);
	  profile[i] = float(s2-s1);
	  // if (i<20) cout << "x=" << x << "  s2=" << s2 << " s1=" << s1 << " profile[i]=" << profile[i] << "\n";
     }
     // find maximum
     float gmax =0;
     for (i=0;i<N;++i)
     {
	  if (profile[i] > gmax) gmax=profile[i];
     }
     float factor=1.0/gmax;
     cout << "gmax=" << gmax << " factor=" << factor << "\n";
     // scale to [0,1]
     for (i=0;i<N;++i)
     {
	  profile[i] *= factor;
     }          
}

DistanceProfile::~DistanceProfile()
{
     cout << "DistanceProfile::destructor\n";fflush(stdout);
     delete profile;     
}

float DistanceProfile::GetValue(float dist)
{
     int i = int(dist*N);
     if (i>=N) i=N-1;
     return profile[i];
}

float DistanceProfile::GetMaxConDistance(float minweight)
{
     int i=N-1;
     while((i>0) && (profile[i--]<minweight));
     return float(i)/N;
}

void DistanceProfile::print()
{
     cout << "DistanceProfile::print()  N=" << N << "\n";
     int i;
     for(i=0;i<N;++i) cout << profile[i] << "  ";
     cout << "\n";
}

CircleDistanceProfile::CircleDistanceProfile(int _N, float Radius, bool invers): DistanceProfile(_N)
{
     cout << "CircleDistanceProfile::CircleDistanceProfile, N=" << N << "\n";
     float MaxX=1;
     profile = new float [N];   
     int i;
     double x;     
     for (i=0;i<N;++i)
     {
	  x = i*MaxX/N;
	  if (invers) {
	       if (x>=Radius) profile[i]=1; else profile[i]=0;		       
	  } else {
	       if (x<=Radius) profile[i]=1; else profile[i]=0;
	  }
     }
}
///////////////////////////////////////////////

Recorder::Recorder(char* _FileName): FileName(_FileName)
{
     fs = fopen(FileName.c_str(),"w");  
}

Recorder::~Recorder()
{
     cout << "Recorder Destructor, FileName=" << FileName << "\n";
     fclose(fs);
}

int Recorder::record(float time, float value)
{
     fprintf(fs, "%f  %f\n", time, value);         
}

int Recorder::record(float time, float value, float value2)
{
     fprintf(fs, "%f  %f %f\n", time, value, value2);         
}
////////////////////////////////////////////////


BinRecorder::BinRecorder(int _MacroTimeStep, int _N, int _M, float** _PBuffer, float _dt, const char* _FileName): MacroTimeStep(_MacroTimeStep), N(_N), M(_M), PBuffer(_PBuffer), BaseFileName(_FileName), SimTag("")

{
     int i;
     const char* tmpinfo = "BinRec\n";
     const char* tmpversion = "0.02\n";
     strcpy(Header.finfo, tmpinfo);
     strcpy(Header.version, tmpversion);
     Header.N=N;
     Header.M=M;
     Header.ValueSize=sizeof(**PBuffer);
     Header.FrameSize = Header.N*Header.ValueSize;
     Header.dt = _dt;
     FileName = BaseFileName + SimTag;
     fs = fopen(FileName.c_str(),"w");  
     // Write Header Infos
     fwrite(&Header, sizeof(Header), 1, fs);
     fclose(fs);
     cout << "BinRecorder; MacroTimeStep =" << MacroTimeStep << "  N=" << N << "\n";
     RecBuffer = new float [N*MacroTimeStep];
     RecBuffPointer = RecBuffer;
     RecBuffMax = RecBuffPointer + N*(MacroTimeStep-1)+1;  //??
}

BinRecorder::~BinRecorder()
{
//     fclose(fs);
     delete [] RecBuffer;
}

int BinRecorder::record()
{
     int i;
//      fs = fopen(FileName.c_str(),"a");  
//      for (i=0; i<N;++i) fwrite(PBuffer[i], sizeof(**PBuffer), 1, fs);
//      fclose(fs);
//     cout << "a"; fflush(stdout);
     float test;
     if (RecBuffPointer < RecBuffMax)
	  for (i=0; i<N;++i) *(RecBuffPointer++) = *PBuffer[i];
//	  for (i=0; i<N;++i) test = *PBuffer[i];
     else cout << "ERROR: BinRecorder: RecBuffer overflow\n";
//     cout << "g"; fflush(stdout);
}

int BinRecorder::restart()
{
     fs = fopen(FileName.c_str(),"w");  
     fwrite(&Header, sizeof(Header), 1, fs);
     fclose(fs);
}

int BinRecorder::restart(string _SimTag)
{
     FileName = BaseFileName + _SimTag;
     fs = fopen(FileName.c_str(),"w");  
     fwrite(&Header, sizeof(Header), 1, fs);
     fclose(fs);
}

int BinRecorder::WriteToFile()
{
     fs = fopen(FileName.c_str(),"a");  
     fwrite(RecBuffer, N*MacroTimeStep*sizeof(*RecBuffer), 1, fs);
     fclose(fs);
     RecBuffPointer = RecBuffer;
}

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

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

SimpleTextProgressBar::SimpleTextProgressBar(int MaxSteps): Final(MaxSteps), Status(0), StepSize(0.1)
{
}

int SimpleTextProgressBar::Next(int StepNr)
{
	while ((Status*Final) <= StepNr) 
     {
	  cout << "."; fflush(stdout);
	  Status += StepSize;
     }
}

int SimpleTextProgressBar::Reset(int MaxSteps)
{
     Final=MaxSteps;
     Status=0;
}


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

CSimStandardOptions::CSimStandardOptions(AnyOptionWrapper* _aowrap)
: myAnyWrap(_aowrap), NTrials(1), NSteps(10), TestNSteps(10), LoadWeights(false), InputStrength(1), 
DataDirectory("."), SimFinishButton(false), SaveInitialWeights(true)
{
    if (myAnyWrap) {
        std::string HomeDir(getenv("HOME"));
        std::string DefaultDataDir(HomeDir+"/data/tmp/");
        myAnyWrap->setOption("DataDirectory", 'a', &DataDirectory, DefaultDataDir.c_str());
        myAnyWrap->setOption("NTrials", 'a', &NTrials, 1);
        myAnyWrap->setOption("NSteps", 'a', &NSteps, 1);
        myAnyWrap->setOption("TestNSteps", 'a', &TestNSteps, 1);
        myAnyWrap->setFlag ( "LoadWeights", 'L', &LoadWeights);
        myAnyWrap->setOption( "InputStrength", 'i', &InputStrength, 1);
        myAnyWrap->setFlag ( "SimFinishButton", 'a', &SimFinishButton);
        myAnyWrap->setFlag("SaveInitialWeights", 'S', &SaveInitialWeights);
    }
}

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


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

// IDLrpc::IDLrpc()
// {
//    CLIENT *pClient;
//    char    cmdBuffer[512];
//    int     result;

// /* Modification by jimu 9-jan-2006 */
// char   *startupstr;
// char   startupcmd[100] = "@/home/frank/idl/test_startup";

// /*
//  * Connect to the server.
//  */
//    if((pClient = IDL_RPCInit(0, (char*)NULL)) == (CLIENT*)NULL){
//       fprintf(stderr, "Can't register with IDL server\n");
//       exit(1);
//    }

// printf("Get Environment Variable");fflush(stdout);
// startupstr = getenv("IDL_STARTUP");
// printf("%s\n",startupstr);fflush(stdout);
// strcat(startupcmd, startupstr);
// printf("IDL startup command: %s\n",startupcmd);fflush(stdout);
// result = IDL_RPCExecuteStr(pClient, startupcmd);
// result = IDL_RPCExecuteStr(pClient, "print, 'The value of !MyTestVar is:  '");
// result = IDL_RPCExecuteStr(pClient, "help, !MyTestVar");



// }

// IDLrpc::~IDLrpc()
// {

// }


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