BlendStopping.cpp

#include "BlendStopping.h"

extern StreamOutput cout;

BlendStopping::BlendStopping()
{
    cout<<"Constructor BlendStopping"<<endl;

    // init member vars
    m_newPrediction = 0;
    m_nPredictors = 0;
    m_nClass = 0;
    m_nDomain = 0;
    m_nTrain = 0;
    m_blendingRegularization = 0;
    m_classificationError = 0;
    m_singlePrediction = 0;
    m_prediction = 0;
    m_A = 0;
    m_b = 0;
    m_x = 0;
    m_xBest = 0;
    m_data = 0;
    m_paramsAUC = 0;
    m_optimizeRMSE = 0;
    m_globalWeights = 0;
}

BlendStopping::BlendStopping ( Algorithm *data )
{
    cout<<"Constructor BlendStopping"<<endl;
    // init member vars
    m_newPrediction = 0;
    m_nPredictors = 0;
    m_nClass = 0;
    m_nDomain = 0;
    m_nTrain = 0;
    m_blendingRegularization = 0;
    m_classificationError = 0;
    m_singlePrediction = 0;
    m_prediction = 0;
    m_A = 0;
    m_b = 0;
    m_x = 0;
    m_xBest = 0;
    m_data = 0;
    m_paramsAUC = 0;
    m_optimizeRMSE = 0;

    m_data = data;
    m_nClass = m_data->m_nClass;
    m_nDomain = m_data->m_nDomain;
    m_nTrain = m_data->m_nTrain;
    m_globalWeights = m_data->m_enableGlobalBlendingWeights;
}

BlendStopping::BlendStopping ( Algorithm *data, string algorithmFullPrediction )
{
    cout<<"Constructor BlendStopping"<<endl;
    // init member vars
    m_newPrediction = 0;
    m_nPredictors = 0;
    m_nClass = 0;
    m_nDomain = 0;
    m_nTrain = 0;
    m_blendingRegularization = 0;
    m_classificationError = 0;
    m_singlePrediction = 0;
    m_prediction = 0;
    m_A = 0;
    m_b = 0;
    m_x = 0;
    m_xBest = 0;
    m_data = 0;
    m_paramsAUC = 0;
    m_optimizeRMSE = 0;

    m_data = data;
    m_nClass = m_data->m_nClass;
    m_nDomain = m_data->m_nDomain;
    m_nTrain = m_data->m_validationType=="ValidationSet"? m_data->m_validSize : m_data->m_nTrain;
    m_globalWeights = m_data->m_enableGlobalBlendingWeights;

    m_algorithmFullPrediction = algorithmFullPrediction;
    string directory = m_data->m_datasetPath + "/" + m_data->m_fullPredPath + "/";

    vector<string> files = m_data->m_algorithmNameList;//Data::getDirectoryFileList(directory);

    // add a prediction to blend, if the file exists
    string additionalFullPrediction = m_data->m_datasetPath + "/" + m_data->m_fullPredPath + "/" + m_algorithmFullPrediction;
    fstream f ( additionalFullPrediction.c_str(), ios::in );
    if ( f.is_open() )
    {
        cout<<"ADD:"<<additionalFullPrediction<<" ";
        files.push_back ( additionalFullPrediction );
        m_algorithmFullPrediction = "";
    }
    f.close();

    m_nPredictors = 0;
    for ( int i=0;i<files.size();i++ )
    {
        if ( files[i].at ( files[i].size()-1 ) != '.' && files[i].find ( ".dat" ) == files[i].length()-4 )
        {
            //cout<<endl<<"File: "<<files[i]<<endl<<endl;
            m_usedFiles.push_back ( files[i] );
            m_nPredictors++;
        }
    }

    // add the constant predictor
    m_nPredictors++;

    // add the new prediction
    if ( m_algorithmFullPrediction=="tune" )
    {
        m_nPredictors++;
        cout<<"Number of predictors for blendStopping: "<<m_nPredictors<<" (+1 const, +1 new)"<<endl;
    }
    else
        cout<<"Number of predictors for blendStopping: "<<m_nPredictors<<" (+1 const)"<<endl;

    // allocate per class a linear equation system
    // find for every target class the best linear combination
    m_singlePrediction = new REAL[m_nTrain];
    m_prediction = new REAL[m_nTrain*m_nClass*m_nDomain];
    m_A = new REAL*[m_nClass*m_nDomain];
    m_b = new REAL*[m_nClass*m_nDomain];
    m_x = new REAL*[m_nClass*m_nDomain];
    m_xBest = new REAL*[m_nClass*m_nDomain];
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
    {
        m_A[i] = new REAL[m_nPredictors * m_nTrain];
        m_b[i] = new REAL[m_nTrain];
        m_x[i] = new REAL[m_nPredictors];
        m_xBest[i] = new REAL[m_nPredictors];
        for ( int j=0;j<m_nPredictors;j++ )
        {
            m_x[i][j] = 0.0;
            m_xBest[i][j] = 0.0;
        }
    }

    // set target values
    for ( int i=0;i<m_nTrain;i++ )
        for ( int j=0;j<m_nClass*m_nDomain;j++ )
        {
            REAL t = m_data->m_trainTargetOrig[i*m_nClass*m_nDomain + j];
            if(m_data->m_validationType == "ValidationSet")
                t = m_data->m_validTarget[i*m_nClass*m_nDomain + j];
            
            m_b[j][i] = t;
        }

    // set constant predictor as first predictor
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
        for ( int j=0;j<m_nTrain;j++ )
            m_A[i][j] = 1.0;

    // load the fullPredictors
    cout<<endl;
    for ( int i=0;i<m_usedFiles.size();i++ )
    {
        fstream f ( m_usedFiles[i].c_str(), ios::in );
        if ( f.is_open() == false )
        {
            cout<<"Cannot open:"<<m_usedFiles[i]<<endl;
            assert ( false );
        }
        REAL* cache0 = new REAL[m_nTrain*m_nClass*m_nDomain];
        f.read ( ( char* ) cache0, sizeof ( REAL ) *m_nTrain*m_nClass*m_nDomain );
        f.close();
        double rmse0 = 0.0, err;
        for ( int j=0;j<m_nClass*m_nDomain;j++ )
        {
            for ( int k=0;k<m_nTrain;k++ )
            {
                m_A[j][k + ( i+1 ) *m_nTrain] = cache0[j + k*m_nClass*m_nDomain];
                err = m_A[j][k + ( i+1 ) *m_nTrain] - m_b[j][k];
                rmse0 += err * err;
            }
        }
        cout<<"File:"<<m_usedFiles[i]<<"  RMSE:"<<sqrt ( rmse0/ ( double ) ( m_nClass*m_nTrain*m_nDomain ) ) <<endl;
        if ( cache0 )
            delete[] cache0;
        cache0 = 0;
    }
    cout<<endl;

    // init random prediction
    if ( m_algorithmFullPrediction=="tune" ) // HACK
    {
        for ( int i=0;i<m_nClass*m_nDomain;i++ )
            for ( int j=0;j<m_nTrain;j++ )
                m_A[i][ ( m_nPredictors - 1 ) *m_nTrain + j] = 0.0;//NumericalTools::getNormRandomNumber(0.0,0.01);

        // the pointers to the prediction for outside use
        m_newPrediction = new REAL*[m_nClass*m_nDomain];
        for ( int i=0;i<m_nClass*m_nDomain;i++ )
            m_newPrediction[i] = m_A[i] + ( m_nPredictors - 1 ) *m_nTrain;
    }

    if ( m_data->m_errorFunction=="AUC" )
        m_paramsAUC = new double[m_nClass*m_nPredictors*m_nDomain];
}

BlendStopping::~BlendStopping()
{
    cout<<"destructor BlendStopping"<<endl;
    if ( m_singlePrediction )
        delete[] m_singlePrediction;
    m_singlePrediction = 0;
    if ( m_prediction )
        delete[] m_prediction;
    m_prediction = 0;
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
    {
        if ( m_A )
        {
            if ( m_A[i] )
                delete[] m_A[i];
            m_A[i] = 0;
        }
        if ( m_b )
        {
            if ( m_b[i] )
                delete[] m_b[i];
            m_b[i] = 0;
        }
        if ( m_x )
        {
            if ( m_x[i] )
                delete[] m_x[i];
            m_x[i] = 0;
        }
        if ( m_xBest )
        {
            if ( m_xBest[i] )
                delete[] m_xBest[i];
            m_xBest[i] = 0;
        }
    }
    if ( m_A )
        delete[] m_A;
    m_A = 0;
    if ( m_b )
        delete[] m_b;
    m_b = 0;
    if ( m_x )
        delete[] m_x;
    m_x = 0;
    if ( m_xBest )
        delete[] m_xBest;
    m_xBest = 0;
    if ( m_newPrediction )
        delete[] m_newPrediction;
    m_newPrediction = 0;
    if ( m_paramsAUC )
        delete[] m_paramsAUC;
    m_paramsAUC = 0;
}

void BlendStopping::saveTmpBestWeights()
{
    cout<<"[SB]";
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
        for ( int j=0;j<m_nPredictors;j++ )
            m_xBest[i][j] = m_x[i][j];
}

double BlendStopping::calcBlending ( char* fname )
{
    cout<<" [CalcBlend] ";
    double rmse = 0.0, mae = 0.0, err;
    double* rmseClass = new double[m_nClass*m_nDomain];

    // search for optimal lambda in ridge regression
    if ( m_data->m_blendingEnableCrossValidation == true && m_data->m_blendingAlgorithm != "Average" && m_data->m_blendingAlgorithm != "TakeLast" )
    {
        if ( m_data->m_errorFunction=="AUC" && Framework::getDatasetType() && m_data->m_blendingAlgorithm=="NelderMead" )
        {
            // init weights with constants
            for ( int i=0;i<m_nClass*m_nDomain;i++ )
                for ( int j=0;j<m_nPredictors;j++ )
                    m_x[i][j] = m_xBest[i][j] = ( j+1==m_nPredictors ) ? 0.1 : 1.0;

            for ( int i=0;i<m_nClass*m_nDomain;i++ )
            {
                for ( int j=0;j<m_nPredictors;j++ )
                {
                    char buf[1024];
                    sprintf ( buf,"w-t%d-p%d",i,j );
                    m_paramsAUC[i*m_nPredictors+j] = m_x[i][j];
                    addDoubleParameter ( & ( m_paramsAUC[i*m_nPredictors+j] ), buf, -100.0, 100.0 );
                }
            }

            // start the nelder-mead searcher
            setDebug ( 0 );
            //if(fname)
            //    setDebug(1);
            NelderMeadSearch ( 500 );
            //expSearcher(0,100,3,2,0.99,true,false);

            int cnt = 0;
            for ( int i=0;i<m_nClass*m_nDomain;i++ )
            {
                for ( int j=0;j<m_nPredictors;j++ )
                {
                    char buf[1024];
                    sprintf ( buf,"w-t%d-p%d",i,j );

                    m_x[i][j] = m_paramsAUC[i*m_nPredictors+j];
                    m_xBest[i][j] = m_paramsAUC[i*m_nPredictors+j];
                    removeDoubleParameter ( buf );
                }
            }
        }
        else
        {
            setDebug ( 0 );
            addDoubleParameter ( &m_blendingRegularization, "lambda", 1e-10, 1e3 );

            // start the structured searcher
            int minTuninigEpochs = 0, maxTuninigEpochs = 30;
            m_optimizeRMSE = 1;
            expSearcher ( minTuninigEpochs, maxTuninigEpochs, 3, 1, 0.8, true, false );
            m_optimizeRMSE = 0;

            removeDoubleParameter ( "lambda" );

            for ( int i=0;i<m_nClass*m_nDomain;i++ )
                for ( int j=0;j<m_nPredictors;j++ )
                    m_xBest[i][j] = m_x[i][j];

            if ( fname == 0 )
            {
                //cout<<"lambda:"<<m_blendingRegularization<<" ";
                //return m_expSearchErrorBest;
            }
            else
                cout<<endl<<"Blending-CV-error:"<<m_expSearchErrorBest<<endl;
        }
    }

    cout<<"lambda:"<<m_blendingRegularization<<" ";

    for ( int i=0;i<m_nClass*m_nDomain;i++ )
    {
        rmseClass[i] = 0.0;

        if ( m_data->m_blendingAlgorithm=="Average" )
        {
            // weight all algorithms equal (constant-prediction: weight=0)
            m_x[i][0] = 0.0;
            REAL normalizer = 1.0 / ( double ) ( m_nPredictors-1 );
            for ( int j=1;j<m_nPredictors;j++ )
                m_x[i][j] = normalizer;
        }
        else if ( m_data->m_blendingAlgorithm=="TakeLast" )
        {
            // take only the last algorithm in the ensemble
            for ( int j=0;j<m_nPredictors;j++ )
                m_x[i][j] = 0.0;
            m_x[i][m_nPredictors-1] = 1.0;
        }
        else if ( m_data->m_blendingAlgorithm=="LinearRegression" )
        {
            if ( m_globalWeights ) // equal weights for all classes and domains
            {
                REAL* A = new REAL[m_nTrain*m_nPredictors*m_nClass*m_nDomain], *b = new REAL[m_nTrain*m_nClass*m_nDomain];
                for ( int j=0;j<m_nPredictors;j++ )
                    for ( int k=0;k<m_nClass*m_nDomain;k++ )
                        for ( int l=0;l<m_nTrain;l++ )
                            A[l + k*m_nTrain + j*m_nTrain*m_nClass*m_nDomain] = m_A[k][l+j*m_nTrain];
                for ( int k=0;k<m_nClass*m_nDomain;k++ )
                    for ( int l=0;l<m_nTrain;l++ )
                        b[l + k*m_nTrain] = m_b[k][l];
                m_numTools.RidgeRegressionMultisolutionSinglecallGELSS ( A, b, m_x[i], m_nTrain*m_nClass*m_nDomain, m_nPredictors, 1, m_blendingRegularization, false );
                delete[] A;
                delete[] b;
            }
            else
                m_numTools.RidgeRegressionMultisolutionSinglecallGELSS ( m_A[i], m_b[i], m_x[i], m_nTrain, m_nPredictors, 1, m_blendingRegularization, false );
        }
        else if ( m_data->m_blendingAlgorithm=="LinearRegressionNonNeg" )
        {
            if ( m_globalWeights ) // equal weights for all classes and domains
            {
                REAL* A = new REAL[m_nTrain*m_nPredictors*m_nClass*m_nDomain], *b = new REAL[m_nTrain*m_nClass*m_nDomain];
                for ( int j=0;j<m_nPredictors;j++ )
                    for ( int k=0;k<m_nClass*m_nDomain;k++ )
                        for ( int l=0;l<m_nTrain;l++ )
                            A[l + k*m_nTrain + j*m_nTrain*m_nClass*m_nDomain] = m_A[k][l+j*m_nTrain];
                for ( int k=0;k<m_nClass*m_nDomain;k++ )
                    for ( int l=0;l<m_nTrain;l++ )
                        b[l + k*m_nTrain] = m_b[k][l];
                m_numTools.RidgeRegressionNonNegSinglecall ( A, b, m_x[i], m_nTrain, m_nPredictors, m_blendingRegularization, 1e-10, 1000, false );
                delete[] A;
                delete[] b;
            }
            else
                m_numTools.RidgeRegressionNonNegSinglecall ( m_A[i], m_b[i], m_x[i], m_nTrain, m_nPredictors, m_blendingRegularization, 1e-10, 1000, false );

            // in case of zero weights
            REAL sumX = 0.0;
            for ( int j=1;j<m_nPredictors;j++ )
                sumX += m_x[i][j];
            if ( sumX == 0.0 )
            {
                if ( i > 0 )
                    for ( int j=0;j<m_nPredictors;j++ )
                        m_x[i][j] = m_x[i-1][j];
                else
                {
                    for ( int j=0;j<m_nPredictors;j++ )
                        m_x[i][j] = 0.0;
                    m_x[i][m_nPredictors-1] = 1.0;
                }
            }
        }
        else if ( m_data->m_blendingAlgorithm=="NelderMead" )
        {
            //for(int j=0;j<m_nPredictors;j++)
            //    m_x[i][j] = m_xBest[i][j];
            //m_x[i][m_nPredictors-1] = 1.0;
        }

        // calc RMSE with BLAS
        CBLAS_GEMM ( CblasColMajor, CblasNoTrans, CblasNoTrans, m_nTrain, 1, m_nPredictors , 1.0, m_A[i], m_nTrain, m_x[i], m_nPredictors, 0.0, m_singlePrediction, m_nTrain );  // Ab = A'b

        if ( m_data->m_enablePostBlendClipping )
        {
            IPPS_THRESHOLD ( m_singlePrediction, m_singlePrediction, m_nTrain, m_data->m_negativeTarget, ippCmpLess );  // clip negative
            IPPS_THRESHOLD ( m_singlePrediction, m_singlePrediction, m_nTrain, m_data->m_positiveTarget, ippCmpGreater );  // clip positive
        }

        memcpy ( ( char* ) ( m_prediction+m_nTrain*i ), m_singlePrediction, sizeof ( REAL ) *m_nTrain );
        for ( int j=0;j<m_nTrain;j++ )
        {
            rmse += ( m_singlePrediction[j] - m_b[i][j] ) * ( m_singlePrediction[j] - m_b[i][j] );
            mae += fabs ( m_singlePrediction[j] - m_b[i][j] );
            rmseClass[i] += ( m_singlePrediction[j] - m_b[i][j] ) * ( m_singlePrediction[j] - m_b[i][j] );
        }
    }

    //for ( int i=0;i<m_nClass*m_nDomain;i++ )
    //    for ( int j=0;j<m_nPredictors;j++ )
    //        m_xBest[i][j] = m_x[i][j];

    if ( fname )
        printWeights ( false );

    if ( fname )
    {
        cout<<"[Write train prediction:"<<fname<<"]"<<" nSamples:"<<m_nTrain<<endl;
        fstream f ( fname, ios::out );
        if ( f.is_open() == false )
            assert ( false );
        REAL* tmpOut = new REAL[m_nTrain*m_nClass*m_nDomain];
        for ( int i=0;i<m_nTrain;i++ )
        {
            int realIndex = m_data->m_mixDatasetIndices[i];
            if(m_data->m_validationType=="ValidationSet")
                realIndex = i;
            for ( int j=0;j<m_nClass*m_nDomain;j++ )
                tmpOut[j*m_nTrain + realIndex] = m_prediction[j*m_nTrain + i];
        }
        for ( int i=0;i<m_nTrain;i++ )
            for ( int j=0;j<m_nClass*m_nDomain;j++ )
            {
                float predictionSP = tmpOut[j*m_nTrain + i];
                f.write ( ( char* ) &predictionSP,sizeof ( float ) );
            }
        f.close();
        if ( tmpOut )
            delete[] tmpOut;
        tmpOut = 0;
    }

    // only if the dataset is for classification
    if ( Framework::getDatasetType() )
    {
        int* wrongLabelCnt = new int[m_nDomain];
        for ( int d=0;d<m_nDomain;d++ )
            wrongLabelCnt[d] = 0;

        // in all domains
        for ( int d=0;d<m_nDomain;d++ )
        {
            // classification error
            for ( int i=0;i<m_nTrain;i++ )
            {
                REAL max = -1e10;
                int indBest = -1;
                for ( int j=0;j<m_nClass;j++ )
                {
                    if ( max < m_prediction[d*m_nClass*m_nTrain + j*m_nTrain + i] )
                    {
                        max = m_prediction[d*m_nClass*m_nTrain + j*m_nTrain + i];
                        indBest = j;
                    }
                }
                if (m_data->m_validationType == "ValidationSet")
                {
                    if ( indBest != m_data->m_validLabel[d + i*m_nDomain] )
                        wrongLabelCnt[d]++;
                }
                else
                {
                    if ( indBest != m_data->m_trainLabelOrig[d + i*m_nDomain] )
                        wrongLabelCnt[d]++;
                }
            }
        }

        int nWrong = 0;
        for ( int d=0;d<m_nDomain;d++ )
        {
            nWrong += wrongLabelCnt[d];
            //if(m_nDomain > 1)
            //    cout<<"["<<(double)wrongLabelCnt[d]/(double)m_nTrain<<"] ";
        }
        m_classificationError = ( double ) nWrong/ ( ( double ) m_nTrain*m_nDomain );
        cout<<" [classErr:"<<100.0*m_classificationError<<"%] ";

        delete[] wrongLabelCnt;
    }

    if ( rmseClass )
        delete[] rmseClass;
    rmseClass = 0;

    rmse = sqrt ( rmse/ ( ( double ) m_nTrain * ( double ) m_nClass * ( double ) m_nDomain ) );
    mae = mae/ ( ( double ) m_nTrain * ( double ) m_nClass * ( double ) m_nDomain );

    if ( m_data->m_errorFunction=="AUC" && Framework::getDatasetType() )
    {
        cout<<"[rmse:"<<rmse<<"] ";

        // calc area under curve
        REAL auc = getAUC ( m_prediction, m_data->m_trainLabelOrig, m_data->m_nClass, m_data->m_nDomain, m_data->m_nTrain );
        return auc;
    }
    else if ( m_data->m_errorFunction=="MAE" )
        return mae;

    return rmse;
}

double BlendStopping::getClassificationError()
{
    return m_classificationError;
}

void BlendStopping::saveBlendingWeights ( string path, bool saveBest )
{
    printWeights ( saveBest );

    char buf[1024];
    sprintf ( buf, "blendingWeights_%02d.dat", m_nPredictors );
    string name = path + "/" + string ( buf );
    cout<<"Save blending weights: "<<name<<endl;
    fstream f ( name.c_str(), ios::out );
    f.write ( ( char* ) &m_nClass, sizeof ( int ) );
    f.write ( ( char* ) &m_nDomain, sizeof ( int ) );
    f.write ( ( char* ) &m_nPredictors, sizeof ( int ) );

    // write the solution weights
    REAL** xPtr = (saveBest==true?m_xBest:m_x);
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
        f.write ( ( char* ) xPtr[i], sizeof ( REAL ) *m_nPredictors );

    f.close();
}

void BlendStopping::loadBlendingWeights ( string path, int nPredictors )
{
    char buf[1024];
    sprintf ( buf, "blendingWeights_%02d.dat", nPredictors );
    string name = path + "/" + string ( buf );

    cout<<"Load blending weights: "<<name<<endl;
    fstream f ( name.c_str(), ios::in );
    if ( f.is_open() == false )
        assert ( false );
    f.read ( ( char* ) &m_nClass, sizeof ( int ) );
    f.read ( ( char* ) &m_nDomain, sizeof ( int ) );
    f.read ( ( char* ) &m_nPredictors, sizeof ( int ) );

    // alloc mem
    if ( m_x == 0 )
    {
        m_x = new REAL*[m_nClass*m_nDomain];
        for ( int i=0;i<m_nClass*m_nDomain;i++ )
            m_x[i] = new REAL[m_nPredictors];
    }

    // read the solution weights
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
        f.read ( ( char* ) m_x[i], sizeof ( REAL ) *m_nPredictors );

    printWeights ( false );

    f.close();
}

void BlendStopping::predictEnsembleOutput ( REAL** predictions, REAL* output )
{
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
        output[i] = 0.0;
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
    {
        for ( int j=0;j<m_nPredictors;j++ )
            output[i] += m_x[i][j] * predictions[j][i];

        // no clipping, because of ambiguty of max. values (to determine classID)
        if ( m_data->m_enablePostBlendClipping )
            output[i] = NumericalTools::clipValue ( output[i], m_data->m_negativeTarget, m_data->m_positiveTarget );
    }
}

void BlendStopping::printWeights ( bool printBest )
{
    cout<<"Blending weights (row: classes, col: predictors[1.col=const predictor])"<<endl;
    REAL** xPtr = printBest==true?m_xBest : m_x;
    for ( int i=0;i<m_nClass*m_nDomain;i++ )
    {
        if ( cout.m_enableOutput )
        {
            for ( int j=0;j<m_nPredictors;j++ )
                printf ( "%1.3f\t",xPtr[i][j] );
            printf ( "\n" );
        }
    }
}

void BlendStopping::setRegularization ( REAL reg )
{
    cout<<"Blending regularization: "<<reg<<endl;
    m_blendingRegularization = reg;
}


double BlendStopping::calcRMSEonProbe()
{
    if ( m_optimizeRMSE )
    {
        int *randomIndex = new int[m_nTrain];
        REAL* tmpA = new REAL[m_nTrain * m_nPredictors];
        REAL* tmpB = new REAL[m_nTrain];

        // fill random 0..NCross-1  numbers in array
        for ( int i=0;i<m_nTrain;i++ )
        {
            if(m_data->m_validationType=="Bagging")
                randomIndex[i] = rand() % m_data->m_nCross;  // this is random
            else
                randomIndex[i] = m_data->m_crossIndex[i];  // this is taken from the predefined CV-split
        }
        
        double rmse = 0.0;
        int cnt = 0;

        for ( int cross=0;cross<m_data->m_nCross;cross++ )
        {
            for ( int i=0;i<m_nClass*m_nDomain;i++ )
            {
                // clear and save the values of A and b
                for ( int j=0;j<m_nTrain;j++ )
                {
                    if ( randomIndex[j] == cross )
                    {
                        for ( int k=0;k<m_nPredictors;k++ )
                        {
                            tmpA[j + k*m_nTrain] = m_A[i][j + k*m_nTrain];
                            m_A[i][j + k*m_nTrain] = 0.0;
                        }
                        tmpB[j] = m_b[i][j];
                        m_b[i][j] = 0.0;
                    }
                }

                // solve linear equation system
                if ( m_data->m_blendingAlgorithm=="LinearRegression" || m_data->m_blendingAlgorithm=="NelderMead" )
                {
                    if ( m_globalWeights ) // equal weights for all classes and domains
                    {
                        REAL* A = new REAL[m_nTrain*m_nPredictors*m_nClass*m_nDomain], *b = new REAL[m_nTrain*m_nClass*m_nDomain];
                        for ( int j=0;j<m_nPredictors;j++ )
                            for ( int k=0;k<m_nClass*m_nDomain;k++ )
                                for ( int l=0;l<m_nTrain;l++ )
                                    A[l + k*m_nTrain + j*m_nTrain*m_nClass*m_nDomain] = m_A[k][l+j*m_nTrain];
                        for ( int k=0;k<m_nClass*m_nDomain;k++ )
                            for ( int l=0;l<m_nTrain;l++ )
                                b[l + k*m_nTrain] = m_b[k][l];
                        m_numTools.RidgeRegressionMultisolutionSinglecallGELSS ( A, b, m_x[i], m_nTrain*m_nClass*m_nDomain, m_nPredictors, 1, m_blendingRegularization, false );
                        delete[] A;
                        delete[] b;
                    }
                    else
                        m_numTools.RidgeRegressionMultisolutionSinglecallGELSS ( m_A[i], m_b[i], m_x[i], m_nTrain, m_nPredictors, 1, m_blendingRegularization, false );
                }
                else if ( m_data->m_blendingAlgorithm=="LinearRegressionNonNeg" )
                {
                    if ( m_globalWeights ) // equal weights for all classes and domains
                    {
                        REAL* A = new REAL[m_nTrain*m_nPredictors*m_nClass*m_nDomain], *b = new REAL[m_nTrain*m_nClass*m_nDomain];
                        for ( int j=0;j<m_nPredictors;j++ )
                            for ( int k=0;k<m_nClass*m_nDomain;k++ )
                                for ( int l=0;l<m_nTrain;l++ )
                                    A[l + k*m_nTrain + j*m_nTrain*m_nClass*m_nDomain] = m_A[k][l+j*m_nTrain];
                        for ( int k=0;k<m_nClass*m_nDomain;k++ )
                            for ( int l=0;l<m_nTrain;l++ )
                                b[l + k*m_nTrain] = m_b[k][l];
                        m_numTools.RidgeRegressionNonNegSinglecall ( A, b, m_x[i], m_nTrain, m_nPredictors, m_blendingRegularization, 1e-10, 1000, false );
                        delete[] A;
                        delete[] b;
                    }
                    else
                        m_numTools.RidgeRegressionNonNegSinglecall ( m_A[i], m_b[i], m_x[i], m_nTrain, m_nPredictors, m_blendingRegularization, 1e-10, 1000, false );
                }

                // restore the values of A and b
                for ( int j=0;j<m_nTrain;j++ )
                {
                    if ( randomIndex[j] == cross )
                    {
                        for ( int k=0;k<m_nPredictors;k++ )
                            m_A[i][j + k*m_nTrain] = tmpA[j + k*m_nTrain];
                        m_b[i][j] = tmpB[j];
                    }
                }

                // calc RMSE with BLAS
                CBLAS_GEMM ( CblasColMajor, CblasNoTrans, CblasNoTrans, m_nTrain, 1, m_nPredictors , 1.0, m_A[i], m_nTrain, m_x[i], m_nPredictors, 0.0, m_singlePrediction, m_nTrain );  // Ab = A'b
                for ( int j=0;j<m_nTrain;j++ )
                {
                    if ( randomIndex[j] == cross ) // predict only ratings in this fold
                    {
                        if ( m_data->m_enablePostBlendClipping )
                        {
                            if ( m_singlePrediction[j] < m_data->m_negativeTarget )
                                m_singlePrediction[j] = m_data->m_negativeTarget;
                            if ( m_singlePrediction[j] > m_data->m_positiveTarget )
                                m_singlePrediction[j] = m_data->m_positiveTarget;
                        }
                        m_prediction[j + m_nTrain*i] = m_singlePrediction[j];
                        rmse += ( m_singlePrediction[j] - m_b[i][j] ) * ( m_singlePrediction[j] - m_b[i][j] );
                        cnt++;
                    }
                }

            }
        }
        if ( randomIndex )
            delete[] randomIndex;
        randomIndex = 0;
        if ( tmpA )
            delete[] tmpA;
        tmpA = 0;
        if ( tmpB )
            delete[] tmpB;
        tmpB = 0;
        rmse = sqrt ( rmse/ ( double ) cnt );
        return rmse;
    }
    if ( m_optimizeRMSE == 0 ) // optimize AUC
    {
        // init blending weights with params
        for ( int i=0;i<m_nClass*m_nDomain;i++ )
            for ( int j=0;j<m_nPredictors;j++ )
                m_x[i][j] = m_paramsAUC[j + i*m_nPredictors];

        // calc predictions
        for ( int i=0;i<m_nClass*m_nDomain;i++ )
            for ( int j=0;j<m_nTrain;j++ )
            {
                REAL out = 0.0;
                for ( int k=0;k<m_nPredictors;k++ )
                    out += m_A[i][j+k*m_nTrain] * m_x[i][k];
                m_prediction[j+i*m_nTrain] = out;
            }

        // calc area under curve
        REAL auc = getAUC ( m_prediction, m_data->m_trainLabelOrig, m_data->m_nClass, m_data->m_nDomain, m_data->m_nTrain );
        return auc;
    }

    assert ( false );
    return 0;
}

double BlendStopping::calcRMSEonBlend()
{
    assert ( false );
    return 0.0;
}

void BlendStopping::saveBestPrediction()
{

}