nnrbm.cpp

#include "nnrbm.h"

extern StreamOutput cout;

NNRBM::NNRBM()
{
    // init member vars
    m_nrTargets = 0;
    m_nrInputs = 0;
    m_nrExamplesTrain = 0;
    m_nrExamplesProbe = 0;
    m_inputsTrain = 0;
    m_inputsProbe = 0;
    m_targetsTrain = 0;
    m_targetsProbe = 0;
    m_initWeightFactor = 0;
    m_globalEpochs = 0;
    m_RPROP_etaPos = 0;
    m_RPROP_etaNeg = 0;
    m_RPROP_updateMin = 0;
    m_RPROP_updateMax = 0;
    m_learnRate = 0;
    m_learnRateMin = 0;
    m_learnrateDecreaseRate = 0;
    m_learnrateDecreaseRateEpoch = 0;
    m_momentum = 0;
    m_weightDecay = 0;
    m_minUpdateBound = 0;
    m_batchSize = 0;
    m_scaleOutputs = 0;
    m_offsetOutputs = 0;
    m_maxEpochs = 0;
    m_useBLAS = 0;
    m_enableRPROP = 0;
    m_normalTrainStopping = 0;
    m_nrLayer = 0;
    m_neuronsPerLayer = 0;
    m_nrWeights = 0;
    m_nrOutputs = 0;
    m_nrLayWeights = 0;
    m_outputs = 0;
    m_outputsTmp = 0;
    m_derivates = 0;
    m_d1 = 0;
    m_weights = 0;
    m_weightsTmp0 = 0;
    m_weightsTmp1 = 0;
    m_weightsTmp2 = 0;
    m_weightsBatchUpdate = 0;
    m_weightsOld = 0;
    m_weightsOldOld = 0;
    m_deltaW = 0;
    m_deltaWOld = 0;
    m_adaptiveRPROPlRate = 0;
    m_enableL1Regularization = 0;
    m_errorFunctionMAE = 0;
    m_activationFunctionPerLayer = 0;
    m_enableAutoencoder = 0;
    m_nrLayWeightOffsets = 0;
    m_sumSquaredError = 0.0;
    m_sumSquaredErrorSamples = 0.0;
    m_rbmLearnrateWeights = 0.0;
    m_rbmLearnrateBiasVis = 0.0;
    m_rbmLearnrateBiasHid = 0.0;
    m_rbmWeightDecay = 0.0;
    m_rbmMaxEpochs = 0;
}

NNRBM::~NNRBM()
{
    if ( m_neuronsPerLayer )
        delete[] m_neuronsPerLayer;
    m_neuronsPerLayer = 0;
    if ( m_nrLayWeights )
        delete[] m_nrLayWeights;
    m_nrLayWeights = 0;
    if ( m_outputs )
        delete[] m_outputs;
    m_outputs = 0;
    if ( m_outputsTmp )
        delete[] m_outputsTmp;
    m_outputsTmp = 0;
    if ( m_derivates )
        delete[] m_derivates;
    m_derivates = 0;
    if ( m_d1 )
        delete[] m_d1;
    m_d1 = 0;
    if ( m_weights )
        delete[] m_weights;
    m_weights = 0;
    if ( m_weightsTmp0 )
        delete[] m_weightsTmp0;
    m_weightsTmp0 = 0;
    if ( m_weightsTmp1 )
        delete[] m_weightsTmp1;
    m_weightsTmp1 = 0;
    if ( m_weightsTmp2 )
        delete[] m_weightsTmp2;
    m_weightsTmp2 = 0;
    if ( m_weightsBatchUpdate )
        delete[] m_weightsBatchUpdate;
    m_weightsBatchUpdate = 0;
    if ( m_weightsOld )
        delete[] m_weightsOld;
    m_weightsOld = 0;
    if ( m_weightsOldOld )
        delete[] m_weightsOldOld;
    m_weightsOldOld = 0;
    if ( m_deltaW )
        delete[] m_deltaW;
    m_deltaW = 0;
    if ( m_deltaWOld )
        delete[] m_deltaWOld;
    m_deltaWOld = 0;
    if ( m_adaptiveRPROPlRate )
        delete[] m_adaptiveRPROPlRate;
    m_adaptiveRPROPlRate = 0;
    if ( m_activationFunctionPerLayer )
        delete[] m_activationFunctionPerLayer;
    m_activationFunctionPerLayer = 0;
    if ( m_nrLayWeightOffsets )
        delete[] m_nrLayWeightOffsets;
    m_nrLayWeightOffsets = 0;
}

void NNRBM::enableErrorFunctionMAE ( bool en )
{
    m_errorFunctionMAE = en;
    cout<<"errorFunctionMAE:"<<m_errorFunctionMAE<<endl;
}

void NNRBM::setNrTargets ( int n )
{
    m_nrTargets = n;
    cout<<"nrTargets: "<<m_nrTargets<<endl;
}


void NNRBM::setNrInputs ( int n )
{
    m_nrInputs = n;
    cout<<"nrInputs: "<<m_nrInputs<<endl;
}
void NNRBM::setNrExamplesTrain ( int n )
{
    m_nrExamplesTrain = n;
    //cout<<"nrExamplesTrain: "<<m_nrExamplesTrain<<endl;
}
void NNRBM::setNrExamplesProbe ( int n )
{
    m_nrExamplesProbe = n;
    cout<<"nrExamplesProbe: "<<m_nrExamplesProbe<<endl;
}
void NNRBM::setTrainInputs ( REAL* inputs )
{
    m_inputsTrain = inputs;
    //cout<<"inputsTrain: "<<m_inputsTrain<<endl;
}

void NNRBM::setTrainTargets ( REAL* targets )
{
    m_targetsTrain = targets;
    //cout<<"targetsTrain: "<<m_targetsTrain<<endl;
}

void NNRBM::setProbeInputs ( REAL* inputs )
{
    m_inputsProbe = inputs;
    //cout<<"inputsProbe: "<<m_inputsProbe<<endl;
}

void NNRBM::setProbeTargets ( REAL* targets )
{
    m_targetsProbe = targets;
    //cout<<"targetsProbe: "<<m_targetsProbe<<endl;
}

void NNRBM::setInitWeightFactor ( REAL factor )
{
    m_initWeightFactor = factor;
    cout<<"initWeightFactor: "<<m_initWeightFactor<<endl;
}

void NNRBM::setLearnrate ( REAL learnrate )
{
    m_learnRate = learnrate;
    cout<<"learnRate: "<<m_learnRate<<endl;
}

void NNRBM::setLearnrateMinimum ( REAL learnrateMin )
{
    m_learnRateMin = learnrateMin;
    cout<<"learnRateMin: "<<m_learnRateMin<<endl;
}

void NNRBM::setLearnrateSubtractionValueAfterEverySample ( REAL learnrateDecreaseRate )
{
    m_learnrateDecreaseRate = learnrateDecreaseRate;
    cout<<"learnrateDecreaseRate: "<<m_learnrateDecreaseRate<<endl;
}


void NNRBM::setLearnrateSubtractionValueAfterEveryEpoch ( REAL learnrateDecreaseRate )
{
    m_learnrateDecreaseRateEpoch = learnrateDecreaseRate;
    cout<<"learnrateDecreaseRateEpoch: "<<m_learnrateDecreaseRateEpoch<<endl;
}

void NNRBM::setMomentum ( REAL momentum )
{
    m_momentum = momentum;
    cout<<"momentum: "<<m_momentum<<endl;
}

void NNRBM::setWeightDecay ( REAL weightDecay )
{
    m_weightDecay = weightDecay;
    cout<<"weightDecay: "<<m_weightDecay<<endl;
}

void NNRBM::setBatchSize ( int size )
{
    m_batchSize = size;
    cout<<"batchSize: "<<m_batchSize<<endl;
}

void NNRBM::setMinUpdateErrorBound ( REAL minUpdateBound )
{
    m_minUpdateBound = minUpdateBound;
    cout<<"minUpdateBound: "<<m_minUpdateBound<<endl;
}

void NNRBM::setMaxEpochs ( int epochs )
{
    m_maxEpochs = epochs;
    cout<<"maxEpochs: "<<m_maxEpochs<<endl;
}

void NNRBM::setRPROPPosNeg ( REAL etaPos, REAL etaNeg )
{
    m_RPROP_etaPos = etaPos;
    m_RPROP_etaNeg = etaNeg;
    cout<<"RPROP_etaPos: "<<m_RPROP_etaPos<<"  RPROP_etaNeg: "<<m_RPROP_etaNeg<<endl;
}

void NNRBM::setRPROPMinMaxUpdate ( REAL min, REAL max )
{
    m_RPROP_updateMin = min;
    m_RPROP_updateMax = max;
    cout<<"RPROP_updateMin: "<<m_RPROP_updateMin<<"  RPROP_updateMax: "<<m_RPROP_updateMax<<endl;
}

void NNRBM::setScaleOffset ( REAL scale, REAL offset )
{
    m_scaleOutputs = scale;
    m_offsetOutputs = offset;
    cout<<"scaleOutputs: "<<m_scaleOutputs<<"   offsetOutputs: "<<m_offsetOutputs<<"  [transformation: output = outputNN * scale + offset]"<<endl;
}

void NNRBM::setNormalTrainStopping ( bool en )
{
    m_normalTrainStopping = en;
    cout<<"normalTrainStopping: "<<m_normalTrainStopping<<endl;
}

void NNRBM::setL1Regularization ( bool en )
{
    m_enableL1Regularization = en;
    cout<<"enableL1Regularization: "<<m_enableL1Regularization<<endl;
}

void NNRBM::enableRPROP ( bool en )
{
    m_enableRPROP = en;
    cout<<"enableRPROP: "<<m_enableRPROP<<endl;
}

void NNRBM::useBLASforTraining ( bool enable )
{
    m_useBLAS = enable;
    cout<<"useBLAS: "<<m_useBLAS<<endl;
}

void NNRBM::setNNStructure ( int nrLayer, int* neuronsPerLayer, bool lastLinearLayer, int* layerType )
{
    m_nrLayer = nrLayer;
    cout<<"nrLayer: "<<m_nrLayer<<endl;

    cout<<"#layers: "<<m_nrLayer<<" ("<< ( m_nrLayer-1 ) <<" hidden layer, 1 output layer)"<<endl;

    // alloc space for structure variables
    m_neuronsPerLayer = new int[m_nrLayer+1];
    m_neuronsPerLayer[0] = m_nrInputs;  // number of inputs
    for ( int i=0;i<m_nrLayer-1;i++ )
        m_neuronsPerLayer[1+i] = neuronsPerLayer[i];
    m_neuronsPerLayer[m_nrLayer] = m_nrTargets;  // one output

    cout<<"Neurons    per Layer: ";
    for ( int i=0;i<m_nrLayer+1;i++ )
        cout<<m_neuronsPerLayer[i]<<" ";
    cout<<endl;

    cout<<"Outputs    per Layer: ";
    for ( int i=0;i<m_nrLayer+1;i++ )
        cout<<m_neuronsPerLayer[i]+1<<" ";
    cout<<endl;

    cout<<"OutOffsets per Layer: ";
    int cnt=0;
    for ( int i=0;i<m_nrLayer+1;i++ )
    {
        cout<<cnt<<" ";
        cnt += m_neuronsPerLayer[i]+1;
    }
    cout<<endl;

    cout<<"Act. function per Layer(0=sig, 1=lin, 2=tanh, 3=softmax): ";
    m_activationFunctionPerLayer = new int[m_nrLayer+1];
    for ( int i=0;i<m_nrLayer+1;i++ )
        m_activationFunctionPerLayer[i] = 0;
    if ( layerType )
    {
        for ( int i=0;i<m_nrLayer+1;i++ )
            m_activationFunctionPerLayer[i] = layerType[i];
    }
    if ( m_enableAutoencoder )
        m_activationFunctionPerLayer[m_nrLayer/2] = 1;
    if ( lastLinearLayer )
        m_activationFunctionPerLayer[m_nrLayer] = 1;
    for ( int i=0;i<m_nrLayer+1;i++ )
        cout<<m_activationFunctionPerLayer[i]<<" ";
    cout<<endl;

    // init the total number of weights and outputs
    m_nrWeights = 0;
    m_nrOutputs = m_neuronsPerLayer[0] + 1;
    m_nrLayWeights = new int[m_nrLayer+1];
    m_nrLayWeightOffsets = new int[m_nrLayer+2];
    m_nrLayWeights[0] = 0;
    for ( int i=0;i<m_nrLayer;i++ )
    {
        m_nrLayWeights[i+1] = m_neuronsPerLayer[i+1] * ( m_neuronsPerLayer[i]+1 );  // +1 for input bias
        m_nrWeights += m_nrLayWeights[i+1];
        m_nrOutputs += m_neuronsPerLayer[i+1] + 1;  // +1 for input bias
    }

    // print it
    cout<<"Weights       per Layer: ";
    for ( int i=0;i<m_nrLayer+1;i++ )
        cout<<m_nrLayWeights[i]<<" ";
    cout<<endl;

    cout<<"WeightOffsets per Layer: ";
    m_nrLayWeightOffsets[0] = 0;
    for ( int i=0;i<m_nrLayer+1;i++ )
    {
        cout<<m_nrLayWeightOffsets[i]<<" ";
        m_nrLayWeightOffsets[i+1] = m_nrLayWeightOffsets[i] + m_nrLayWeights[i];
    }
    cout<<endl;

    cout<<"nrOutputs="<<m_nrOutputs<<"  nrWeights="<<m_nrWeights<<endl;

    // allocate the inner calculation structure
    m_outputs = new REAL[m_nrOutputs];
    m_outputsTmp = new REAL[m_nrTargets];
    m_derivates = new REAL[m_nrOutputs];
    m_d1 = new REAL[m_nrOutputs];

    for ( int i=0;i<m_nrOutputs;i++ ) // init as biases
    {
        m_outputs[i] = 1.0;
        m_derivates[i] = 0.0;
        m_d1[i] = 0.0;
    }

    // allocate weights and temp vars
    m_weights = new REAL[m_nrWeights];
    m_weightsTmp0 = new REAL[m_nrWeights];
    m_weightsTmp1 = new REAL[m_nrWeights];
    m_weightsTmp2 = new REAL[m_nrWeights];
    m_weightsBatchUpdate = new REAL[m_nrWeights];
    m_weightsOld = new REAL[m_nrWeights];
    m_weightsOldOld = new REAL[m_nrWeights];
    m_deltaW = new REAL[m_nrWeights];

    m_deltaWOld = new REAL[m_nrWeights];
    m_adaptiveRPROPlRate = new REAL[m_nrWeights];
    for ( int i=0;i<m_nrWeights;i++ )
    {
        m_deltaWOld[i] = 0.0;
        m_adaptiveRPROPlRate[i] = m_learnRate;
    }
    for ( int i=0;i<m_nrWeights;i++ )
        m_weights[i] = m_weightsOld[i] = m_deltaW[i] = m_weightsTmp0[i] = m_weightsTmp1[i] = m_weightsTmp2[i] = 0.0;

    // this should be implemented (LeCun suggest such a linear factor in the activation function)
    //m_linFac = 0.01;
    //cout<<"linFac="<<m_linFac<<" (no active, just tanh used)"<<endl;
}

REAL NNRBM::getInitWeight ( int fanIn )
{
    double nr = 2.0* ( rand() / ( double ) RAND_MAX-0.5 );  // -1 .. +1
    return ( 1.0/sqrt ( ( double ) fanIn ) ) * nr;
}

void NNRBM::initNNWeights ( time_t seed )
{
    srand ( seed );
    cout<<"init weights "<<endl;
    REAL factor = m_initWeightFactor;
    int cnt = 0;
    for ( int i=0;i<m_nrLayer;i++ ) // through all layers
    {
        int n = m_neuronsPerLayer[i+1];
        int nprev = m_neuronsPerLayer[i] + 1;  // +1 for bias
        for ( int j=0;j<n;j++ ) // all neurons per layer
        {
            for ( int k=0;k<nprev;k++ ) // all weights from this neuron
            {
                m_weights[cnt] = m_weightsOld[i] = m_weightsOldOld[i] = getInitWeight ( nprev ) * factor;
                cnt++;
            }
        }
    }

    // check the number
    if ( cnt != m_nrWeights )
        assert ( false );
}

void NNRBM::setRBMLearnParams ( REAL learnrateWeights, REAL learnrateBiasVis, REAL learnrateBiasHid, REAL weightDecay, int maxEpoch )
{
    m_rbmLearnrateWeights = learnrateWeights;
    m_rbmLearnrateBiasVis = learnrateBiasVis;
    m_rbmLearnrateBiasHid = learnrateBiasHid;
    m_rbmWeightDecay = weightDecay;
    m_rbmMaxEpochs = maxEpoch;
}

int NNRBM::getWeightIndex ( int layer, int neuron, int weight )
{
    if ( layer == 0 )
        assert ( false );

    int nrNeur = m_neuronsPerLayer[layer];
    int nrNeurPrev = m_neuronsPerLayer[layer-1];
    if ( neuron >= nrNeur )
    {
        cout<<"neuron:"<<neuron<<" nrNeur:"<<nrNeur<<endl;
        assert ( false );
    }
    if ( weight >= nrNeurPrev )
    {
        cout<<"weight:"<<weight<<" nrNeurPrev:"<<nrNeurPrev<<endl;
        assert ( false );
    }

    int ind = m_nrLayWeightOffsets[layer];
    if ( layer == 1 ) // input layer
        ind += 1 + weight + neuron* ( nrNeurPrev + 1 );
    else
        ind += weight + neuron* ( nrNeurPrev + 1 );

    if ( ind >= m_nrWeights )
    {
        cout<<"ind:"<<ind<<" m_nrWeights:"<<m_nrWeights<<endl;
        assert ( false );
    }

    return ind;
}

int NNRBM::getBiasIndex ( int layer, int neuron )
{
    if ( layer == 0 )
        assert ( false );

    int nrNeur = m_neuronsPerLayer[layer];
    int nrNeurPrev = m_neuronsPerLayer[layer-1];
    if ( neuron >= nrNeur )
    {
        cout<<"neuron:"<<neuron<<" nrNeur:"<<nrNeur<<endl;
        assert ( false );
    }
    int ind = m_nrLayWeightOffsets[layer];
    if ( layer == 1 ) // input layer
        ind += neuron* ( nrNeurPrev + 1 );
    else
        ind += nrNeurPrev + neuron* ( nrNeurPrev + 1 );

    if ( ind >= m_nrWeights )
    {
        cout<<"ind:"<<ind<<" m_nrWeights:"<<m_nrWeights<<endl;
        assert ( false );
    }

    return ind;
}

int NNRBM::getOutputIndex ( int layer, int neuron )
{
    if ( layer == 0 || layer > m_nrLayer )
        assert ( false );

    if ( neuron >= m_neuronsPerLayer[layer] )
        assert ( false );

    int ind = 0;
    for ( int i=0;i<layer;i++ )
        ind += m_neuronsPerLayer[i] + 1;

    return ind + neuron;
}

vector<int> NNRBM::getEncoder()
{
    vector<int> v;
    cout<<"Encoder:"<<endl;
    int midLayer = m_nrLayer/2 + 1;
    int weightCnt = 0;
    for ( int i=0;i<midLayer;i++ )
    {
        v.push_back ( m_neuronsPerLayer[i] );
        weightCnt += m_nrLayWeights[i];
        cout<<"neurons|weights: "<<m_neuronsPerLayer[i]<<"|"<<m_nrLayWeights[i]<<endl;
    }
    v.push_back ( weightCnt );
    return v;
}

void NNRBM::rbmPretraining ( REAL* input, REAL* target, int nSamples, int nTarget, int firstNLayer, int crossRun, int nFinetuningEpochs, double finetuningLearnRate )
{
    cout<<endl<<endl<<"=== set structure for fine tuning of the unrolled autoencoder ==="<<endl;
    m_nnAuto = new NNRBM();
    m_nnAuto->setInitWeightFactor ( 0.0 );
    m_nnAuto->setScaleOffset ( m_scaleOutputs, m_offsetOutputs );
    m_nnAuto->setNrTargets ( m_nrInputs );
    m_nnAuto->setNrInputs ( m_nrInputs );
    m_nnAuto->setNrExamplesTrain ( nSamples );
    m_nnAuto->setTrainInputs ( input );
    m_nnAuto->setTrainTargets ( input );

    // learn parameters
    m_nnAuto->setInitWeightFactor ( m_initWeightFactor );
    if ( finetuningLearnRate != 0.0 )
        m_nnAuto->setLearnrate ( finetuningLearnRate );
    else
        m_nnAuto->setLearnrate ( m_learnRate );
    m_nnAuto->setLearnrateMinimum ( m_learnRateMin );
    m_nnAuto->setLearnrateSubtractionValueAfterEverySample ( m_learnrateDecreaseRate );
    m_nnAuto->setLearnrateSubtractionValueAfterEveryEpoch ( m_learnrateDecreaseRateEpoch );
    m_nnAuto->setMomentum ( m_momentum );
    m_nnAuto->setWeightDecay ( 0.0 );
    m_nnAuto->setMinUpdateErrorBound ( m_minUpdateBound );
    m_nnAuto->setBatchSize ( m_batchSize );
    m_nnAuto->setMaxEpochs ( m_maxEpochs );
    m_nnAuto->setL1Regularization ( m_enableL1Regularization );
    m_nnAuto->enableErrorFunctionMAE ( m_errorFunctionMAE );
    m_nnAuto->setNormalTrainStopping ( true );
    m_nnAuto->useBLASforTraining ( m_useBLAS );


    int nLayer = m_nrLayer;
    if ( firstNLayer != -1 )
        nLayer = firstNLayer;

    int* layerConfig = new int[2* ( nLayer-1 ) ];
    for ( int i=0;i<2* ( nLayer-1 ) - 1;i++ )
    {
        if ( i < nLayer-1 )
            layerConfig[i] = m_neuronsPerLayer[i+1];
        else
            layerConfig[i] = m_neuronsPerLayer[2* ( nLayer-1 ) - 1 - i];
    }
    m_nnAuto->m_enableAutoencoder = true;
    m_nnAuto->setNNStructure ( 2* ( nLayer-1 ), layerConfig );
    m_nnAuto->initNNWeights ( 0 );


    cout<<endl<<"=== RBM pretraining ("<<nSamples<<" samples) ==="<<endl;
    int weightOffset = 0;
    REAL* batchdata = new REAL[nSamples*m_neuronsPerLayer[0]];
    for ( int j=0;j<nSamples*m_neuronsPerLayer[0];j++ )
        batchdata[j] = input[j];

    // all layers from output to input
    for ( int i=0;i<nLayer-1;i++ )
    {
        int nVis = m_neuronsPerLayer[i];
        int nHid = m_neuronsPerLayer[i+1];
        cout<<endl<<"layer:"<<i<<" nVis:"<<nVis<<" nHid:"<<nHid<<endl;

        // probabilities, biases and other
        REAL* batchdataNext = new REAL[nSamples*nHid];
        REAL* vishid = new REAL[nVis*nHid];
        REAL* visbias = new REAL[nVis];
        REAL* hidbias = new REAL[nHid];
        REAL* poshidprobs = new REAL[nHid];
        REAL* neghidprobs = new REAL[nHid];
        REAL* hidbiasinc = new REAL[nHid];
        REAL* visbiasinc = new REAL[nVis];
        REAL* poshidstates = new REAL[nHid];
        REAL* negdata = new REAL[nVis];
        REAL* neghidact = new REAL[nHid];
        REAL* negvisact = new REAL[nVis];
        REAL* poshidact = new REAL[nHid];
        REAL* posvisact = new REAL[nVis];
        REAL* probs = new REAL[nHid];
        REAL constant = 0.0;
        //REAL constant = -1.0;
        for ( int j=0;j<nVis;j++ )
        {
            visbias[j] = constant;
            visbiasinc[j] = constant;
            negdata[j] = constant;
        }
        for ( int j=0;j<nHid;j++ )
        {
            hidbias[j] = constant;
            hidbiasinc[j] = constant;
            poshidprobs[j] = constant;
            neghidprobs[j] = constant;
            poshidstates[j] = constant;
        }
        for ( int j=0;j<nVis*nHid;j++ )
        {
            vishid[j] = constant;
        }

        // init weights
        unsigned int seed = rand();
        cout<<"seed:"<<seed<<endl;
        //GAUSS_DISTRIBUTION(vishid, nVis*nHid, 0.0, 1.0/(double)sqrt((double)nVis), &seed);
        GAUSS_DISTRIBUTION ( vishid, nVis*nHid, 0.0, 0.1, &seed );
        time_t t0 = time ( 0 );


        // layer-wise training
        for ( int epoch=0;epoch<m_rbmMaxEpochs;epoch++ )
        {
            double errsum = 0.0;
            int progress = nSamples/10 + 1, errcnt = 0;

            for ( int sample=0;sample<nSamples;sample++ )
            {
                // data vector
                REAL* data = batchdata + nVis * sample;
                REAL min = 1e10, max = -1e10, mean = 0.0;
                for ( int k=0;k<nVis;k++ )
                {
                    if ( min > data[k] )
                        min = data[k];
                    if ( max < data[k] )
                        max = data[k];
                    mean += data[k];
                }
                mean /= ( double ) nVis;

                if ( min > 1.0 || min < 0.0 )
                {
                    cout<<"min:"<<min<<endl;
                    assert ( false );
                }
                if ( max > 1.0 || max < 0.0 )
                {
                    cout<<"max:"<<max<<endl;
                    assert ( false );
                }


                // ===== POSITIVE PHASE =====
                // this is like forward calculation
                // after calc the weighted sum per neuron, the values of the neurons get sampled

                // poshidprobs = 1./(1 + exp(-data*vishid - repmat(hidbiases,numcases,1)));

                // calc: poshidprobs = vishid * data

                //for(int k=0;k<nHid;k++)
                //{
                //    REAL sum = 0.0;
                //    for(int l=0;l<nVis;l++)
                //        sum += data[l] * vishid[l + k*nVis];
                //    poshidprobs[k] = 1.0 / (1.0 + exp(-(sum+hidbias[k])));
                //}

                CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, nHid, nVis, 1.0, vishid, nVis, data, 1, 0.0, poshidprobs, 1 );

                //for(int k=0;k<nHid;k++)
                //    poshidprobs[k] = 1.0 / (1.0 + exp(-(poshidprobs[k]+hidbias[k])));
                V_ADD ( nHid, poshidprobs, hidbias, poshidprobs );  // poshidprobs[k]+hidbias[k]
                V_MULCI ( -1.0, poshidprobs, nHid );             // *(-1)
                V_EXP ( poshidprobs, nHid );                     // exp(x)
                V_ADDCI ( 1.0, poshidprobs, nHid );              // +1
                V_INV ( nHid, poshidprobs, poshidprobs );        // inv(x)

                // batchposhidprobs(:,:,batch)=poshidprobs;

                // posprods    = data' * poshidprobs;
                //for(int k=0;k<nHid;k++)
                //    for(int l=0;l<nVis;l++)
                //        posprods[l+k*nVis] = data[l] * poshidprobs[k];

                // poshidact   = sum(poshidprobs);

                //for(int k=0;k<nHid;k++)
                //    poshidact[k] = poshidprobs[k];
                V_COPY ( poshidprobs, poshidact, nHid );   // poshidact = poshidprobs

                // posvisact = sum(data);
                //for(int k=0;k<nVis;k++)
                //    posvisact[k] = data[k];
                V_COPY ( data, posvisact, nVis );   // posvisact = data

                // poshidstates = poshidprobs > rand(numcases,numhid);
                for ( int k=0;k<nHid;k++ )
                    probs[k] = ( double ) rand() / ( double ) RAND_MAX;
                //UNIFORM_DISTRIBUTION(probs, nHid, 0.0, 1.0, seed);

                for ( int k=0;k<nHid;k++ )
                {
                    if ( poshidprobs[k] > probs[k] )
                        poshidstates[k] = 1.0;
                    else
                        poshidstates[k] = 0.0;
                }


                // ===== NEGATIVE PHASE =====
                // reconstruct the input sample, based on binary representation in the hidden layer
                // use the reconstruction error as signal for learning

                //negdata = 1./(1 + exp(-poshidstates*vishid' - repmat(visbiases,numcases,1)));
                //for(int k=0;k<nVis;k++)
                //{
                //    REAL sum = 0.0;
                //    for(int l=0;l<nHid;l++)
                //        sum += poshidstates[l] * vishid[k + l*nVis];
                //    negdata[k] = sum;//1.0 / (1.0 + exp(-(sum+visbias[k])));
                //}

                CBLAS_GEMV ( CblasRowMajor, CblasTrans, nHid, nVis, 1.0, vishid, nVis, poshidstates, 1, 0.0, negdata, 1 );

                //for(int k=0;k<nVis;k++)
                //    negdata[k] = 1.0 / (1.0 + exp(-(negdata[k]+visbias[k])));
                V_ADD ( nVis, negdata, visbias, negdata );  // negdata[k]+visbias[k]
                V_MULCI ( -1.0, negdata, nVis );             // *(-1)
                V_EXP ( negdata, nVis );                     // exp(x)
                V_ADDCI ( 1.0, negdata, nVis );              // +1
                V_INV ( nVis, negdata, negdata );        // inv(x)

                //neghidprobs = 1./(1 + exp(-negdata*vishid - repmat(hidbiases,numcases,1)));
                //for(int k=0;k<nHid;k++)
                //{
                //    REAL sum = 0.0;
                //    for(int l=0;l<nVis;l++)
                //        sum += negdata[l] * vishid[l + k*nVis];
                //    neghidprobs[k] = 1.0 / (1.0 + exp(-(sum+hidbias[k])));
                //}

                CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, nHid, nVis, 1.0, vishid, nVis, negdata, 1, 0.0, neghidprobs, 1 );

                //for(int k=0;k<nHid;k++)
                //    neghidprobs[k] = 1.0 / (1.0 + exp(-(neghidprobs[k]+hidbias[k])));
                V_ADD ( nHid, neghidprobs, hidbias, neghidprobs );  // neghidprobs[k]+hidbias[k]
                V_MULCI ( -1.0, neghidprobs, nHid );             // *(-1)
                V_EXP ( neghidprobs, nHid );                     // exp(x)
                V_ADDCI ( 1.0, neghidprobs, nHid );              // +1
                V_INV ( nHid, neghidprobs, neghidprobs );        // inv(x)

                //negprods  = negdata'*neghidprobs;
                //for(int k=0;k<nHid;k++)
                //    for(int l=0;l<nVis;l++)
                //        negprods[l+k*nVis] = negdata[l] * neghidprobs[k];

                //neghidact = sum(neghidprobs);
                //for(int k=0;k<nHid;k++)
                //    neghidact[k] = neghidprobs[k];
                V_COPY ( neghidprobs, neghidact, nHid );   // neghidact = neghidprobs

                //negvisact = sum(negdata);
                //REAL negvisact = 0.0;
                //for(int k=0;k<nVis;k++)
                //    negvisact[k] = negdata[k];
                V_COPY ( negdata, negvisact, nVis );   // negvisact = negdata

                //err= sum(sum( (data-negdata).^2 ));
                double err = 0.0;
                for ( int k=0;k<nVis;k++ )
                    err += ( data[k] - negdata[k] ) * ( data[k] - negdata[k] );

                if ( isnan ( err ) || isinf ( err ) )
                    cout<<endl;

                errsum = err + errsum;
                errcnt += nVis;


                // ===== UPDATE WEIGHTS =====
                //vishidinc = momentum*vishidinc + epsilonw*( (posprods-negprods)/numcases - weightcost*vishid);
                //visbiasinc = momentum*visbiasinc + (epsilonvb/numcases)*(posvisact-negvisact);
                //hidbiasinc = momentum*hidbiasinc + (epsilonhb/numcases)*(poshidact-neghidact);
                //vishid = vishid + vishidinc;
                //visbiases = visbiases + visbiasinc;
                //hidbiases = hidbiases + hidbiasinc;

                for ( int k=0;k<nHid;k++ )
                    for ( int l=0;l<nVis;l++ )
                        vishid[l+k*nVis] += m_rbmLearnrateWeights* ( ( data[l] * poshidprobs[k] - negdata[l] * neghidprobs[k] ) - m_rbmWeightDecay*vishid[l+k*nVis] );

                //for(int k=0;k<nHid*nVis;k++)
                //    vishid[k] += m_rbmLearnrateWeights*(posprods[k]-negprods[k] - m_rbmWeightDecay*vishid[k]);

                for ( int k=0;k<nVis;k++ )
                    visbias[k] += m_rbmLearnrateBiasVis * ( posvisact[k]-negvisact[k] );
                for ( int k=0;k<nHid;k++ )
                    hidbias[k] += m_rbmLearnrateBiasHid * ( poshidact[k]-neghidact[k] );

            }
            for ( int sample=0;sample<nSamples;sample++ )
            {
                REAL* data = batchdata + nVis * sample;
                // save for next batch
                CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, nHid, nVis, 1.0, vishid, nVis, data, 1, 0.0, poshidprobs, 1 );
                for ( int k=0;k<nHid;k++ )
                    poshidprobs[k] = 1.0 / ( 1.0 + exp ( - ( poshidprobs[k]+hidbias[k] ) ) );
                REAL* batchPtr = batchdataNext + sample * nHid;
                for ( int k=0;k<nHid;k++ )
                    batchPtr[k] = poshidprobs[k];
            }
            cout<<"epoch:"<<epoch+1<<"/"<<m_rbmMaxEpochs<<"  errsum:"<<errsum<<"  errcnt:"<<errcnt<<"  mse:"<<errsum/ ( double ) errcnt<<"  "<<time ( 0 )-t0<<"[s]"<<endl;
            t0 = time ( 0 );
        }

        int nrWeightsHere = m_nrLayWeights[i];
        weightOffset += nrWeightsHere;
        int offset = weightOffset;
        cout<<"weight offset:"<<offset<<"  offsetNext:"<<weightOffset<<endl;
        REAL* weightPtr = m_weights + offset;

        // copy weights to net
        for ( int j=0;j<nHid;j++ )
            for ( int k=0;k<nVis;k++ )
                m_weights[getWeightIndex ( i + 1, j, k ) ] = vishid[k + j*nVis];
        // copy weights to nnAuto
        REAL* nnAutoWeights = m_nnAuto->getWeightPtr();
        for ( int j=0;j<nHid;j++ )
            for ( int k=0;k<nVis;k++ )
            {
                nnAutoWeights[m_nnAuto->getWeightIndex ( i + 1, j, k ) ] = vishid[k + j*nVis];
                nnAutoWeights[m_nnAuto->getWeightIndex ( 2* ( nLayer-1 ) - i, k, j ) ] = vishid[k + j*nVis];
            }

        // copy biases to net
        for ( int j=0;j<nHid;j++ )
            m_weights[getBiasIndex ( i + 1, j ) ] = hidbias[j];
        // copy biases to nnAuto
        for ( int j=0;j<nHid;j++ )
            nnAutoWeights[m_nnAuto->getBiasIndex ( i + 1, j ) ] = hidbias[j];
        for ( int j=0;j<nVis;j++ )
            nnAutoWeights[m_nnAuto->getBiasIndex ( 2* ( nLayer-1 ) - i, j ) ] = visbias[j];

        /*
        // check it forward calc
        int noutOffset = 0;
        for(int j=0;j<i;j++)
            noutOffset += m_neuronsPerLayer[j] + 1;
        int sample = rand() % nSamples;
        cout<<"sample:"<<sample<<endl;
        // input and data vector
        REAL* data = batchdata + nVis * sample;
        REAL* in = input + m_neuronsPerLayer[0] * sample;
        // print NN outputs
        cout<<"===NN Outputs==="<<endl;
        forwardCalculationBLAS(in);
        int of = m_neuronsPerLayer[i]+1 + noutOffset;
        for(int j=0;j<nHid;j++)
            cout<<m_outputs[j+of]<<" ";
        cout<<endl;
        // print RBM output
        REAL* batchPtr = batchdataNext + sample * nHid;
        cout<<"===Poshidprobs==="<<endl;
        for(int k=0;k<nHid;k++)
            cout<<batchPtr[k]<<" ";
        cout<<endl<<endl;
        */


        /*
        if(i==0) // input layer
        {
            // copy biases
            for(int j=0;j<nHid;j++)
            {
                int ind = j*(nVis+1);
                assert(ind+offset<m_nrWeights);
                weightPtr[ind] = hidbias[j];
            }
            // copy weights
            for(int j=0;j<nHid;j++)
                for(int k=0;k<nVis;k++)
                {
                    int ind = 1 + k + j*(nVis+1);
                    assert(ind+offset<m_nrWeights);
                    weightPtr[ind] = vishid[k + j*nVis];
                }
        }
        else
        {
            // copy biases
            for(int j=0;j<nHid;j++)
            {
                int ind = nVis + j*(nVis+1);
                assert(ind+offset<m_nrWeights);
                weightPtr[ind] = hidbias[j];
            }
            // copy weights
            for(int j=0;j<nHid;j++)
                for(int k=0;k<nVis;k++)
                {
                    int ind = k + j*(nVis+1);
                    assert(ind+offset<m_nrWeights);
                    weightPtr[ind] = vishid[k + j*nVis];
                }
        }
        */

        /*
        int checkEnd = offset + m_nrLayWeights[i+1];
        cout<<"check end:"<<checkEnd<<endl;
        for(int j=0;j<checkEnd;j++)
        if(m_weights[j] == (float)0.01)
        {
            cout<<"j:"<<j<<endl;
            assert(false);
        }
        */

        /*
        // check it forward calc
        noutOffset = 0;
        for(int j=0;j<i;j++)
            noutOffset += m_neuronsPerLayer[j] + 1;
        //for(int sample=0;sample<nSamples&&sample<3;sample++)
        cout<<"sample:"<<sample<<endl;
        {
            // input and data vector
            REAL* data = batchdata + nVis * sample;
            REAL* in = input + m_neuronsPerLayer[0] * sample;

            // print NN outputs
            cout<<"===NN Outputs==="<<endl;
            forwardCalculationBLAS(in);
            int of = m_neuronsPerLayer[i]+1 + noutOffset;
            for(int j=0;j<nHid;j++)
                cout<<m_outputs[j+of]<<" ";
            cout<<endl;

            // print RBM output
            REAL* batchPtr = batchdataNext + sample * nHid;
            cout<<"===Poshidprobs==="<<endl;
            for(int k=0;k<nHid;k++)
                cout<<batchPtr[k]<<" ";

            cout<<endl<<endl;
        }
        */

        delete[] batchdata;
        batchdata = new REAL[nSamples*nHid];
        for ( int j=0;j<nSamples*nHid;j++ )
            batchdata[j] = batchdataNext[j];
        delete[] batchdataNext;
        delete[] vishid;
        delete[] visbias;
        delete[] hidbias;
        delete[] poshidprobs;
        delete[] neghidprobs;
        delete[] hidbiasinc;
        delete[] visbiasinc;
        delete[] poshidstates;
        delete[] negdata;
        delete[] neghidact;
        delete[] negvisact;
        delete[] poshidact;
        delete[] posvisact;
        delete[] probs;
    }

    delete[] batchdata;

    // check if all weights have values
    for ( int i=0;i<m_nnAuto->getNrWeights();i++ )
    {
        REAL* p = m_nnAuto->getWeightPtr();
        if ( p[i] == 0.0 || isnan ( p[i] ) )
        {
            cout<<"p["<<i<<"]:"<<p[i]<<endl;
            assert ( false );
        }
    }

    cout<<"Weights summary from neural net:"<<endl;
    for ( int i=0;i<m_nrLayer;i++ )
    {
        int n = m_neuronsPerLayer[i+1];
        int nIn = m_neuronsPerLayer[i];
        cout<<"#neur:"<<n<<" #neurPrev:"<<nIn<<" ";

        // weights
        REAL min = 1e10, max = -1e10;
        double sum = 0.0, cnt = 0.0, sum2 = 0.0;

        for ( int j=0;j<n;j++ )
            for ( int k=0;k<nIn;k++ )
            {
                REAL v = m_weights[getWeightIndex ( i + 1, j, k ) ];
                if ( min > v )
                    min = v;
                if ( max < v )
                    max = v;
                sum += v;
                sum2 += v*v;
                cnt += 1.0;
            }
        cout<<" weights: cnt|min|max|mean|std: "<<cnt<<"|"<<min<<"|"<<max<<"|"<<sum/cnt<<"|"<<sqrt ( sum2/cnt- ( sum/cnt ) * ( sum/cnt ) ) <<endl;

        min = 1e10;
        max = -1e10;
        sum = 0.0;
        cnt = 0.0;
        sum2 = 0.0;
        // biases
        for ( int j=0;j<n;j++ )
        {
            REAL v = m_weights[getBiasIndex ( i + 1, j ) ];
            if ( min > v )
                min = v;
            if ( max < v )
                max = v;
            sum += v;
            sum2 += v*v;
            cnt += 1.0;
        }
        cout<<"                       biases: cnt|min|max|mean|std: "<<cnt<<"|"<<min<<"|"<<max<<"|"<<sum/cnt<<"|"<<sqrt ( sum2/cnt- ( sum/cnt ) * ( sum/cnt ) ) <<endl;

    }

    // dummy print of middle layer
    /*if ( crossRun == 0 )
    {
        fstream f0 ( "tmp/r0.txt",ios::out );
        REAL* p = new REAL[m_nrInputs];
        for ( int i=0;i<nSamples;i++ )
        {
            m_nnAuto->predictSingleInput ( input + i*m_nrInputs, p );

            for ( int j=0;j<nTarget;j++ )
                f0<<target[j+i*nTarget]<<" ";

            int midLayer = nLayer - 1;
            for ( int j=0;j<m_nnAuto->m_neuronsPerLayer[midLayer];j++ )
                f0<<m_nnAuto->m_outputs[m_nnAuto->getOutputIndex ( midLayer, j ) ]<<" ";
            f0<<endl;

        }
        delete[] p;
        f0.close();
    }*/

    // finetuning of the autoencoder with backprop
    if ( nFinetuningEpochs != -1 )
    {
        cout<<"=== begin: finetuning ==="<<endl;
        m_nnAuto->m_normalTrainStopping = false;
        m_nnAuto->m_maxEpochs = nFinetuningEpochs;
        if ( nFinetuningEpochs > 0 )
            m_nnAuto->trainNN();
        cout<<"=== end: finetuning ==="<<endl;
        /*
        cout<<endl<<"Adjust bias and weights in coding layer to met the mean=0 and std=1 condition"<<endl;
        int nCodeNeurons = m_neuronsPerLayer[nLayer-1];
        cout<<"#codeNeurons:"<<nCodeNeurons<<endl;
        double* mean = new double[nCodeNeurons];
        double* std = new double[nCodeNeurons];
        REAL* p = new REAL[m_nrInputs];
        for(int i=0;i<nCodeNeurons;i++)
        {
            mean[i] = 0.0;
            std[i] = 0.0;
        }
        for(int i=0;i<nSamples;i++)
        {
            m_nnAuto->predictSingleInput(input + i*m_nrInputs, p);
            for(int j=0;j<nCodeNeurons;j++)
            {
                REAL v = m_nnAuto->m_outputs[m_nnAuto->getOutputIndex(nLayer-1, j)];
                mean[j] += v;
                std[j] += v*v;
            }
        }
        for(int i=0;i<nCodeNeurons;i++)
        {
            mean[i] /= (double)nSamples;
            std[i] /= (double)nSamples;
            std[i] = sqrt(std[i] - mean[i]*mean[i]);
        }
        cout<<endl<<"[no correction]:"<<endl;
        for(int i=0;i<nCodeNeurons;i++)
            cout<<"codeLayer "<<i<<":  mean:"<<mean[i]<<"  std:"<<std[i]<<endl;

        // correct the coding weights
        // copy weights : from nnAuto to net
        for(int j=0;j<nCodeNeurons;j++)
            for(int k=0;k<m_neuronsPerLayer[nLayer-2];k++)
                m_nnAuto->m_weights[m_nnAuto->getWeightIndex(nLayer-1,j,k)] /= std[j];
        // copy biases : from nnAuto to net
        //for(int j=0;j<nCodeNeurons;j++)
        //    m_nnAuto->m_weights[getBiasIndex(nLayer-1,j)] -= mean[j]/std[j];

        for(int i=0;i<nCodeNeurons;i++)
        {
            mean[i] = 0.0;
            std[i] = 0.0;
        }
        for(int i=0;i<nSamples;i++)
        {
            m_nnAuto->predictSingleInput(input + i*m_nrInputs, p);
            for(int j=0;j<nCodeNeurons;j++)
            {
                REAL v = m_nnAuto->m_outputs[m_nnAuto->getOutputIndex(nLayer-1, j)];
                mean[j] += v;
                std[j] += v*v;
            }
        }
        for(int i=0;i<nCodeNeurons;i++)
        {
            mean[i] /= (double)nSamples;
            std[i] /= (double)nSamples;
            std[i] = sqrt(std[i] - mean[i]*mean[i]);
        }
        cout<<endl<<"[with corrected std]:"<<endl;
        for(int i=0;i<nCodeNeurons;i++)
            cout<<"codeLayer "<<i<<":  mean:"<<mean[i]<<"  std:"<<std[i]<<endl;

        // correct the coding weights
        // copy weights : from nnAuto to net
        //for(int j=0;j<nCodeNeurons;j++)
        //    for(int k=0;k<m_neuronsPerLayer[nLayer-2];k++)
        //        m_nnAuto->m_weights[m_nnAuto->getWeightIndex(nLayer-1,j,k)] /= std[j];
        // copy biases : from nnAuto to net
        for(int j=0;j<nCodeNeurons;j++)
            m_nnAuto->m_weights[getBiasIndex(nLayer-1,j)] -= mean[j];

        for(int i=0;i<nCodeNeurons;i++)
        {
            mean[i] = 0.0;
            std[i] = 0.0;
        }
        for(int i=0;i<nSamples;i++)
        {
            m_nnAuto->predictSingleInput(input + i*m_nrInputs, p);
            for(int j=0;j<nCodeNeurons;j++)
            {
                REAL v = m_nnAuto->m_outputs[m_nnAuto->getOutputIndex(nLayer-1, j)];
                mean[j] += v;
                std[j] += v*v;
            }
        }
        for(int i=0;i<nCodeNeurons;i++)
        {
            mean[i] /= (double)nSamples;
            std[i] /= (double)nSamples;
            std[i] = sqrt(std[i] - mean[i]*mean[i]);
        }
        cout<<endl<<"[with corrected std and mean]:"<<endl;
        for(int i=0;i<nCodeNeurons;i++)
            cout<<"codeLayer "<<i<<":  mean:"<<mean[i]<<"  std:"<<std[i]<<endl;

        delete[] p;
        delete[] mean;
        delete[] std;
        */

        cout<<endl<<"Copy the autoencoder weights to the net"<<endl;
        cout<<"#weights:"<<m_nrLayWeightOffsets[nLayer]<<endl;
        for ( int i=0;i<nLayer;i++ )
        {
            cout<<"layer "<<i<<": #neurons:"<<m_neuronsPerLayer[i]<<"  #neuronsAuto:"<<m_nnAuto->m_neuronsPerLayer[i]<<endl;
            if ( i > 0 )
            {
                // copy weights : from nnAuto to net
                for ( int j=0;j<m_neuronsPerLayer[i];j++ )
                    for ( int k=0;k<m_neuronsPerLayer[i-1];k++ )
                        m_weights[getWeightIndex ( i,j,k ) ] = m_nnAuto->m_weights[m_nnAuto->getWeightIndex ( i,j,k ) ];
                // copy biases : from nnAuto to net
                for ( int j=0;j<m_neuronsPerLayer[i];j++ )
                    m_weights[getBiasIndex ( i,j ) ] = m_nnAuto->m_weights[m_nnAuto->getBiasIndex ( i,j ) ];
            }
        }

    }

    // dummy print of middle layer
    /*if ( crossRun == 0 )
    {
        fstream f0 ( "tmp/r1.txt",ios::out );
        REAL* p = new REAL[m_nrInputs];
        for ( int i=0;i<nSamples;i++ )
        {
            m_nnAuto->predictSingleInput ( input + i*m_nrInputs, p );

            for ( int j=0;j<nTarget;j++ )
                f0<<target[j+i*nTarget]<<" ";

            int midLayer = nLayer - 1;
            for ( int j=0;j<m_nnAuto->m_neuronsPerLayer[midLayer];j++ )
                f0<<m_nnAuto->m_outputs[m_nnAuto->getOutputIndex ( midLayer, j ) ]<<" ";
            f0<<endl;

        }
        delete[] p;
        f0.close();
    }
    */
}

void NNRBM::printMiddleLayerToFile(string fname, REAL* input, REAL* target, int nSamples, int nTarget)
{
    fstream f0 ( fname.c_str(),ios::out );
    REAL* p = new REAL[m_nrInputs];
    for ( int i=0;i<nSamples;i++ )
    {
        m_nnAuto->predictSingleInput ( input + i*m_nrInputs, p );
        
        for ( int j=0;j<nTarget;j++ )
            f0<<target[j+i*nTarget]<<" ";
        
        int midLayer = m_nrLayer - 1;
        for ( int j=0;j<m_nnAuto->m_neuronsPerLayer[midLayer];j++ )
            f0<<m_nnAuto->m_outputs[m_nnAuto->getOutputIndex ( midLayer, j ) ]<<" ";
        
        //for ( int j=0;j<m_nrInputs;j++ )
        //    f0<<p[j]<<" ";
        
        f0<<endl;
    }
    delete[] p;
    f0.close();
}

void NNRBM::backpropBLAS ( REAL* input, REAL* target, bool* updateLayer )
{
    REAL sum0, d1;

    int outputOffset = m_nrOutputs - m_neuronsPerLayer[m_nrLayer] - 1;  // -1 for bias and output neuron
    int n0 = m_neuronsPerLayer[m_nrLayer-1];
    int outputOffsetPrev = outputOffset - n0 - 1;
    int outputOffsetNext = outputOffset;

    int weightOffset = m_nrWeights - m_nrLayWeights[m_nrLayer];
    int weightOffsetNext, nP1;

    REAL *deltaWPtr, *derivatesPtr, *weightsPtr, *outputsPtr, *d1Ptr, *d1Ptr0, targetConverted, error;

    // ================== the output neuron:  d(j)=(b-o(j))*Aj' ==================
    for ( int i=0;i<m_nrTargets;i++ )
    {
        REAL out = m_outputs[outputOffset+i];

        REAL errorTrain = out * m_scaleOutputs + m_offsetOutputs - target[i];
        m_sumSquaredError += errorTrain * errorTrain;
        m_sumSquaredErrorSamples += 1.0;

        targetConverted = ( target[i] - m_offsetOutputs ) / m_scaleOutputs;

        // cross-entropy error: E = - target * log(out/target)
        // dE/dout = - target * (1/out) * out'
        //if(out < 1e-6)
        //    out = 1e-6;
        //error = - targetConverted / out;

        /*REAL eps = 1e-6;
        out = (out < eps)? eps : out;
        out = (out > 1.0-eps)? 1.0-eps : out;
        error = (out-targetConverted)/(out*(1.0-out));*/

        error = out - targetConverted;

        if ( m_errorFunctionMAE )
            error = error > 0.0? 1.0 : -1.0;

        d1 = error * m_derivates[outputOffset+i];
        m_d1[outputOffset+i] = d1;

        if ( updateLayer )
            if ( updateLayer[m_nrLayer] == false )
                continue;

        deltaWPtr = m_deltaW + weightOffset + i* ( n0+1 );
        if ( m_nrLayer==1 )
        {
            outputsPtr = input - 1;
            deltaWPtr[0] = d1;
            for ( int j=1;j<n0+1;j++ )
                deltaWPtr[j] = d1 * outputsPtr[j];
        }
        else
        {
            outputsPtr = m_outputs + outputOffsetPrev;
            for ( int j=0;j<n0+1;j++ )
                deltaWPtr[j] = d1 * outputsPtr[j];
        }
    }

    // ================== all other neurons in the net ==================
    outputOffsetNext = outputOffset;  // next to current
    outputOffset = outputOffsetPrev;  // current to prev
    n0 = m_neuronsPerLayer[m_nrLayer-2];
    outputOffsetPrev -= n0 + 1;  // prev newnrInputs_
    weightOffset -= m_nrLayWeights[m_nrLayer-1];  // offset to weight pointer
    weightOffsetNext = m_nrWeights - m_nrLayWeights[m_nrLayer];

    for ( int i=m_nrLayer-1;i>0;i-- ) // all layers from output to input
    {
        int n = m_neuronsPerLayer[i];
        int nNext = m_neuronsPerLayer[i+1];
        int nPrev = m_neuronsPerLayer[i-1];
        nP1 = n+1;

        d1Ptr0 = m_d1 + outputOffsetNext;
        derivatesPtr = m_derivates + outputOffset;
        weightsPtr = m_weights + weightOffsetNext;
        d1Ptr = m_d1 + outputOffset;
        deltaWPtr = m_deltaW + weightOffset;
        if ( i==1 )
            outputsPtr = input;
        else
            outputsPtr = m_outputs + outputOffsetPrev;

        // d(j) = SUM(d(k)*w(k,j))
        CBLAS_GEMV ( CblasRowMajor, CblasTrans, nNext, n, 1.0, weightsPtr, nP1, d1Ptr0, 1, 0.0, d1Ptr, 1 );  // d1(j) =W_T*d1(k)
        V_MUL ( n, d1Ptr, derivatesPtr, d1Ptr );

        bool updateWeights = true;
        if ( updateLayer )
            if ( updateLayer[i] == false )
                updateWeights = false;

        // every neuron in the layer calc weight update
        for ( int j=0;j<n;j++ )
        {
            if ( updateWeights )
            {
                if ( i==1 )
                {
                    V_COPY ( outputsPtr, deltaWPtr+1, nPrev );
                    deltaWPtr[0] = 1.0;
                }
                else
                    V_COPY ( outputsPtr, deltaWPtr, nPrev+1 );
                V_MULC ( deltaWPtr, d1Ptr[j], deltaWPtr, nPrev+1 );
            }
            deltaWPtr += nPrev+1;
        }

        outputOffsetNext = outputOffset;  // next to current
        outputOffset = outputOffsetPrev;  // current to prev
        n0 = m_neuronsPerLayer[i-2];
        outputOffsetPrev -= n0 + 1;  // prev new
        weightOffset -= m_nrLayWeights[i-1];  // offset to weight pointer
        weightOffsetNext -= m_nrLayWeights[i];
    }
}

void NNRBM::backprop ( REAL* input, REAL* target, bool* updateLayer )
{
    REAL sum0, d1;

    int outputOffset = m_nrOutputs - m_neuronsPerLayer[m_nrLayer] - 1;  // -1 for bias and output neuron
    int n0 = m_neuronsPerLayer[m_nrLayer-1];
    int outputOffsetPrev = outputOffset - n0 - 1;
    int outputOffsetNext = outputOffset;

    int weightOffset = m_nrWeights - m_nrLayWeights[m_nrLayer];
    int weightOffsetNext, nP1;

    REAL *deltaWPtr, *derivatesPtr, *weightsPtr, *outputsPtr, *d1Ptr, *d1Ptr0, targetConverted, error;

    // ================== the output neuron:  d(j)=(b-o(j))*Aj' ==================
    for ( int i=0;i<m_nrTargets;i++ )
    {
        double out = m_outputs[outputOffset+i];

        REAL errorTrain = out * m_scaleOutputs + m_offsetOutputs - target[i];
        m_sumSquaredError += errorTrain * errorTrain;
        m_sumSquaredErrorSamples += 1.0;

        targetConverted = ( target[i] - m_offsetOutputs ) / m_scaleOutputs;
        error = out - targetConverted;

        if ( m_errorFunctionMAE )
            error = error > 0.0? 1.0 : -1.0;
        d1 = error * m_derivates[outputOffset+i];
        m_d1[outputOffset+i] = d1;

        if ( updateLayer )
            if ( updateLayer[m_nrLayer] == false )
                continue;

        deltaWPtr = m_deltaW + weightOffset + i* ( n0+1 );
        if ( m_nrLayer==1 )
        {
            outputsPtr = input - 1;
            deltaWPtr[0] = d1;
            for ( int j=1;j<n0+1;j++ )
                deltaWPtr[j] = d1 * outputsPtr[j];
        }
        else
        {
            outputsPtr = m_outputs + outputOffsetPrev;
            for ( int j=0;j<n0+1;j++ )
                deltaWPtr[j] = d1 * outputsPtr[j];
        }

    }

    // ================== all other neurons in the net ==================
    outputOffsetNext = outputOffset;  // next to current
    outputOffset = outputOffsetPrev;  // current to prev
    n0 = m_neuronsPerLayer[m_nrLayer-2];
    outputOffsetPrev -= n0 + 1;  // prev newnrInputs_
    weightOffset -= m_nrLayWeights[m_nrLayer-1];  // offset to weight pointer
    weightOffsetNext = m_nrWeights - m_nrLayWeights[m_nrLayer];

    for ( int i=m_nrLayer-1;i>0;i-- ) // all layers from output to input
    {
        int n = m_neuronsPerLayer[i];
        int nNext = m_neuronsPerLayer[i+1];
        int nPrev = m_neuronsPerLayer[i-1];
        nP1 = n+1;

        d1Ptr0 = m_d1 + outputOffsetNext;
        derivatesPtr = m_derivates + outputOffset;
        weightsPtr = m_weights + weightOffsetNext;
        d1Ptr = m_d1 + outputOffset;
        deltaWPtr = m_deltaW + weightOffset;
        if ( i==1 )
            outputsPtr = input - 1;
        else
            outputsPtr = m_outputs + outputOffsetPrev;

        bool updateWeights = true;
        if ( updateLayer )
            if ( updateLayer[i] == false )
                updateWeights = false;

        for ( int j=0;j<n;j++ ) // every neuron in the layer
        {
            // calc d1
            sum0 = 0.0;
            for ( int k=0;k<nNext;k++ ) // all neurons in the next layer:  d(j)=Aj'*Sum(k,d(k)*w(k,j))
                sum0 += d1Ptr0[k] * weightsPtr[k*nP1];
            sum0 *= *derivatesPtr;
            d1Ptr[j] = sum0;

            if ( updateWeights )
            {
                // weight updates
                if ( i==1 )
                {
                    deltaWPtr[0] = sum0;
                    for ( int k=1;k<nPrev+1;k++ )
                        deltaWPtr[k] = sum0 * outputsPtr[k];
                }
                else
                {
                    for ( int k=0;k<nPrev+1;k++ )
                        deltaWPtr[k] = sum0 * outputsPtr[k];
                }
            }

            deltaWPtr += nPrev+1;
            weightsPtr++;
            derivatesPtr++;
        }

        outputOffsetNext = outputOffset;  // next to current
        outputOffset = outputOffsetPrev;  // current to prev
        n0 = m_neuronsPerLayer[i-2];
        outputOffsetPrev -= n0 + 1;  // prev new
        weightOffset -= m_nrLayWeights[i-1];  // offset to weight pointer
        weightOffsetNext -= m_nrLayWeights[i];
    }

}

void NNRBM::forwardCalculationBLAS ( REAL* input )
{
    int outputOffset = m_neuronsPerLayer[0]+1, outputOffsetPrev = 0;
    REAL tmp0, tmp1, sum0;
    REAL *outputPtr, *ptr0, *ptr1, *weightPtr = m_weights;

    for ( int i=0;i<m_nrLayer;i++ ) // to all layer
    {
        int n = m_neuronsPerLayer[i+1];
        int nprev = m_neuronsPerLayer[i] + 1;
        int inputOffset = 0;
        int activationType = m_activationFunctionPerLayer[i+1];
        ptr0 = m_outputs + outputOffset;
        ptr1 = m_derivates + outputOffset;
        if ( i==0 )
        {
            outputPtr = input;
            inputOffset = 1;
        }
        else
            outputPtr = m_outputs + outputOffsetPrev;

        // WeightMatrix*InputVec = Outputs
        CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, n, nprev - inputOffset, 1.0, weightPtr + inputOffset, nprev, outputPtr, 1, 0.0, ptr0, 1 );
        if ( inputOffset )
        {
            for ( int j=0;j<n;j++ )
                ptr0[j] += weightPtr[j * nprev];
        }

        // activation fkt
        if ( activationType == 0 ) // sigmoid: f(x)=1/(1+exp(-x))
        {
            for ( int j=0;j<n;j++ )
                ptr0[j] = -ptr0[j];
            V_EXP ( ptr0, n );
            for ( int j=0;j<n;j++ )
            {
                ptr0[j] = 1.0 / ( 1.0 + ptr0[j] );  // sigmoid
                ptr1[j] = ptr0[j] * ( 1.0 - ptr0[j] );  // derivative
            }
        }
        else if ( activationType == 1 ) // linear: f(x)=x,  df(x)/dx=1
        {
            for ( int j=0;j<n;j++ )
                ptr1[j] = 1.0;
        }
        else if ( activationType == 2 ) // tanh: f(x)=tanh(x)
        {
            V_TANH ( n, ptr0, ptr0 );  // m_outputs = tanh(m_outputs)
            V_SQR ( n, ptr0, ptr1 );  // ptr1[j] = ptr0[j]*ptr0[j]
            for ( int j=0;j<n;j++ )
                ptr1[j] = 1.0 - ptr1[j];
        }
        else if ( activationType == 3 ) // softmax: f(x)=exp(x)/sum(exp(xi))  (where xi are all outputs on this layer)
        {
            V_EXP ( ptr0, n );

            REAL sumOut = 0.0, invSumOut;
            for ( int j=0;j<n;j++ )
                sumOut += ptr0[j];

            invSumOut = 1.0 / sumOut;

            for ( int j=0;j<n;j++ )
                ptr1[j] = ptr0[j] * ( sumOut - ptr0[j] ) * invSumOut * invSumOut;  // derivative

            for ( int j=0;j<n;j++ )
                ptr0[j] *= invSumOut;  // output

            /*
            REAL sumOut = 0.0, invSumOut;
            for(int j=0;j<n;j++)
                sumOut += exp(ptr0[j]);
            invSumOut = 1.0 / sumOut;
            invSumOut = invSumOut * invSumOut;
            for(int j=0;j<n;j++)
                ptr1[j] = invSumOut * ( sumOut * exp(ptr0[j]) - exp(ptr0[j])*exp(ptr0[j]) );  // derivative
            // normalize output: outputs are probabilites
            for(int j=0;j<n;j++)
                ptr0[j] = exp(ptr0[j]) / sumOut;
            */
        }
        else
            assert ( false );

        // update index
        weightPtr += n*nprev;
        outputOffset += n+1;  // this points to first neuron in current layer
        outputOffsetPrev += nprev;  // this points to first neuron in previous layer

    }

}

void NNRBM::forwardCalculation ( REAL* input )
{
    int outputOffset = m_neuronsPerLayer[0] + 1, outputOffsetPrev = 0;
    REAL tmp0, tmp1, sum0;
    REAL *outputPtr, *ptr0, *ptr1, *weightPtr = m_weights;
    for ( int i=0;i<m_nrLayer;i++ ) // to all layer
    {
        int n = m_neuronsPerLayer[i+1];
        int nprev = m_neuronsPerLayer[i] + 1;
        int loopOffset = i==0? 1 : 0;
        int activationType = m_activationFunctionPerLayer[i+1];
        ptr0 = m_outputs + outputOffset;
        ptr1 = m_derivates + outputOffset;
        if ( i==0 )
            outputPtr = input - loopOffset;
        else
            outputPtr = m_outputs + outputOffsetPrev;
        
        for ( int j=0;j<n;j++ ) // all neurons in this layer
        {
            sum0 = i==0? weightPtr[0] : 0.0;  // dot product sum, for inputlayer: init with bias
            for ( int k=loopOffset;k<nprev;k++ ) // calc dot product
                sum0 += weightPtr[k] * outputPtr[k];
            weightPtr += nprev;
            if ( activationType == 0 ) // sigmoid: f(x)=1/(1+exp(-x))
            {
                tmp0 = exp ( -sum0 );
                ptr1[j] = tmp0/ ( ( 1.0 + tmp0 ) * ( 1.0 + tmp0 ) );  // derivation
                ptr0[j] = 1.0 / ( 1.0 + tmp0 );  // activation function
            }
            else if ( activationType == 1 ) // linear: f(x)=x
            {
                ptr0[j] = sum0;
                ptr1[j] = 1.0;
            }
            else if ( activationType == 2 ) // tanh: f(x)=tanh(x)
            {
                tmp0 = tanh ( sum0 );
                ptr0[j] = tmp0;
                ptr1[j] = ( 1.0 - tmp0*tmp0 );
            }
            else
                assert ( false );
        }
        outputOffset += n+1;  // this points to first neuron in current layer
        outputOffsetPrev += nprev;  // this points to first neuron in previous layer
    }
}

void NNRBM::printAutoencoderWeightsToJavascript(string fname)
{
    fstream autoF(fname.c_str(),ios::out);
    autoF<<"function evalAuto(vIn,vOut)"<<endl;
    autoF<<"{"<<endl;
    
    int outputOffset = m_neuronsPerLayer[0] + 1, outputOffsetPrev = 0;
    REAL tmp0, tmp1, sum0;
    REAL *weightPtr = m_weights;
    for ( int i=0;i<m_nrLayer;i++ ) // to all layer
    {
        int n = m_neuronsPerLayer[i+1];
        int nprev = m_neuronsPerLayer[i] + 1;
        int loopOffset = i==0? 1 : 0;
        int activationType = m_activationFunctionPerLayer[i+1];
        
        if(i >= m_nrLayer/2)
        {
            autoF<<endl<<"// layer: "<<i<<endl;
            autoF<<"outVecPrev = new Array("<<nprev<<");"<<endl;
            autoF<<"for(k=0;k<"<<nprev<<";k++)"<<endl;
            if(i == m_nrLayer/2)
                autoF<<"  outVecPrev[k]=vIn[k];"<<endl;
            else
                autoF<<"  outVecPrev[k]=outVec[k];"<<endl;
            autoF<<"outVec = new Array("<<n+1<<");"<<endl;
        }
        
        for ( int j=0;j<n;j++ ) // all neurons in this layer
        {
            if(i >= m_nrLayer/2)
            {
                autoF<<"// neuron: "<<j<<endl;
                autoF<<"w=[";
                for ( int k=loopOffset;k<nprev;k++ )
                    autoF<<weightPtr[k]<<((k+1<nprev)?",":"");
                autoF<<"];"<<endl;
                autoF<<"sum=0.0;"<<endl;
                autoF<<"for(k=0;k<"<<nprev<<";k++)"<<endl;
                autoF<<"  sum +=w[k]*outVecPrev[k];"<<endl;
                autoF<<"outVec["<<j<<"]=1.0/(1.0+Math.exp(-sum));"<<endl;
                
            }
            weightPtr += nprev;
        }
        
        if(i >= m_nrLayer/2)
        {
            autoF<<"outVec["<<n<<"]=1.0;"<<endl;
            
            if(i==m_nrLayer-1)
            {
                autoF<<"//========= copy to output ========="<<endl;
                autoF<<"for(k=0;k<"<<n<<";k++)"<<endl;
                autoF<<"  vOut[k]=outVec[k];"<<endl;
            }
        }
        
        outputOffset += n+1;  // this points to first neuron in current layer
        outputOffsetPrev += nprev;  // this points to first neuron in previous layer
    }
    
    autoF<<"}"<<endl;
    autoF.close();
}

int NNRBM::trainNN()
{
    cout<<"Train the NN with "<<m_nrExamplesTrain<<" samples"<<endl;
    double rmseMin = 1e10, lastUpdate = 1e10, rmseProbeOld = 1e10, rmseTrain = 1e10, rmseProbe = 1e10;
    time_t t0 = time ( 0 );
    while ( 1 )
    {
        rmseProbeOld = rmseProbe;
        rmseTrain = getRMSETrain();
        rmseProbe = getRMSEProbe();
        lastUpdate = rmseProbeOld - rmseProbe;

        cout<<"e:"<<m_globalEpochs<<"  rmseTrain:"<<rmseTrain<<"  rmseProbe:"<<rmseProbe<<"  "<<flush;

        if ( m_normalTrainStopping )
        {
            if ( rmseProbe < rmseMin )
            {
                rmseMin = rmseProbe;
                saveWeights();
            }
            else
            {
                cout<<"rmse rises."<<endl;
                return m_globalEpochs;
            }
            if ( m_minUpdateBound > fabs ( lastUpdate ) )
            {
                cout<<"min update too small (<"<<m_minUpdateBound<<")."<<endl;
                return m_globalEpochs;
            }
        }
        if ( m_maxEpochs == m_globalEpochs )
        {
            cout<<"max epochs reached."<<endl;
            return m_globalEpochs;
        }

        trainOneEpoch();

        cout<<"lRate:"<<m_learnRate<<"  ";
        cout<<time ( 0 )-t0<<"[s]"<<endl;
        t0 = time ( 0 );
    }
    return -1;
}

void NNRBM::trainOneEpoch ( bool* updateLayer )
{
    int batchCnt = 0;
    V_ZERO ( m_weightsBatchUpdate, m_nrWeights );
    m_sumSquaredError = 0.0;
    m_sumSquaredErrorSamples = 0.0;

    for ( int i=0;i<m_nrExamplesTrain;i++ )
    {
        if ( updateLayer != 0 )
            V_ZERO ( m_deltaW, m_nrWeights );  // clear the weight update, it might be that not every layer gets updates

        REAL* inputPtr = m_inputsTrain + i * m_nrInputs;
        REAL* targetPtr = m_targetsTrain + i * m_nrTargets;

        // forward
        if ( m_useBLAS )
            forwardCalculationBLAS ( inputPtr );
        else
            forwardCalculation ( inputPtr );

        // backward: calc weight update
        if ( m_useBLAS )
            backpropBLAS ( inputPtr, targetPtr, updateLayer );
        else
            backprop ( inputPtr, targetPtr, updateLayer );

        // accumulate the weight updates
        if ( m_batchSize > 1 )
            V_ADD ( m_nrWeights, m_deltaW, m_weightsBatchUpdate, m_weightsBatchUpdate );

        batchCnt++;

        // if batch size is reached, or the last element in training list
        if ( batchCnt >= m_batchSize || i == m_nrExamplesTrain - 1 )
        {
            // batch init
            batchCnt = 0;
            if ( m_batchSize > 1 )
            {
                V_COPY ( m_weightsBatchUpdate, m_deltaW, m_nrWeights ); // deltaW = weightsBatchUpdate
                V_ZERO ( m_weightsBatchUpdate, m_nrWeights );
            }

            if ( m_enableRPROP )
            {
                // weight update:
                // deltaW = {  if dE/dW>0  then  -adaptiveRPROPlRate
                //             if dE/dW<0  then  +adaptiveRPROPlRate
                //             if dE/dW=0  then   0  }
                // learnrate adaption:
                // adaptiveRPROPlRate = {  if (dE/dW_old * dE/dW)>0  then  adaptiveRPROPlRate*RPROP_etaPos
                //                         if (dE/dW_old * dE/dW)<0  then  adaptiveRPROPlRate*RPROP_etaNeg
                //                         if (dE/dW_old * dE/dW)=0  then  adaptiveRPROPlRate  }
                REAL dW, dWOld, sign, update, prod;
                for ( int j=0;j<m_nrWeights;j++ )
                {
                    dW = m_deltaW[j];
                    dWOld = m_deltaWOld[j];
                    prod = dW * dWOld;
                    sign = dW > 0.0? 1.0 : -1.0;
                    if ( prod > 0.0 )
                    {
                        m_adaptiveRPROPlRate[j] *= m_RPROP_etaPos;
                        if ( m_adaptiveRPROPlRate[j] > m_RPROP_updateMax )
                            m_adaptiveRPROPlRate[j] = m_RPROP_updateMax;
                        update = sign * m_adaptiveRPROPlRate[j];
                        m_weights[j] -= update + m_weightDecay * m_weights[j];   // weight update and weight decay
                        m_deltaWOld[j] = dW;
                    }
                    else if ( prod < 0.0 )
                    {
                        m_adaptiveRPROPlRate[j] *= m_RPROP_etaNeg;
                        if ( m_adaptiveRPROPlRate[j] < m_RPROP_updateMin )
                            m_adaptiveRPROPlRate[j] = m_RPROP_updateMin;
                        m_deltaWOld[j] = 0.0;
                    }
                    else // prod == 0.0
                    {
                        update = sign * m_adaptiveRPROPlRate[j];
                        m_weights[j] -= update + m_weightDecay * m_weights[j];   // weight update and weight decay
                        m_deltaWOld[j] = dW;
                    }
                }
            }
            else  // stochastic gradient decent (batch-size: m_batchSize)
            {
                //=========== slower stochastic updates (without vector libraries) ===========
                // update weights + weight decay
                // formula: weights -= eta * (dE(w)/dw + lambda * w)
                //if(m_momentum > 0.0)
                //{
                //    for(int j=0;j<m_nrWeights;j++)
                //        m_weightsOldOld[j] = m_weightsOld[j];
                //    for(int j=0;j<m_nrWeights;j++)
                //        m_weightsOld[j] = m_weights[j];
                //}
                //
                //if(m_momentum > 0.0)
                //{
                //    for(int j=0;j<m_nrWeights;j++)
                //    {
                //        m_weightsTmp0[j] = (1.0 - m_momentum)*m_weightsTmp0[j] + m_momentum*m_weightsTmp1[j];
                //        m_weightsTmp1[j] = m_weightsTmp0[j];
                //        m_weights[j] -= m_weightsTmp0[j];
                //    }
                //}
                //for(int j=0;j<m_nrWeights;j++)
                //    m_weights[j] -= m_learnRate * (m_deltaW[j] + m_weightDecay * m_weights[j]);


                // update weights + weight decay(L2 reg.)
                // formula: weights = weights - eta * (dE(w)/dw + lambda * weights)
                //          weights = weights - (eta*deltaW + eta*lambda*weights)
                for ( int j=0;j<m_nrWeights;j++ )
                    if ( isnan ( m_deltaW[j] ) || isinf ( m_deltaW[j] ) )
                    {
                        cout<<endl<<"value:"<<m_deltaW[j]<<endl<<"j:"<<j<<" i:"<<i<<endl;
                        assert ( false );
                    }
                V_MULC ( m_deltaW, m_learnRate, m_weightsTmp0, m_nrWeights );     // tmp0 = learnrate * deltaW

                // if weight decay enabled
                if ( m_weightDecay > 0.0 )
                {
                    if ( m_enableL1Regularization )
                    {
                        // update weights + L1 reg.
                        // formula: weights = weights - eta * (dE(w)/dw + lambda*sign(w))
                        //          weights = weights - (eta*deltaW + eta*lambda*sign(w))
                        for ( int j=0;j<m_nrWeights;j++ )
                            m_weightsTmp2[j] = 1.0;
                        for ( int j=0;j<m_nrWeights;j++ )
                            if ( m_weights[j]<0.0 )
                                m_weightsTmp2[j] = -1.0;
                        //REAL c = m_weightDecay * m_learnRate;
                        //for(int j=0;j<m_nrWeights;j++)
                        //    m_weightsTmp0[j] += c * m_weightsTmp2[j];
                        CBLAS_AXPY ( m_nrWeights, m_weightDecay * m_learnRate, m_weightsTmp2, 1, m_weightsTmp0, 1 );  // tmp0 = reg*learnrate*weights+tmp0
                    }
                    else
                        //saxpy(n, a, x, incx, y, incy)
                        //y := a*x + y
                        CBLAS_AXPY ( m_nrWeights, m_weightDecay * m_learnRate, m_weights, 1, m_weightsTmp0, 1 );  // tmp0 = reg*learnrate*weights+tmp0
                }

                // if momentum is used
                if ( m_momentum > 0.0 )
                {
                    V_MULC ( m_weightsTmp0, 1.0 - m_momentum, m_weightsTmp0, m_nrWeights ); // tmp0 = tmp0 * (1 - momentum)   [actual update]
                    V_MULC ( m_weightsTmp1, m_momentum, m_weightsTmp1, m_nrWeights );    // tmp1 = tmp1 * momentum         [last update]

                    // sum updates
                    V_ADD ( m_nrWeights, m_weightsTmp0, m_weightsTmp1, m_weightsTmp0 ); // tmp0 = tmp0 + tmp1

                    V_COPY ( m_weightsTmp0, m_weightsTmp1, m_nrWeights );               // tmp1 = tmp0
                }

                // standard weight update in the NNRBM
                V_SUB ( m_nrWeights, m_weights, m_weightsTmp0, m_weights );       // weights = weights - tmp0
            }
        }

        // make the learnrate smaller (per sample)
        m_learnRate -= m_learnrateDecreaseRate;
        if ( m_learnRate < m_learnRateMin )
            m_learnRate = m_learnRateMin;

    }

    // make the learnrate smaller (per epoch)
    m_learnRate -= m_learnrateDecreaseRateEpoch;
    if ( m_learnRate < m_learnRateMin )
        m_learnRate = m_learnRateMin;

    // epoch counter
    m_globalEpochs++;
}

void NNRBM::setGlobalEpochs ( int e )
{
    m_globalEpochs = e;
}

void NNRBM::printLearnrate()
{
    cout<<"lRate:"<<m_learnRate<<" "<<flush;
}

void NNRBM::saveWeights()
{
    // nothing here
}

void NNRBM::predictSingleInput ( REAL* input, REAL* output )
{
    REAL* inputPtr = input;

    // forward
    if ( m_useBLAS )
        forwardCalculationBLAS ( inputPtr );
    else
        forwardCalculation ( inputPtr );

    // output correction
    REAL* outputPtr = m_outputs + m_nrOutputs - m_nrTargets - 1;
    for ( int i=0;i<m_nrTargets;i++ )
        output[i] = outputPtr[i] * m_scaleOutputs + m_offsetOutputs;
}

REAL NNRBM::calcRMSE ( REAL* inputs, REAL* targets, int examples )
{
    double rmse = 0.0;
    for ( int i=0;i<examples;i++ )
    {
        REAL* inputPtr = inputs + i * m_nrInputs;
        REAL* targetPtr = targets + i * m_nrTargets;

        predictSingleInput ( inputPtr, m_outputsTmp );

        for ( int j=0;j<m_nrTargets;j++ )
            rmse += ( m_outputsTmp[j] - targetPtr[j] ) * ( m_outputsTmp[j] - targetPtr[j] );
    }
    rmse = sqrt ( rmse/ ( double ) ( examples*m_nrTargets ) );
    return rmse;
}

REAL NNRBM::getRMSETrain()
{
    return calcRMSE ( m_inputsTrain, m_targetsTrain, m_nrExamplesTrain );
}

REAL NNRBM::getRMSEProbe()
{
    return calcRMSE ( m_inputsProbe, m_targetsProbe, m_nrExamplesProbe );
}

REAL* NNRBM::getWeightPtr()
{
    return m_weights;
}

void NNRBM::setWeights ( REAL* w )
{
    cout<<"Set new weights"<<endl;
    for ( int i=0;i<m_nrWeights;i++ )
        m_weights[i] = w[i];
}

int NNRBM::getNrWeights()
{
    return m_nrWeights;
}