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

#include "vconnection.hpp"
#include "vlearn.hpp"
#include "layer.hpp"
#include "matrix4d.hpp"
#include "chunkfile/chunkfile/chunkfile.h"
#include <algorithm>

const T_Delays DELETED_SYNAPSE=-1;

VecConnection::VecConnection():  NSynapses(0), learnobj(0), FileTypeString("VecConnection_1.0")
{

}

VecConnection::VecConnection(layer* SL, layer* TL, csimInputChannel _InputChannel, bool _nonself) 
     : Connection(SL,TL,_InputChannel,_nonself), NSynapses(0),
       learnobj(0), FileTypeString("VecConnection_1.0")
{
     Name="VecConnection";
     WeightFileName ="vweights.dat";
     // initialize all pointers with 0
     // delays(SourceLayer->N, maximumDelay, 2);
     delays.resize(ns);
     PreSynNr.resize(nt);
}

VecConnection::~VecConnection()
{
     Dout(dc::con, "VecConnection Destructor");fflush(stdout);
     if (learnobj != 0) delete learnobj;
}

VecConnectionInfo VecConnection::GetConnectionInfo()
{
     VecConnectionInfo info;
     info.Dmax = Dmax;
     info.maximumDelay = maximumDelay;
     info.MaxWeight = maxWeight;
     info.TargetLayer = TargetLayer;
     info.SourceLayer = SourceLayer;
     info.PSynWeights = &SynWeights;
     info.PSynSourceNr = &SynSourceNr;
     info.PSynTargetNr = &SynTargetNr;
     info.PSynDelays = &SynDelays;
     info.PPreSynNr = &PreSynNr;
     info.Pdelays = &delays;
         
     return info;
}



int VecConnection::proceede(int TotalTime)
{     
     int t = TotalTime % MacroTimeStep;
     int i,j,k, mi, ipre;
     int CurDelay;

     // calculate input for target layer  
     k=SourceLayer->N_firings;
//      if (rec) rec->record(dt*TotalTime, s[Observe_s][Observe_m]); 
//      if (BinRec) BinRec->record();

     while (t-(SourceLayer->firings[--k][0]) < maximumDelay)     // Nur Spikes, die nicht laenger als maximumDelay in der Vergangeheit liegen werden beruecksichtigt. firins-array besteht aus firings[SpikeNr][0]:zeitpunkte, firnigs[SpikeNr][1]:Neuronennummer
     {
      CurDelay = t-SourceLayer->firings[k][0];
      ipre = SourceLayer->firings[k][1];
      for (vector <T_NSynapses>::iterator it=delays[ipre][CurDelay].begin();it != delays[ipre][CurDelay].end(); ++it)
      {
           InputPointer[SynTargetNr[(*it)]] += SynWeights[(*it)];
      }
     }    
     if (learn == true) {
      learnobj->proceede(TotalTime);
     }
}

int VecConnection::prepare(int step)
{    
     SimElement::prepare(step);
     if ((learn == true) && (learnobj != 0)) {
      learnobj->prepare();
      if (RewireOn) {
           Rewire(RewireThreshold, RewireMaxConnectivity);
      }
      if (AutoSave) Save();   
     }   
}

int VecConnection::PushBackNewSynapse(int source, int target, float weight, int delay) 
{
     SynSourceNr.push_back(source);
     SynTargetNr.push_back(target);
     SynWeights.push_back(weight);
     SynDelays.push_back(delay);
     PreSynNr[target].push_back(NSynapses);     
     ++NSynapses;
}


int VecConnection::ScrambleSynTargets()
{
// Targets der Synapsen zufaellig vertauschen (fuer Reviewer)
     // Scramble Targets
     for (vector<T_NNeurons>::iterator it=SynTargetNr.begin(); it != SynTargetNr.end();++it) {
      (*it) =gsl_rng_uniform_int(gslr, nt);
     }

     // rebuild delays array and PreSynNr array
     SetupDelaysArray();
     SetupPreSynNrArray();     
}

///////////// Connect-Routinen um verschiedene Typen von Verbindungsmatrizen  zu erzeugen
// ConnectFull
// ConnectRandomIncomming
// ConnectSelf
// ConnectGaussian


int VecConnection::ConnectFull(float MaxWeight, float maxDelay, float minDelay, bool RandomWeights)
{
     Dout(dc::con, "Connection::ConnectFull");

     SetMinMaxDelay(maxDelay, minDelay);      
     Strength=MaxWeight;
     
    if (ns != nt) NonSelf = false;   // NonSelf makes sense only for layers with same size

     int i,j;
     for (i=0;i<ns;i++) {
      for (j=0;j<nt;j++) 
      {
           if ((!NonSelf) || (i != j)) {
            // if (RandomWeights)
            PushBackNewSynapse(i, j, MaxWeight, 0);
           }
      }
     }

     Dout(dc::con, "connected, NSynapses=" << NSynapses);
     SetRandomDelays();
     SetupDelaysArray();
//      SetupPresynapticInfo();

     Dout(dc::con, "Connection initialized");       
}

int VecConnection::ConnectRandomIncomming(double _Connectivity, float _InitialWeights, float maxDelay, float minDelay, bool RandomDelays)
{
     if ((_Connectivity <= 1) && (_Connectivity>=0)) {
      int NConnections = int(round(float(ns)*_Connectivity));
      ConnectRandomIncomming(NConnections, _InitialWeights, maxDelay, minDelay, RandomDelays);
     } else {
      cerr << "fatal ERROR in connection::ConnectRandom2: _Connectivity should be [0..1] (exit simulation)\n";
      exit (1);
     }
}


int VecConnection::ConnectRandomIncomming(int NIncommingCon, float _InitialWeights, float maxDelay, float minDelay, bool RandomDelays)
{
     Dout(dc::con, "VecConnection::ConnectRandomIncomming");

     SetMinMaxDelay(maxDelay, minDelay);
      

     if (ns != nt) NonSelf = false;   // NonSelf makes sense only for layers with same size

     int i;

     float InitialWeights=1;

     connectivity=float(NIncommingCon)/float(ns);

     vector <int> SourceList;
     SourceList.reserve(ns);
     for (int i=0;i<ns;++i) {
      SourceList.push_back(i);
     }
     long SynCounter=0;
     Dout(dc::con, "  nt=" << nt << "");
     Dout(dc::con, "  NIncommingCon=" << NIncommingCon << "");
     for (int tar=0;tar<nt;++tar) {
      vector <int> AvailSrcNeurons (SourceList);
      for (int IncomCon=0;IncomCon<NIncommingCon;++IncomCon) {
           int rndNum=gsl_rng_uniform_int(gslr, AvailSrcNeurons.size());
           int SourceNum=AvailSrcNeurons[rndNum];
           PushBackNewSynapse(SourceNum, tar, _InitialWeights, 0);
           AvailSrcNeurons.erase(AvailSrcNeurons.begin()+rndNum);
      }
     }
     Dout(dc::con, "* connected, NSynapses=" << NSynapses << ""); fflush(stdout);
     SetRandomDelays();
     SetupDelaysArray();

     Dout(dc::con, "  Tar0"<<PreSynNr[0].size()<<"");
     Dout(dc::con, "  Tar1"<<PreSynNr[1].size()<<"");
     Dout(dc::con, "  Tar2"<<PreSynNr[2].size()<<"");
     Dout(dc::con, "  NSynapses=" << NSynapses <<"");

     Dout(dc::con, "   VecConnection initialized\n");       
}

// verbinde jedes Neuron der Source-Schicht mit dem entsprechenden
// Neuron in der Target-Schicht, welches die gleiche Nummer hat
int VecConnection::ConnectSelf(float MaxWeight, float maxDelay, float minDelay)
{
    Strength=MaxWeight;
    Dout(dc::con, "VecConnection::ConnectSelf");
    
    SetMinMaxDelay(maxDelay, minDelay);
    
    int CurDelay;
    int n_connections = min(ns,nt);
    Dout(dc::con, "  number of connections =" << n_connections << "");

    // generate delays
    vector<T_Delays> Delays(n_connections, minimumDelay);  
    if (DelayDiff>=0) {
        Dout(dc::con, "generate Delase for each synapse");
        for (int i=0;i<n_connections;i++) {
            Delays[i] = minimumDelay + gsl_rng_uniform_int(gslr, DelayDiff+1);
        }
    }
    
    for (int i=0;i<n_connections;i++)
    {
        PushBackNewSynapse(i,i,MaxWeight, Delays[i]);
    }
    SetupDelaysArray();
    Dout(dc::con, "  VecConnection initialized\n"); fflush(stdout);
}




// ConnectGaussian: Verbinde Neuronen in Abhaengigkeit von ihrem Abstand 
// in einer 2D-Ebene, Verbindungsstaerke faellt Gauss-Foermig ab
int VecConnection::ConnectGaussian(float Sigma, float MaxWeight, float maxDelay, float minDelay, bool Cyclic)
{    
    Dout(dc::con, "VecConnection::ConnectGaussian");
    Strength=MaxWeight;

    vector<vector2d> TargetPositions;
    vector<vector2d> SourcePositions;
    try
    {
        TargetPositions = TargetLayer->getPositions();
        SourcePositions = SourceLayer->getPositions();
    }
    catch (std::bad_alloc) {
        cerr << "ERROR: Target or Source layer dead\n";
        throw;
    }
    
    Dout(dc::con, "sizeof TargetPositions=" << TargetPositions.size());
         
    if (ns != nt) NonSelf = false;   // NonSelf makes sense only for layers with same size
    
    int i,j,k,exists, r;
    
      // double loop (ns, nt) goes throug every possible synapse and checks wether to connect or not
    
    int CurDelay=0;
    SetMinMaxDelay(maxDelay, minDelay);
    
    float Distance;
    vector2d basis(1,1);
    SimpleTextProgressBar pgbar(ns);
    float MaxConDistance = 3*Sigma;
    for (i=0;i<ns;i++) {
        pgbar.Next(i);
        for (j=0;j<nt;j++)
        {
            
           // calculate Distance between neurons
            if (Cyclic) Distance = (TargetPositions[j]).CyclicDistance(SourcePositions[i], basis);
            else Distance = (TargetPositions[j] - SourcePositions[i]).abs() ;
            if ((Distance < MaxConDistance) && (!NonSelf || (i!=j))) {
                CurDelay =  minimumDelay + int(Distance*DelayDiff/MaxConDistance);
                if (CurDelay>=maximumDelay)
                {
                    CurDelay = maximumDelay;
                    cerr << "ERROR: delay too high\n";
                }
                PushBackNewSynapse(i,j,MaxWeight*gauss(Distance,Sigma), CurDelay);
            }
        }
    }
    SetupDelaysArray();
    Dout(dc::con, "  NSynapses=" << NSynapses);
    Dout(dc::con, "  connected"); fflush(stdout);
    Dout(dc::con, "  VecConnection initialized\n");
}


void VecConnection::ConnectCircular(const RangeConnectionParameters& Paras)
{
    Dout(dc::con, "VecConnection::ConnectCircular");
    
    vector<vector2d> TargetPositions;
    vector<vector2d> SourcePositions;
    try
    {
        TargetPositions = TargetLayer->getPositions();
        SourcePositions = SourceLayer->getPositions();
    }
    catch (std::bad_alloc) {
        cerr << "ERROR: Target or Source layer dead\n";
        throw;
    }
    
    if (!TargetLayer->isPositionNormalized() || !SourceLayer->isPositionNormalized()) {
        throw std::runtime_error("Source or Target positions not normalized");
    }
    
    if (ns != nt) NonSelf = false;   // NonSelf makes sense only for layers with same size
       
    int CurDelay=0;
    SetMinMaxDelay(Paras.MaxDelay, Paras.MinDelay);
    
    float Distance;
    vector2d basis(1,1);
    SimpleTextProgressBar pgbar(ns);
          
    // double loop (ns, nt) goes throug every possible synapse
    // and checks wether to connect or not    
    for (int i=0;i<ns;i++) {
        pgbar.Next(i);
        for (int j=0;j<nt;j++)
        {
           // calculate Distance between neurons
            if (Paras.Cyclic) {
                Distance = (TargetPositions[j]).CyclicDistance(SourcePositions[i], basis);
            } else {
                Distance = (TargetPositions[j] - SourcePositions[i]).abs() ;
            }
            
            if ((Distance < Paras.Range) && (!NonSelf || (i!=j))) {
                if (Paras.Connectivity>=1 || (gsl_rng_uniform(gslr) <= Paras.Connectivity)) {
                    CurDelay =  minimumDelay + int(Distance*DelayDiff/Paras.Range);
                    if (CurDelay>=maximumDelay)
                    {
                        CurDelay = maximumDelay;
                        cerr << "ERROR: delay too high\n";
                    }
                    PushBackNewSynapse(i,j,Paras.Strength, CurDelay);
                }
            }
        }
    }
    
    SetupDelaysArray();
    Dout(dc::con, "  NSynapses=" << NSynapses);
    Dout(dc::con, "  connected"); 
    Dout(dc::con, "  VecConnection initialized\n");    
}

///////////// ENDE Connect-Routinen //////////////////




int VecConnection::SetNewWeights(float IncommingConnectivity, float InitialWeights)
{
    Dout(dc::con, "VecConnection::SetNewWeights"); fflush(stdout);
    vector <int> SourceNeuronList (ns);
    for (int s=0;s<ns;++s) SourceNeuronList[s]=s;
    
    for (int tar=0; tar<nt;++tar) {
        int NNewWeights = int(ns*IncommingConnectivity) - PreSynNr[tar].size();
//        Dout(dc::con, "NNewWeights=" << NNewWeights << ""); fflush(stdout);
//        Dout(dc::con, "ns" << ns << ""); fflush(stdout);
//        Dout(dc::con, "IncommingConnectivity=" << IncommingConnectivity << ""); fflush(stdout);
//        Dout(dc::con, "PreSynNr[tar].size()=" << PreSynNr[tar].size() << ""); fflush(stdout);
        if (NNewWeights >0) {
            vector <int> AvailNeur(SourceNeuronList);
           // remove all presynaptic Neurons from available neurons list
            for (vector<T_NSynapses>::iterator it=PreSynNr[tar].begin();it!=PreSynNr[tar].end();++it) {
                AvailNeur.erase(remove(AvailNeur.begin(), 
                                       AvailNeur.end(), 
                                       SynSourceNr[*it]), 
                                AvailNeur.end());
            }
            for (int nw=0;nw<NNewWeights;++nw) {
            // Dout(dc::con, "AvailNeur=" << AvailNeur.size() << ""); fflush(stdout);
                int rndNum=gsl_rng_uniform_int(gslr, AvailNeur.size());
                int SourceNum=AvailNeur[rndNum];
                T_NSynapses NewSynNr = SynSourceNr.size();
                SynSourceNr.push_back(SourceNum);
                SynTargetNr.push_back(tar);
                SynWeights.push_back(InitialWeights);
            // random delay??       
                SynDelays.push_back( minimumDelay + gsl_rng_uniform_int(gslr, DelayDiff+1)); 
            //! @ToDo: Check wether +1 is correct
                PreSynNr[tar].push_back(NewSynNr);
                delays[SynSourceNr[NewSynNr]][SynDelays[NewSynNr]].push_back(NewSynNr);
                ++NSynapses;
                AvailNeur.erase(AvailNeur.begin()+rndNum);
            }
        }
    }
}

// vector<int>* _NNewWeights gibt fuer jedes Target-Neuron an, wieviele neue eingehende
// Synapsen gesetzt werden
int VecConnection::SetNewWeights(vector<int>* _NNewWeights, float InitialWeights)
{
     Dout(dc::con, "VecConnection::SetNewWeights");
     vector <int> SourceNeuronList (ns);
     for (int s=0;s<ns;++s) SourceNeuronList[s]=s;
     
     for (int tar=0; tar<nt;++tar) {
      int NNewWeights = (*_NNewWeights)[tar];
      // Dout(dc::con, "NNewWeights=" << NNewWeights << "");fflush(stdout);
      if (NNewWeights >0) {
           vector <int> AvailNeur(SourceNeuronList);
           // remove all presynaptic Neurons from available neurons list
           for (vector<T_NSynapses>::iterator it=PreSynNr[tar].begin();it!=PreSynNr[tar].end();++it) {
            AvailNeur.erase(remove(AvailNeur.begin(),AvailNeur.end(),SynSourceNr[*it]), AvailNeur.end());
           }
           for (int nw=0;nw<NNewWeights;++nw) {
            Dout(dc::con, "AvailNeur=" << AvailNeur.size() << ""); fflush(stdout);
            int rndNum=gsl_rng_uniform_int(gslr, AvailNeur.size());
            int SourceNum=AvailNeur[rndNum];
            int NewSynNr = SynSourceNr.size();
            SynSourceNr.push_back(SourceNum);
            SynTargetNr.push_back(tar);
            SynWeights.push_back(InitialWeights);
            // random delay??       
            SynDelays.push_back( minimumDelay + gsl_rng_uniform_int(gslr, DelayDiff+1)); 
            //ToDo: Check wether +1 is correct
            PreSynNr[tar].push_back(NewSynNr);
            delays[SynSourceNr[NewSynNr]][SynDelays[NewSynNr]].push_back(NewSynNr);
            ++NSynapses;
            AvailNeur.erase(AvailNeur.begin()+rndNum);
           }
      }
     }
}


/** calculate maximumDelay and minimumDelay (in time steps) from maxDelay, minDelay (in mm/s)
    @param maxDelay maximum delay in mm/s
    @param minDelay minimum delay in mm/s
*/
int VecConnection::SetMinMaxDelay(float maxDelay, float minDelay)
{
    maximumDelay = int(maxDelay/dt);
    if (maximumDelay == 0) maximumDelay = 1;
     // equivalent to Dmax; compare allways with "<" (not with "<=")
    minimumDelay = int(minDelay/dt);
    if (maximumDelay >= Dmax) {
        cerr << "ERROR: SetMinMaxDelay: Parameter Error: maximumDelay="
            << (int) maximumDelay << " > Dmax=" << (int) Dmax << "\n";
        throw RequestedDelayTooLarge(maximumDelay);        
    }
    if (minimumDelay >= maximumDelay) {
        cerr << "ERROR: SetMinMaxDelay: Parameter Error:"
            << " minimumDelay=" << (int)minimumDelay << " >= maximumDelay =" << (int)maximumDelay << "\n";
        minimumDelay = maximumDelay-1;
        cerr << "setting minimumDelay to " << (int)minimumDelay << "\n";
    }
    
    DelayDiff = maximumDelay-minimumDelay-1;
    Dout(dc::con, "MaximumDelay=" << (int)maximumDelay
         << "  MinimumDelay=" << (int)minimumDelay
         << "  DelayDiff= " << (int)DelayDiff);
}




int VecConnection::SetRandomDelays()
{
     Dout(dc::con, "SetRandomDelays"); fflush(stdout);
     int NSynapses = SynWeights.size();
     for (int SynNr=0;SynNr<NSynapses;++SynNr) {
      // SynDelays[SynNr] = minimumDelay + getrandom(DelayDiff+1);
      SynDelays[SynNr] = minimumDelay + gsl_rng_uniform_int(gslr, DelayDiff+1); 
     }
}

int VecConnection::SetupPreSynNrArray()
{
    Dout(dc::con, "VecConnection::SetupPreSynNrArray");
    Dout(dc::con, "PreSynNr.size()=" << PreSynNr.size()
         << " NSynapses=" << NSynapses << " size()=" << SynDelays.size());

    if (PreSynNr.size() != nt) {
        PreSynNr.resize(nt);
    }
    
    if (SynWeights.size() != NSynapses) {
        cerr << "ERROR in SetupPreSynNrArray(): there are deleted Synapses \n";
        cerr << "              or did you forget 'NSynapses=SynCounter;' in your Connect routine?\n";
        exit(1);
    }
     // clear PreSynNrArray
    for (vector< vector<T_NSynapses> >::iterator it=PreSynNr.begin();it!=PreSynNr.end();++it) {
        (*it).clear();
    }
    
    for (T_NSynapses i=0;i<NSynapses;++i) {
        if (SynDelays[i] == -1) {
            cerr << "ERROR in SetupPreSynNrArray(): there are deleted Synapses \n";
            exit(1);
        } else {
            PreSynNr[SynTargetNr[i]].push_back(i);
        }
    }
}

 
int VecConnection::SetupDelaysArray()
{
    Dout(dc::con, "SetupDelaysArray"); fflush(stdout);
    
    if (SynWeights.size() != NSynapses) {
        cerr << "ERROR in SetupDelaysArray(): there are deleted Synapses \n";
        exit(1);
    }

    delays.resize(ns);
    
    Dout(dc::con, "maximumDelay=" << maximumDelay << "");
    for (int i=0;i<ns;++i) delays[i].resize(maximumDelay);
    
    Dout(dc::con, "Clear delays Array" << "");
    for (vector<vector<vector<T_NSynapses> > >::iterator ita=delays.begin();ita!=delays.end();++ita) {
        for (vector<vector<T_NSynapses> >::iterator itb=(*ita).begin();itb!=(*ita).end();++itb) {
            (*itb).clear();
        }
    }
    
#ifdef DEBUG
    int CurDelay=0;
    int CurSource=0;
    Dout(dc::con, "DEBUGGING mode ");
    for (int SynNr=0;SynNr<NSynapses;++SynNr) {
        CurDelay = SynDelays[SynNr];
        CurSource=SynSourceNr[SynNr];
        if (CurDelay < delays[CurSource].capacity()) {
            delays[CurSource][CurDelay].push_back(SynNr);
        } else {
            cerr << "ERROR in SetupDelays() \n";
            exit (1);
        }
    }
#else //DEBUG
    Dout(dc::con, "Rebuild delays Array" << ""); fflush(stdout);
    
    for (int SynNr=0;SynNr<NSynapses;++SynNr) {
//    cout << "[" 
//         << CountSyn++ << ":" 
//         << delays[SynSourceNr[SynNr]][SynDelays[SynNr]].size() << "," 
//         << SynSourceNr[SynNr] << "," 
//         << SynTargetNr[SynNr] << "," 
//         << SynDelays[SynNr] << ",c" 
//         << delays[SynSourceNr[SynNr]][SynDelays[SynNr]].capacity() <<  "] .\n "; fflush(stdout);
        if (SynDelays[SynNr] >= 0) {
            delays[SynSourceNr[SynNr]][SynDelays[SynNr]].push_back(SynNr);
        } else {
            cerr << "WARNING in SetupDelays(): there are deleted synapses,  \n";
        }
    }
#endif //DEBUG
    Dout(dc::con, "done rebuilding delays Array" << ""); fflush(stdout);
}

/**
 * delete a single synapse
 * @param SynNr is the number of the synapse (index in SynSourceNr, 
 *              SynTargetNr, SynDelays, SynWeights
 * @return 
 */
int VecConnection::DeleteSynapse(int SynNr)
{
     // Dout(dc::con, "Deleting Synapse Nr:" << SynNr << " weight=" << SynWeights[SynNr] << "  TargetNr=" << SynTargetNr[SynNr] << "");
     int DelSourceNr=SynSourceNr[SynNr];
     int DelTargetNr = SynTargetNr[SynNr];
     int DelDelay = SynDelays[SynNr];

     if (DelDelay <0) {
      cerr << "Deleting Synapse which is already deleted\n";
      exit (1);
     }

     // delete synapse in delays vector
     for (vector <T_NSynapses>::iterator it=delays[DelSourceNr][DelDelay].begin(); it != delays[DelSourceNr][DelDelay].end(); ++it) 
     {
      if (*it == SynNr) {
           delays[DelSourceNr][DelDelay].erase(it);
           break;
      }
     }
    

     // delete synapse in PreSynNr vector
     for (vector<T_NSynapses>::iterator it=PreSynNr[DelTargetNr].begin(); it != PreSynNr[DelTargetNr].end(); ++it) {
      if (*it == SynNr) {
           PreSynNr[DelTargetNr].erase(it);
           break;
      }
     }
     
     // set synapse to -1
     SynSourceNr[SynNr] = -1;
     SynTargetNr[SynNr] = -1;
     SynDelays[SynNr] = -1;
     SynWeights[SynNr] = -1;
     --NSynapses;
}


/**
 * set all valid synapses to new WeightValue
 * @param WeightValue new weight value
 * @return 
 */
int VecConnection::SetAllWeights(float WeightValue)
{
    int NWeights=SynWeights.size();
    for (int i=0; i<NWeights; ++i) {
        if (SynDelays[i] != -1) { // invalid (deleted) synapses have delay -1
            SynWeights[i]=WeightValue;
        }
    }
}

int VecConnection::DeleteLowWeights(float Threshold)
{
     vector <long> DeleteList;
     for (int SourceNr=0; SourceNr<ns;++SourceNr) {
      int NDelays=delays[SourceNr].size();
      for (int Delay=0;Delay<NDelays; ++Delay) {
           for (vector <T_NSynapses>::iterator it=delays[SourceNr][Delay].begin(); it != delays[SourceNr][Delay].end(); ++it) {
            if (SynWeights[*it ] < Threshold) {
             DeleteList.push_back(*it);
             if (SynDelays[*it] <0) {
                  cerr << "Hier stimmt was nicht\n";
                  cerr << "SynNr=" << (*it) << "\n";
                  exit (1);
             }
            }
           }
      }
     }
     
     Dout(dc::con, "DeleteSynapses");fflush(stdout);
     for (vector <long>::iterator it=DeleteList.begin(); it != DeleteList.end(); ++it) {
      DeleteSynapse(*it);
     }
     int NDeletedSynapses = DeleteList.size();
     Dout(dc::con, "Deleted " << NDeletedSynapses << "synapses");

     Dout(dc::con, "Done DeleteSynapses");fflush(stdout);
     if (SynWeights.size() > 20*NSynapses) {
      CleanupArrays();
     }
     Dout(dc::con, "DeleteThreshold = " << Threshold << "");
     return NDeletedSynapses;
}

int VecConnection::CleanupArrays()
{
     Dout(dc::con, "Renewing delays array and PreSynNr array");
     // remove deleted synapses
     int c=0;
     while (c < SynDelays.size()) {
      if (SynDelays[c] == -1) {
           SynDelays.erase(SynDelays.begin()+c);
           SynSourceNr.erase(SynSourceNr.begin()+c);
           SynTargetNr.erase(SynTargetNr.begin()+c);
           SynWeights.erase(SynWeights.begin()+c);
      } else {
           ++c;
      }
     }    
     SetupDelaysArray();
     SetupPreSynNrArray();     
}


void VecConnection::SetName(const char* _name)
{
     SimElement::SetName(_name);
     WeightFileName = Name + "weights.dat";     
     Dout(dc::con, "VecConnection, WeightFileName = " << WeightFileName << "");
}


int VecConnection::Save()
{
     Save(WeightFileName.c_str());
}

int VecConnection::Save(int nr)
{
     Save((WeightFileName + stringify(nr)).c_str());
}

void VecConnection::Save_VecConnection_2_0(const string& SaveWeightFileName)
{
    FileFormat VecFileFormat("VecConnection", 2, 0);
    ChunkFileWriter MyFileWriter(SaveWeightFileName, VecFileFormat);
    SourceTargetDim MyVecConHeader(ns,
                                   SourceNx,
                                   SourceNy,
                                   nt,
                                   TargetNx,
                                   TargetNy);
				   
    
    MyFileWriter.write("VecHeader", MyVecConHeader);    
    MyFileWriter.writeVector("SourceNr", SynSourceNr);
    MyFileWriter.writeVector("TargetNr", SynTargetNr);
    MyFileWriter.writeVector("Delays", SynDelays);
    MyFileWriter.writeVector("Weights", SynWeights);
    
    MyFileWriter.close();
}

void VecConnection::Load_VecConnection_2_0(const string & FileName)
{
    Dout(dc::con, "VecConnection::Load_VecConnection_2_0, FileName=" << FileName);
    ChunkFileReader FileReader(FileName, FileFormat("VecConnection", 2, 0));
    Dout(dc::con, "allocate SourceTargetDim");
    SourceTargetDim VecConHeader;
    Dout(dc::con, "read VecHeader");
    FileReader.read("VecHeader", VecConHeader);
    ns = VecConHeader.NSource;
    nt = VecConHeader.NTarget;
    SourceNx = VecConHeader.NSx;
    SourceNy = VecConHeader.NSy;
    TargetNx = VecConHeader.NTx;
    TargetNy = VecConHeader.NTy;

    Dout(dc::con, "read vectors");
    FileReader.readAndCastVector("SourceNr", SynSourceNr);
    FileReader.readAndCastVector("TargetNr", SynTargetNr);
    FileReader.readAndCastVector("Delays", SynDelays);
    FileReader.readVector("Weights", SynWeights);

    FileReader.close();

    Dout(dc::con, "Load_VecConnection_2_0: Data loaded");
    //Dout(dc::con, "Delays=" << stringifyVector(OriginalDelays));
    vector<T_Delays>::iterator it = max_element(SynDelays.begin(),SynDelays.end());
    maximumDelay = (*it)+1;
    if (maximumDelay >= Dmax) {
        throw RequestedDelayTooLarge(maximumDelay);
    }
    updateNSynapses();
    SetupDelaysArray();
    SetupPreSynNrArray();
}

void VecConnection::Load_VecConnection_2_1(const string & FileName)
{
    Dout(dc::con, "VecConnection::Load_VecConnection_2_1, FileName=" << FileName);
    ChunkFileReader FileReader(FileName, FileFormat("VecConnection", 2, 1));
    Dout(dc::con, "allocate SourceTargetDim");
    SourceTargetDim1 VecConHeader;
    Dout(dc::con, "read VecHeader");
    FileReader.read("VecHeader", VecConHeader);
    ns = VecConHeader.NSource;
    nt = VecConHeader.NTarget;
    SourceNx = VecConHeader.NSx;
    SourceNy = VecConHeader.NSy;
    TargetNx = VecConHeader.NTx;
    TargetNy = VecConHeader.NTy;
    ArrayOrderXFast = VecConHeader.ArrayOrderXFast == 1;
    Dout(dc::con, "XFast=" << ArrayOrderXFast);
    
    Dout(dc::con, "read vectors");
    FileReader.readAndCastVector("SourceNr", SynSourceNr);
    FileReader.readAndCastVector("TargetNr", SynTargetNr);
    FileReader.readAndCastVector("Delays", SynDelays);
    FileReader.readVector("Weights", SynWeights);
    
    FileReader.close();

    Dout(dc::con, "Load_VecConnection_2_1: Data loaded");
    //Dout(dc::con, "Delays=" << stringifyVector(OriginalDelays));
    vector<T_Delays>::iterator it = max_element(SynDelays.begin(),SynDelays.end());
    maximumDelay = (*it)+1;
    if (maximumDelay >= Dmax) {
        throw RequestedDelayTooLarge(maximumDelay);
    }
    updateNSynapses();
    SetupDelaysArray();
    SetupPreSynNrArray();
}


int VecConnection::Save(const string& SaveWeightFileName)
{
    string DirAndFileName = DataDirectory+SaveWeightFileName;
    Save_VecConnection_2_0(DirAndFileName);
}
 
void VecConnection::Save_VecConnection_1_0(const string& DirAndFileName)
{
    if (compiledWithMemsave()) {
          throw NotSupportedWithMemsave("Save_VecConnection_1_0 not available when compiled with -sMEMSAVE=1");
    }
    
    int i,j,k;
    Dout(dc::con, "Save VecCon.file: " << DirAndFileName);
    FILE *fw;
     // save weights
    fw = fopen(DirAndFileName.c_str(),"w");
    
     // save file header
    const char *fts = FileTypeString.c_str();
    char ftsLength = strlen(fts);
    fwrite(&ftsLength, sizeof(ftsLength), 1, fw);
    fwrite(fts, ftsLength, 1, fw);
    
    Dout(dc::con, "SourceNx=" << SourceNx << "  SourceLayer->Nx=" << SourceLayer->Nx);
    
     // save Synapses
    fwrite(&ns, sizeof(ns), 1, fw);
    fwrite(&(SourceNx), sizeof(SourceNx), 1, fw);
    fwrite(&(SourceNy), sizeof(SourceNy), 1, fw);
    fwrite(&nt, sizeof(nt), 1, fw);
    fwrite(&(TargetNx), sizeof(TargetNx), 1, fw);
    fwrite(&(TargetNy), sizeof(TargetNy), 1, fw);
    fwrite(&maximumDelay, sizeof(maximumDelay), 1, fw);
    fwrite(&NSynapses, sizeof(NSynapses), 1, fw);
    int SynSize = SynWeights.size();
    
    for (int SynNr=0;SynNr<SynSize;++SynNr) {
        if (SynSourceNr[SynNr] != -1) {
            fwrite(&SynWeights[SynNr], sizeof(SynWeights[SynNr]), 1, fw);
            fwrite(&SynSourceNr[SynNr], sizeof(SynSourceNr[SynNr]), 1, fw);
            fwrite(&SynTargetNr[SynNr], sizeof(SynTargetNr[SynNr]), 1, fw);
            fwrite(&SynDelays[SynNr], sizeof(SynDelays[SynNr]), 1, fw);
        } // else {Dout(dc::con, "not saving deleted Syn");}
    }
    
    fclose(fw);
    Dout(dc::con, "  saved weights ");
}

int VecConnection::DeleteSynapseArrays()
{
    SynWeights.clear();
    SynSourceNr.clear();
    SynTargetNr.clear();
    SynDelays.clear();
    for (vector<vector<T_NSynapses> >::iterator it=PreSynNr.begin(); it!=PreSynNr.end(); ++it) {
        (*it).clear();  
    }
    for (vector<vector<vector<T_NSynapses> > >::iterator ita=delays.begin(); ita!=delays.end(); ++ita) {
        for (vector< vector<T_NSynapses> >::iterator itb=(*ita).begin();itb!=(*ita).end();++itb) {
            (*itb).clear();
        }
    }
    NSynapses=0;maximumDelay=1;
}

int VecConnection::Load()
{
     return Load(WeightFileName);
}

int VecConnection::Load(const string& FileName)
{
     return Load(FileName, DataDirectory);
}
 

int VecConnection::ReserveSynapses(int _nsynapses) 
{
     SynWeights.reserve(_nsynapses);
     SynTargetNr.reserve(_nsynapses);
     SynSourceNr.reserve(_nsynapses);
     SynDelays.reserve(_nsynapses);
}

/** @brief load VecConnection from data file
@todo design new weight file format, especially with respect
to different data types for Source and Target neurons (T_Delays) and delays (T_Delays)

hints for file format design:
http://www.magicdb.org/filedesign.html
http://decoy.iki.fi/texts/filefd/filefd
*/
int VecConnection::Load(const string& FileName, const string& DirName)
{
    std::string DirAndFileName =  DirName+FileName;
    Dout(dc::con, "DirAndFileName=" << DirAndFileName << "");
    
    const char* DFileName = DirAndFileName.c_str();
    if (!fexist(DFileName)) {
        cerr << "\n\nERROR: connection file " <<  DFileName << " doesn't exist \n\n";
        fflush(stderr);
        return(2);
    } 
        
    Dout(dc::con, "LoadWeightFile: " << DFileName << "  ");
    FileFormat WeightFileFormat = readFileFormat(DirAndFileName);
    Dout(dc::con, "FileFormat = " << WeightFileFormat.print());
    
    //cout << "FileFormat = " << WeightFileFormat.print() << "\n";
    if (WeightFileFormat.isEqual("VecConnection_1.0")) {
        Dout(dc::con, "Load VecConnection_1.0");
        return Load_VecConnection_1_0(DirAndFileName);
    } else if (WeightFileFormat.isEqual(FileFormat("VecConnection",2,0))) {
        Dout(dc::con, "Load VecConnection_2.0");
        Load_VecConnection_2_0(DirAndFileName);
        return 0;
    } else if (WeightFileFormat.isEqual(FileFormat("VecConnection",2,1))) {
        Dout(dc::con, "Load VecConnection_2.1");
        Load_VecConnection_2_1(DirAndFileName);
        return 0;
    }
    return -1;
}


int VecConnection::Load_VecConnection_1_0(const string& FileNameWithDir)
{
    if (compiledWithLowMemConfig()) {
        throw NotSupportedWithMemsave("Load_VecConnection_1_0 not available when compiled with -sMEMSAVE=1");
    }
    int i,j,k;
    FILE *fw;
    int dummy;
    fw = fopen(FileNameWithDir.c_str(), "r");
      
    // Load file header
    unsigned char ftsLength;
    fread(&ftsLength, sizeof(ftsLength), 1, fw);
    int FLength = static_cast<int>(ftsLength);
    Dout(dc::con, "StrLen of file type = " << FLength << "");
    char *LoadedFileType = new char [ftsLength+1];
    fread(LoadedFileType, ftsLength, 1, fw);
    LoadedFileType[ftsLength]=0;
    Dout(dc::con, "loaded file type:" << LoadedFileType << "");
    if (strcmp(LoadedFileType, FileTypeString.c_str())) {
        cerr << "VecConnection.load() tried to load wrong file type\n";
        cerr << "file type was: " << LoadedFileType <<"\n";
        cerr << "file type should have been: " << FileTypeString << "\n";
        return(-1);
    };
    delete[]LoadedFileType;
    DeleteSynapseArrays();
    
    fread(&ns, sizeof(ns), 1, fw);
    Dout(dc::con, "SourceN=" << ns << ""); 
    fread(&SourceNx, sizeof(SourceNx), 1, fw);
    fread(&SourceNy, sizeof(SourceNy), 1, fw);
    fread(&nt, sizeof(nt), 1, fw);
    Dout(dc::con, "TargetN=" << nt << "");
    fread(&TargetNx, sizeof(TargetNx), 1, fw);
    fread(&TargetNy, sizeof(TargetNy), 1, fw);
    fread(&maximumDelay, sizeof(maximumDelay), 1, fw);
    Dout(dc::con, "maximumDelay=" << maximumDelay << ""); 
    Dout(dc::con, "ConDimensions= (" << SourceNx << ", " << SourceNy << ", " << TargetNx << ", " << TargetNy << ")");
    
    if (CheckHeaderConsistency()==false)  {
        return(-1);
    }
    
    delays.resize(ns);
    PreSynNr.resize(nt);
    
    int _nsynapses;
    fread(&_nsynapses, sizeof(_nsynapses), 1, fw);
    Dout(dc::con, "NSynapses=" << _nsynapses<< "");
    ReserveSynapses(_nsynapses);
    Dout(dc::con, "Synapses reserved ");
    float w;
    int tar,src,del;
    for (int i=0;i<_nsynapses;++i) {
        fread(&w,sizeof(w),1,fw);
        fread(&src,sizeof(src),1,fw);
        fread(&tar,sizeof(tar),1,fw);
        fread(&del,sizeof(del),1,fw);
        PushBackNewSynapse(src,tar,w,del);
        if (w>0) {
              // Dout(dc::con, "(x,src,tar,del)=(" <<w<<","<<src<<","<<tar<<","<<del<<")");
        }
    }
    if (_nsynapses != NSynapses) {
        cerr << "ERROR in VecConnection::Load\n";
        exit(1);
    }
    SetupDelaysArray();
    SetupPreSynNrArray();
    Dout(dc::con, "finished loading connection");
    
    return(0);
    

}



int VecConnection::MultiplyTargetWeights(int TargetNr, float Factor)
{
    for (vector<T_NSynapses>::iterator it=PreSynNr[TargetNr].begin(); it!=PreSynNr[TargetNr].end(); ++it) {
        SynWeights[(*it)] *= Factor;
    }
}

int VecConnection::MultiplyAllTargetWeights(float Factor)
{
     for (int t=0; t<nt;++t) MultiplyTargetWeights(t, Factor);
}


float VecConnection::GetWeightSum(int TargetNr, bool quadratic)
{
    float CurWeightSum=0;
    for (vector<T_NSynapses>::iterator it=PreSynNr[TargetNr].begin(); it!=PreSynNr[TargetNr].end(); ++it) {
        CurWeightSum += SynWeights[(*it)];
    }
    return CurWeightSum;
}

/**
 * calculate mean of all weights
 *
 * @return MeanWeight;  or -1 if there are no synapses
 */
float VecConnection::GetMeanWeight()
{
    int NumberOfSynapses = SynDelays.size();
    int SynCount=0;
    float WeightSum=0;
    for (int SynNr=0;SynNr<NumberOfSynapses;++SynNr) {
        if (SynDelays[SynNr] != -1) { // deleted synapses have dalay -1
            WeightSum+=SynWeights[SynNr];
            ++SynCount;
        }
    }
    assert(SynCount==NSynapses);
    float MeanWeight=-1;
    if (SynCount>0) {
        MeanWeight=WeightSum/SynCount;
    }
    return MeanWeight;
}

int VecConnection::WriteSimInfo(fstream &fw)
{
     stringstream sstr;
     sstr << "<Source id=\"" << SourceLayer->IdNumber << "\"/> \n";
     sstr << "<Target id=\"" << TargetLayer->IdNumber << "\"/> \n";
     sstr << "<MaxWeight value=\"" << maxWeight << "\"/> \n";     
     sstr << "<Strength value=\"" << Strength << "\"/> \n";
     sstr << "<InputNumber value=\"" << int(InputChannel) << "\"/> \n";
     if (learnobj) sstr << learnobj->GetSimInfo();
     SimElement::WriteSimInfo(fw, sstr.str());
}

int VecConnection::WriteSimInfo(fstream &fw, const string &ChildInfo)
{
     stringstream sstr;
     sstr << "<Source id=\"" << SourceLayer->IdNumber << "\"/> \n";
     sstr << "<Target id=\"" << TargetLayer->IdNumber << "\"/> \n";
     sstr << ChildInfo;
     sstr << "<MaxWeight value=\"" << maxWeight << "\"/> \n";     
     sstr << "<Strength value=\"" << Strength << "\"/> \n";
     sstr << "<InputNumber value=\"" << int(InputChannel) << "\"/> \n";
     if (learnobj) sstr << learnobj->GetSimInfo();
     SimElement::WriteSimInfo(fw, sstr.str());
}


int VecConnection::Rewire(float minWeight, float maxConnectivity)
{
// first: delete low weights

     Dout(dc::con, "DeleteLowWeights"); fflush(stdout);
     int NFreeWeights = DeleteLowWeights(minWeight);     

// second: set new weights randomly

     Dout(dc::con, "SetNewWeights"); fflush(stdout);
     SetNewWeights(maxConnectivity, InitialWeight);

}


int VecConnection::SetLearn(bool l)
{
     if (learnobj != 0) learn = l;
}

int VecConnection::SetLearnObj(veclearning* lo)
{
     learnobj = lo;
     SetLearn(true);
}

float VecConnection::GetSourceWeights(int SourceNr, vector<float>& WeightMatrix)
{
    WeightMatrix.clear();
    WeightMatrix.resize(nt);
    float MaxWeight=0;
    for(int CurDelay=0;CurDelay<maximumDelay;++CurDelay)
    {
        Dout(dc::con, "(Delay=" << CurDelay << ")  ");
        for (vector <T_NSynapses>::iterator it=delays[SourceNr][CurDelay].begin();
             it != delays[SourceNr][CurDelay].end();++it) 
        {
            int CurTarget=SynTargetNr[(*it)];
            Dout(dc::con, "(Target=" << CurTarget << ")  ");
            float CurWeight=SynWeights[(*it)];
            WeightMatrix[CurTarget]=CurWeight;
            if (CurWeight>MaxWeight) {
                MaxWeight=CurWeight;
            }
        }
    }
    return MaxWeight;
}

/** returns list of synapses form source neuron
 * @param [IN] SourceNr 
 * @param [OUT] SynList is a return vector containing a list of all synapses from source neuron SourceNr
 * @return 
 */
float VecConnection::GetSourceWeights(int SourceNr, vector<Synapse>& SynList, int &MaxDelay)
{
    SynList.clear();
    float MaxWeight=0;
    MaxDelay=0;
    for(int CurDelay=0;CurDelay<maximumDelay;++CurDelay)
    {
        for (vector <T_NSynapses>::iterator it=delays[SourceNr][CurDelay].begin();
             it != delays[SourceNr][CurDelay].end();++it) 
        {           
            int CurTarget=SynTargetNr[(*it)];
            int CurSource=SynSourceNr[(*it)];
            if (SourceNr!=CurSource) {
                cerr << "*******ERROR*****************************************\n";
                cerr << "float GetSourceWeights(int SourceNr, vector<Synapse>& SynList)\n";
                cerr << "korrupt\n";
            }
            float CurWeight=SynWeights[(*it)];
            SynList.push_back(Synapse(CurSource, CurTarget, CurWeight, CurDelay));
            if (CurWeight>MaxWeight) {
                MaxWeight=CurWeight;
            }
            if (CurDelay>MaxDelay) {
                MaxDelay=CurDelay;
            }
        }
    }
    return MaxWeight;   
}


float VecConnection::GetTargetWeights(int TargetNr, vector<float>& WeightMatrix)
{
    WeightMatrix.clear();
    WeightMatrix.resize(ns);
    float MaxWeight=0;
    for (vector<T_NSynapses>::iterator it=PreSynNr[TargetNr].begin();
         it!=PreSynNr[TargetNr].end();++it)
    {
        int CurSynNr=(*it);
        float CurWeight=SynWeights[CurSynNr];
        int SourceNr=SynSourceNr[CurSynNr];
        Dout(dc::con, "SynNr=" << CurSynNr << " SourceNr=" << SourceNr << " CurWeight=" << CurWeight << "");
        WeightMatrix[SourceNr]=CurWeight;
        if (CurWeight>MaxWeight) {
            MaxWeight=CurWeight;
        }
    }
    return MaxWeight;
}

float VecConnection::GetTargetWeights(int TargetNr, vector<Synapse>& SynList, int &MaxDelay)
{
    SynList.clear();
    float MaxWeight=0;
    MaxDelay=0;
    for (vector<T_NSynapses>::iterator it=PreSynNr[TargetNr].begin();
         it!=PreSynNr[TargetNr].end();++it)
    {
        int CurSynNr=(*it);
        float CurWeight=SynWeights[CurSynNr];
        int CurSource=SynSourceNr[CurSynNr];
        int CurTarget = SynTargetNr[CurSynNr];
        int CurDelay = SynDelays[CurSynNr];
        
        if (CurTarget!=TargetNr) {
            cerr << "*******ERROR*****************************************\n";
            cerr << "float GetTargetWeights(int TargetNr, vector<Synapse>& SynList)\n";
            cerr << "korrupt\n";
            
        }
        SynList.push_back(Synapse(CurSource, CurTarget, CurWeight, CurDelay));
        if (CurWeight>MaxWeight) {
            MaxWeight=CurWeight;
        }
        if (CurDelay>MaxDelay) {
            MaxDelay=CurDelay;
        }
    }
    return MaxWeight;
    
}


void VecConnection::GetWeightMatrix4D(Matrix4D& matrix)
{
    matrix.Resize(SourceNx, SourceNy, TargetNx, TargetNy);
    matrix.SetValue(-1);
    for (vector<vector<vector<T_NSynapses> > >::iterator DelSource=delays.begin();
         DelSource!=delays.end();++DelSource) {
             for (vector<vector<T_NSynapses> >::iterator Delay=(*DelSource).begin(); 
                  Delay!=(*DelSource).end();++Delay) {
                      for (vector<T_NSynapses>::iterator Index=(*Delay).begin();
                           Index!=(*Delay).end(); ++Index) {
                               int SourceNr=SynSourceNr[(*Index)];
                               int TargetNr=SynTargetNr[(*Index)];
                               int SourceY = SourceNr/SourceNx;
                               int SourceX = SourceNr-(SourceY*SourceNx);
                               int TargetY = TargetNr/TargetNx;
                               int TargetX = TargetNr-(TargetY*TargetNx);
                               matrix.SetValue(SynWeights[(*Index)], SourceX, SourceY, TargetX, TargetY);
                      }
             }
    }
}

float VecConnection::GetMaxWeight()
{
    float MaxWeight=0;
    for (int SourceNr=0;SourceNr<ns;++SourceNr) {
        for(int CurDelay=0;CurDelay<maximumDelay;++CurDelay)
        {
            for (vector <T_NSynapses>::iterator it=delays[SourceNr][CurDelay].begin();
                 it != delays[SourceNr][CurDelay].end();++it) 
            {           
                int CurTarget=SynTargetNr[(*it)];
                int CurSource=SynSourceNr[(*it)];
                if (SourceNr!=CurSource) {
                    cerr << "*******ERROR*****************************************\n";
                    cerr << "float GetSourceWeights(int SourceNr, vector<Synapse>& SynList)\n";
                    cerr << "korrupt\n";
                }
                float CurWeight=SynWeights[(*it)];
                if (CurWeight>MaxWeight) {
                    MaxWeight=CurWeight;
                }
            }
        }
    }
    return MaxWeight;   
}


float VecConnection::GetMinWeight()
{
    bool MinNotInitialized=true;
    float MinWeight;
    for (int SourceNr=0;SourceNr<ns;++SourceNr) {
        for(int CurDelay=0;CurDelay<maximumDelay;++CurDelay)
        {
            for (vector <T_NSynapses>::iterator it=delays[SourceNr][CurDelay].begin();
                 it != delays[SourceNr][CurDelay].end();++it) 
            {           
                int CurTarget=SynTargetNr[(*it)];
                int CurSource=SynSourceNr[(*it)];
                if (SourceNr!=CurSource) {
                    cerr << "*******ERROR*****************************************\n";
                    cerr << "float GetSourceWeights(int SourceNr, vector<Synapse>& SynList)\n";
                    cerr << "korrupt\n";
                }
                float CurWeight=SynWeights[(*it)];
                if (MinNotInitialized)   MinWeight = CurWeight;
                if (CurWeight<MinWeight) MinWeight = CurWeight;
              
            }
        }
    }
    return MinWeight;   
}


int VecConnection::GetMaxDelay()
{
    int MaxDelay=0;
    for (int SourceNr=0;SourceNr<ns;++SourceNr) {
        for(int CurDelay=MaxDelay;CurDelay<maximumDelay;++CurDelay) {
            if (delays[SourceNr][CurDelay].size()>0) {
                if (CurDelay>MaxDelay) {
                    MaxDelay=CurDelay;
                }
            }
        }
    }
    return MaxDelay;    
}

long VecConnection::calcMemoryConsumption()
{
    long long MemSum=0;
    MemSum += SynWeights.size()*sizeof(float);
    MemSum += SynTargetNr.size()*sizeof(T_NNeurons);
    MemSum += SynSourceNr.size()*sizeof(T_NNeurons);
    MemSum += SynDelays.size()*sizeof(T_Delays);
    
    int PreSynSum = 0;
    PreSynSum += PreSynNr.size()*sizeof(vector<T_NSynapses>);    
    for (int i=0;i<PreSynNr.size();++i) {
        PreSynSum += PreSynNr[i].size()*sizeof(T_NSynapses);
    }
    MemSum += PreSynSum;
    cout << "PreSynSum=" << PreSynSum << "\n";
    cout << "NSynapses = " << NSynapses << " MAXINT=" << INT_MAX << "\n";
    
    int DelaySum=0;
    DelaySum += delays.size()*sizeof(vector<vector<T_NSynapses> >);
    for (int i=0;i<delays.size();++i) {
        DelaySum += delays[i].size()*sizeof(vector<T_NSynapses>);
        for (int j=0;j<delays[i].size();++j) {
            DelaySum += delays[i][j].size()*sizeof(T_NSynapses);
        }
    }
    MemSum += DelaySum;
    cout << "DelaySum=" << DelaySum << "\n";
    return MemSum;
}

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

Connection* LoadConnection(const char* FileName)
{
    // try loading VecConnection
    VecConnection* vcon = new VecConnection();
    int err=vcon->Load(FileName, "");
    if (err==0) {
        return vcon;
    } else  {
        delete vcon;
        vcon=0;
        Dout(dc::con, "Could not load VecConnection file " << FileName);
        Dout(dc::con, "try loading as connection");
        connection* con = new connection();
        err = con->Load(FileName, "");
        if (err==0) {
            return con;
        } else {
            delete con;
            cerr << "ERROR: failed loading connection file " << FileName << "\n";
        }
    }
    return 0;
}


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

VecDepressiveConnection::VecDepressiveConnection(
     layer* SL, layer* TL, 
     csimInputChannel InputNumber, float _TauRec, float _U_se, float InitializeFrequency)
     : VecConnection(SL, TL, InputNumber), 
       U_SE(_U_se),
       TauRec(_TauRec/dt), U_se_fac(1-_U_se)

{     
    E0=InitializeEfficacy(InitializeFrequency);
     Dout(dc::con, "  VecDepressiveConnection::VecDepressiveConnection ");
    Dout(dc::con, "  Initialize Efficacy values to E0=" << E0);
}

int VecDepressiveConnection::proceede(int TotalTime)
{
     int t = TotalTime % MacroTimeStep;
     int i,j,k,mi,ipre;
     // calculate input for target layer  
     k=SourceLayer->N_firings;
     // Nur Spikes, die nicht laenger als maximumDelay in der 
     // Vergangeheit liegen werden beruecksichtigt. firins-array
     // besteht aus firings[SpikeNr][0]:zeitpunkte, 
     // firnigs[SpikeNr][1]:Neuronennummer
     while (t-SourceLayer->firings[--k][0] < maximumDelay)
     {
      int CurDelay = t-SourceLayer->firings[k][0];
      ipre = SourceLayer->firings[k][1];
         for (vector <T_NSynapses>::iterator it=delays[ipre][CurDelay].begin();
           it != delays[ipre][CurDelay].end(); 
           ++it)
      {
           // exponential recovery from last epsp
          if (TauRec >0) {
              
            Efficacy[(*it)] =  (1-(1-Efficacy[(*it)])*exp(-(t-LastEpsp[(*it)])/TauRec)); 
          } else {
              Efficacy[(*it)]=1;
          }
           InputPointer[SynTargetNr[(*it)]] += SynWeights[(*it)]*Efficacy[(*it)];
           Efficacy[(*it)] *= U_se_fac;
           LastEpsp[(*it)] = t;
      }
     }     
     
     if (learn == true) {
      learnobj->proceede(TotalTime);
     }
}

int VecDepressiveConnection::PushBackNewSynapse(int source, int target, float weight, int delay) 
{
     SynSourceNr.push_back(source);
     SynTargetNr.push_back(target);
     SynWeights.push_back(weight);
     SynDelays.push_back(delay);
     Efficacy.push_back(E0);
     LastEpsp.push_back(0);
     PreSynNr[target].push_back(NSynapses);     
     ++NSynapses;
}

int VecDepressiveConnection::ReserveSynapses(int _nsynapses) 
{
     SynWeights.reserve(_nsynapses);
     SynTargetNr.reserve(_nsynapses);
     SynSourceNr.reserve(_nsynapses);
     SynDelays.reserve(_nsynapses);
     Efficacy.reserve(_nsynapses);
     LastEpsp.reserve(_nsynapses);
}


int VecDepressiveConnection::reset(int t)
{
     int nsyn=Efficacy.size();
     if (nsyn!=LastEpsp.size()) {
      Dout(dc::con, "ERROR in VecDepressiveConnection::reset"); exit(1);
     }
     for (int i=0;i<nsyn;++i) {
      Efficacy[i]=E0;
      LastEpsp[i]=1;
     }
}

void VecDepressiveConnection::SetParameter(ParaType p, double value)
{   
    switch (p)
    {
    case PARA_U_SE:
        U_SE=value;
        U_se_fac =1-U_SE;
        break;
    case PARA_TAU_REC:
        TauRec = value/dt;
        break;  
    }
}

int VecDepressiveConnection::prepare(int step)
{    
     VecConnection::prepare(step);
     for (vector <int>::iterator it=LastEpsp.begin(); it!=LastEpsp.end();++it) {
      (*it) -=MacroTimeStep;
     }
}

int VecDepressiveConnection::WriteSimInfo(fstream &fw)
{
     stringstream sstr;
     sstr << "<TauRec value=\"" << TauRec*dt << "\"/> \n";
     sstr << "<U_SE value=\"" << U_SE << "\"/> \n";
     VecConnection::WriteSimInfo(fw, sstr.str());
}

int VecDepressiveConnection::DeleteSynapseArrays()
{
    LastEpsp.clear();
    Efficacy.clear();
    VecConnection::DeleteSynapseArrays();
}


/** initialize Efficacy according to assumed spike frequency
 \f[
    E_0=\frac{1-e^{\frac{-1}{f\tau}}-U_{SE}}{1-e^{\frac{-1}{f\tau}}}
\f]
@param f: frequency 
*/  
float VecDepressiveConnection::InitializeEfficacy(float f)
{
    float E0=1.0;
    float EinsMinExp=0.0;
    float EFact = exp(-1/(f*0.001*TauRec));
    if (f>0.) {
        E0=(U_se_fac-U_se_fac*EFact)/(1-U_se_fac*EFact);
        Dout(dc::con, "E0=" << E0 << "");
    }
    if (E0>0) {
        return E0;
    } else {
        return 0.0;
    }
}

long VecDepressiveConnection::calcMemoryConsumption()
{
    long long MemSum=VecConnection::calcMemoryConsumption();
    MemSum += Efficacy.size()*sizeof(float);
    MemSum += LastEpsp.size()*sizeof(int);
    return MemSum;
}

bool VecConnection::hasEqualSynapses(const VecConnection & OtherVecCon)
{
    bool WeightsAreEqual = vectorsAreEqual(SynWeights, OtherVecCon.SynWeights);
    bool TargetNrAreEqual = vectorsAreEqual(SynTargetNr, OtherVecCon.SynTargetNr);
    bool SourceNrAreEqual = vectorsAreEqual(SynSourceNr, OtherVecCon.SynSourceNr);
    bool DelaysAreEqual = vectorsAreEqual(SynDelays, OtherVecCon.SynDelays);

    Dout(dc::con, "VecConnection::hasEqualSynapses: NSyn=" << NSynapses
         << "  OtherNSyn=" << OtherVecCon.NSynapses
         << "  W=" << WeightsAreEqual
         << "  T=" << TargetNrAreEqual
         << "  S=" << SourceNrAreEqual
         << "  D=" << DelaysAreEqual);
    //Dout(dc::con, stringifyVector(SynDelays));
    //Dout(dc::con, stringifyVector(OtherVecCon.SynDelays));
    return (NSynapses == OtherVecCon.NSynapses)
        && WeightsAreEqual && TargetNrAreEqual && SourceNrAreEqual && DelaysAreEqual ;
}

void VecConnection::updateNSynapses()
{
    NSynapses=0;
    for (vector<T_Delays>::iterator it=SynDelays.begin(); it!=SynDelays.end(); ++it) {
        if ((*it) != DELETED_SYNAPSE) {
            ++NSynapses;
        }
    }
}