nn.cpp

#include "nn.h"

extern StreamOutput cout;

NN::NN()
{
    // 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_sumSquaredError = 0.0;
    m_sumSquaredErrorSamples = 0;
    m_nrLayWeightOffsets = 0;
    m_sumSquaredError = 0.0;
    m_sumSquaredErrorSamples = 0.0;
    m_activationFunctionType = 0;
}

NN::~NN()
{
    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_nrLayWeightOffsets )
        delete[] m_nrLayWeightOffsets;
    m_nrLayWeightOffsets = 0;
}

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

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


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

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

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

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

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

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

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

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


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

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

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

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

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

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

void NN::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 NN::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 NN::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 NN::setNormalTrainStopping ( bool en )
{
    m_normalTrainStopping = en;
    cout<<"normalTrainStopping: "<<m_normalTrainStopping<<endl;
}

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

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

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

void NN::setActivationFunctionType( int type )
{
    if(type==0)
    {
        cout<<"activationFunctionType: tanh"<<endl;
        m_activationFunctionType = 0;
    }
    else if(type==1)
    {
        cout<<"activationFunctionType: sin"<<endl;
        m_activationFunctionType = 1;
    }
    else if(type==2)
    {
        cout<<"activationFunctionType: tanhMod0"<<endl;
        m_activationFunctionType = 2;
    }
    else
        assert(false);
}

int NN::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 NN::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 NN::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;
}

void NN::setNNStructure ( int nrLayer, int* neuronsPerLayer )
{
    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;

    // 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 NN::getInitWeight ( int fanIn )
{
    double nr = 2.0* ( rand() / ( double ) RAND_MAX-0.5 );  // -1 .. +1
    return ( 1.0/sqrt ( ( double ) fanIn ) ) * nr;
}

void NN::initNNWeights ( time_t seed )
{
    srand ( seed );
    cout<<"init weights ";
    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 NN::backpropBLAS ( REAL* input, REAL* target )
{
    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;
        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;
        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 );

        // every neuron in the layer calc weight update
        for ( int j=0;j<n;j++ )
        {
            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 NN::backprop ( REAL* input, REAL* target )
{
    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++;

        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;
        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;

        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;

            // 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 NN::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;
        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];
        }

        if(m_activationFunctionType == 0)
        {
            // activation fkt: f(x)=tanh(x)
            V_TANH ( n, ptr0, ptr0 );  // m_outputs = tanh(m_outputs)
            V_SQR ( n, ptr0, ptr1 );   // m_derivates = tanh(m_outputs) * tanh(m_outputs)
            for ( int j=0;j<n;j++ )
                ptr1[j] = 1.0 - ptr1[j];
        }
        else if(m_activationFunctionType == 1)
        {
            // activation fkt: f(x)=sin(x)+0.01*x
            REAL piHalf = 1.570796326794897;
            for(int j=0;j<n;j++)
            {
                REAL v = ptr0[j], sign = v>0.0? 1.0 : -1.0;
                if(v > -piHalf && v < piHalf)
                {
                    ptr0[j] = sin(v) + v * 0.01;
                    ptr1[j] = cos(v) + sign * 0.01;
                }
                else  // sumWeights is outside a half periode +/-pi/2
                {
                    ptr0[j] = sign + v * 0.01;
                    ptr1[j] = sign * 0.01;
                }
            }
        }
        else if(m_activationFunctionType == 2)
        {
            // activation fkt: f(x)= wenn x>0: f(x)=x^(1+tanh(v)*mul)
            //                       wenn x<0: f(x)=-(-x)^(1+tanh(-v)*mul)
            REAL mul = 0.5;  // 0.25 : swing: [-1.5...+1.5]
            for(int j=0;j<n;j++)
            {
                REAL v = ptr0[j], tanhV = tanh(v), tanhVNeg = -tanhV;
                if(v >= 0.0)
                {
                    ptr0[j] = pow(v,1.0+tanhV*mul); //pow(v,0.3);
                    ptr1[j] = pow(v,tanhV*mul)*(1.0+tanhV*mul)+ptr0[j]*log(v)*mul*(1.0-tanhV*tanhV);
                    if(isnan(ptr1[j]) || isinf(ptr1[j]))
                        ptr1[j] = 1.0;
                }
                else
                {
                    ptr0[j] = -pow(-v,1.0+tanhVNeg*mul); //-pow(-v,0.3);
                    ptr1[j] = -pow(-v,tanhVNeg*mul)*(1.0+tanhVNeg*mul)*(-1.0)+ptr0[j]*log(-v)*mul*(1.0-tanhV*tanhV)*(-1.0);
                    if(isnan(ptr1[j]) || isinf(ptr1[j]))
                        ptr1[j] = -1.0;
                }
            }
        }
        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 NN::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;
        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(m_activationFunctionType == 0)
            {
                // activation fkt: f(x)=tanh(x)
                tmp0 = tanh ( sum0 );
                ptr0[j] = tmp0;
                ptr1[j] = ( 1.0 - tmp0*tmp0 );
            }
            else if(m_activationFunctionType == 1)
            {
                // activation fkt: f(x)=sin(x)+0.01*x
                REAL piHalf = 1.570796326794897;
                REAL v = ptr0[j], sign = v>0.0? 1.0 : -1.0;
                if(v > -piHalf && v < piHalf)
                {
                    ptr0[j] = sin(v) + v * 0.01;
                    ptr1[j] = cos(v) + sign * 0.01;
                }
                else  // sumWeights is outside a half periode +/-pi/2
                {
                    ptr0[j] = sign + v * 0.01;
                    ptr1[j] = sign * 0.01;
                }
            }
            else if(m_activationFunctionType == 2)
            {
                // activation fkt: f(x)= wenn x>0: f(x)=x^(1+tanh(v)*mul)
                //                       wenn x<0: f(x)=-(-x)^(1+tanh(-v)*mul)
                REAL mul = 0.5;  // 0.25 : swing: [-1.5...+1.5]   // 0.5 : swing: [-1.18195...+1.18195]
                REAL v = ptr0[j], tanhV = tanh(v), tanhVNeg = -tanhV;
                if(v >= 0.0)
                {
                    ptr0[j] = pow(v,1.0+tanhV*mul); //pow(v,0.3);
                    ptr1[j] = pow(v,tanhV*mul)*(1.0+tanhV*mul)+ptr0[j]*log(v)*mul*(1.0-tanhV*tanhV);
                    if(isnan(ptr1[j]) || isinf(ptr1[j]))
                        ptr1[j] = 1.0;
                }
                else
                {
                    ptr0[j] = -pow(-v,1.0+tanhVNeg*mul); //-pow(-v,0.3);
                    ptr1[j] = -pow(-v,tanhVNeg*mul)*(1.0+tanhVNeg*mul)*(-1.0)+ptr0[j]*log(-v)*mul*(1.0-tanhV*tanhV)*(-1.0);
                    if(isnan(ptr1[j]) || isinf(ptr1[j]))
                        ptr1[j] = -1.0;
                }
            }
            else
                assert(false);
        }
        outputOffset += n+1;  // this points to first neuron in current layer
        outputOffsetPrev += nprev;  // this points to first neuron in previous layer
    }

}

int NN::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 NN::trainOneEpoch()
{
    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++ )
    {
        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 );
        else
            backprop ( inputPtr, targetPtr );

        // 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)
                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 NN
                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 NN::setGlobalEpochs ( int e )
{
    m_globalEpochs = e;
}

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

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

void NN::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 NN::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 NN::getRMSETrain()
{
    return calcRMSE ( m_inputsTrain, m_targetsTrain, m_nrExamplesTrain );
}

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

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

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

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