InputFeatureSelector.cpp

#include "InputFeatureSelector.h"

extern StreamOutput cout;

InputFeatureSelector::InputFeatureSelector()
{
    cout<<"InputFeatureSelector"<<endl;
}

InputFeatureSelector::~InputFeatureSelector()
{
    cout<<"descructor InputFeatureSelector"<<endl;
}

void InputFeatureSelector::predictWithRidgeRegression ( int* colSelect, int nCol, int totalCol, int nTarget, REAL* ATAfull, REAL* ATbfull, REAL* testFeat, REAL* output, int nPred, REAL* ATA, REAL* ATb, REAL* AbTrans )
{
    // fill ATA matrix
    for ( int i=0;i<nCol;i++ )
        for ( int j=0;j<nCol;j++ )
            ATA[i*nCol + j] = ATAfull[totalCol*colSelect[i] + colSelect[j]];

    // fill ATb matrix
    for ( int i=0;i<nCol;i++ )
        for ( int j=0;j<nTarget;j++ )
            ATb[i*nTarget + j] = ATbfull[colSelect[i]*nTarget + j];

    //int m = nCol, k = nTarget;
    for ( int i=0;i<nTarget;i++ )
        for ( int j=0;j<nCol;j++ )
            AbTrans[j+i*nCol] = ATb[i+j*nTarget];  // copy to colMajor

    int lda = nCol, nrhs = nTarget, ldb = nCol;
    int info;

    //Cholesky factorization of a symmetric (Hermitian) positive-definite matrix
    LAPACK_POTRF ( "U", &nCol, ATA, &lda, &info );

    // Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite matrix
    LAPACK_POTRS ( "U", &nCol, &nrhs, ATA, &lda, AbTrans, &ldb, &info );  // Y=solution vector

    if ( info != 0 )
        cout<<"error equation solver failed to converge: "<<info<<endl;

    // predict test set
    for ( int i=0;i<nTarget;i++ )
    {
        REAL* xPtr = AbTrans + i*nCol;
        for ( int j=0;j<nPred;j++ ) //fold*testSize;j<(fold+1)*testSize;j++)
        {
            REAL sum = 0.0;
            for ( int a=0;a<nCol;a++ )
                sum += xPtr[a] * testFeat[colSelect[a] + j*totalCol];
            output[i+j*nTarget] = sum;
        }
    }
}

void InputFeatureSelector::forwardSelection ( bool* selection, REAL* inputs, int nFeatures, int nExamples, int* labels, REAL* targets, int nClass, int nDomain )
{
    cout<<endl<<"---------------------------------- FORWARD GREEDY FEATURE SELECTION ----------------------------------"<<endl;

    int minFeatures = 1000, maxFeatures = 2000;
    int nCross = 8;
    float percentTrain = 0.5;  // [%] train set

    int testBegin = nExamples * percentTrain;
    cout<<"testBegin:"<<testBegin<<" nExamples:"<<nExamples<<endl;
    int trainSize = testBegin;
    int testSize = nExamples - testBegin;
    cout<<"testSize:"<<testSize<<" trainSize:"<<trainSize<<endl;
    int* trainSplits = new int[nCross+1];
    int* testSplits = new int[nCross+1];
    trainSplits[0] = 0;
    testSplits[0] = testBegin;
    trainSplits[nCross] = testBegin;
    testSplits[nCross] = nExamples;
    REAL trainSlotSize = trainSize / ( REAL ) nCross;
    REAL testSlotSize = testSize / ( REAL ) nCross;
    for ( int i=0;i<nCross;i++ )
    {
        trainSplits[i+1] = trainSlotSize * ( i+1 );
        testSplits[i+1] = testSlotSize * ( i+1 ) + testBegin;
    }

    bool optimizeRMSE = true;
    REAL reg = 1.0e-1;
    cout<<"reg:"<<reg<<endl;

    // precalculated matrices for cross-fold validation on train and testset
    REAL** ATATrainFull = new REAL*[nCross], **ATbTrainFull = new REAL*[nCross];
    REAL** ATATestFull = new REAL*[nCross], **ATbTestFull = new REAL*[nCross];
    for ( int i=0;i<nCross;i++ )
    {
        ATATrainFull[i] = new REAL[nFeatures*nFeatures];
        ATbTrainFull[i] = new REAL[nFeatures*nClass*nDomain];
        ATATestFull[i] = new REAL[nFeatures*nFeatures];
        ATbTestFull[i] = new REAL[nFeatures*nClass*nDomain];
    }

    // mark already used features
    bool* usedFeatures = new bool[nFeatures];
    int* bestFeatures = new int[nFeatures];
    for ( int i=0;i<nFeatures;i++ )
    {
        usedFeatures[i] = false;
        bestFeatures[i] = -1;
    }
    int bestFeaturesCnt = 0;

    // add const (last input dimension)
    cout<<"Add feature nr. "<<nFeatures-1<<" (added const. 1 as last input dimension)"<<endl;
    bestFeatures[0] = nFeatures-1;
    usedFeatures[nFeatures-1] = true;
    //featureIDFile<<bestFeatures[bestFeaturesCnt]<<endl;
    bestFeaturesCnt++;

    cout<<"Precalculate ATA and ATb for every cross-fold (trainset): "<<flush;

    REAL* inputsTmp = new REAL[nExamples*nFeatures];
    REAL* targetsTmp = new REAL[nExamples*nClass*nDomain];
    for ( int i=0;i<nCross;i++ )
    {
        cout<<"."<<flush;

        // ================================== TRAINSET ==================================
        memcpy ( inputsTmp, inputs, sizeof ( REAL ) *nExamples*nFeatures );
        memcpy ( targetsTmp, targets, sizeof ( REAL ) *nExamples*nClass*nDomain );

        // train data: set validation elements to zero
        memset ( inputsTmp + nFeatures * trainSplits[i], 0, sizeof ( REAL ) *nFeatures* ( trainSplits[i+1]-trainSplits[i] ) );
        memset ( targetsTmp + nClass * nDomain * trainSplits[i], 0, sizeof ( REAL ) *nClass*nDomain* ( trainSplits[i+1]-trainSplits[i] ) );

        // train data: set testset elements to zero
        memset ( inputsTmp + nFeatures * testBegin, 0, sizeof ( REAL ) *nFeatures*testSize );
        memset ( targetsTmp + nClass * nDomain * testBegin, 0, sizeof ( REAL ) *nClass*nDomain*testSize );

        int nTrain = 0;
        for ( int j=0;j<nExamples;j++ )
        {
            int cnt = 0;
            for ( int k=0;k<nFeatures;k++ )
                if ( inputsTmp[j*nFeatures+k] == 0.0 )
                    cnt++;
            if ( cnt != nFeatures )
                nTrain++;
        }
        cout<<" [nTrain:"<<nTrain<<"]";

        int n = nExamples, m = nFeatures, k = nClass*nDomain;

        REAL* A = inputsTmp;
        // A'*A
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, ATATrainFull[i], m ); // AA = A'A

        // add lambda
        for ( int j=0;j<m;j++ )
            ATATrainFull[i][j*m + j] += ( ( REAL ) trainSize - ( REAL ) trainSize/ ( REAL ) nCross ) * reg;

        // A'*b
        REAL *b = targetsTmp;
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, k , n, 1.0, A, m, b, k, 0.0, ATbTrainFull[i], k ); // Ab = A'b
    }
    cout<<endl;

    cout<<"Precalculate ATA and ATb for every cross-fold (testset): "<<flush;
    for ( int i=0;i<nCross;i++ )
    {
        cout<<"."<<flush;

        // ================================== TESTSET ==================================
        memcpy ( inputsTmp, inputs, sizeof ( REAL ) *nExamples*nFeatures );
        memcpy ( targetsTmp, targets, sizeof ( REAL ) *nExamples*nClass*nDomain );

        // test data: set validation elements to zero
        memset ( inputsTmp + nFeatures * testSplits[i], 0, sizeof ( REAL ) *nFeatures* ( testSplits[i+1]-testSplits[i] ) );
        memset ( targetsTmp + nClass * nDomain * testSplits[i], 0, sizeof ( REAL ) *nClass*nDomain* ( testSplits[i+1]-testSplits[i] ) );

        // test data: set trainset elements to
        memset ( inputsTmp, 0, sizeof ( REAL ) *nFeatures*trainSize );
        memset ( targetsTmp, 0, sizeof ( REAL ) *nClass*nDomain*trainSize );

        int nTest = 0;
        for ( int j=0;j<nExamples;j++ )
        {
            int cnt = 0;
            for ( int k=0;k<nFeatures;k++ )
                if ( inputsTmp[j*nFeatures+k] == 0.0 )
                    cnt++;
            if ( cnt != nFeatures )
                nTest++;
        }
        cout<<" [nTest:"<<nTest<<"]";

        int n = nExamples, m = nFeatures, k = nClass*nDomain;

        REAL* A = inputsTmp;
        // A'*A
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, ATATestFull[i], m ); // AA = A'A

        // add lambda
        for ( int j=0;j<m;j++ )
            ATATestFull[i][j*m + j] += ( ( REAL ) testSize - ( REAL ) testSize/ ( REAL ) nCross ) * reg;

        // A'*b
        REAL *b = targetsTmp;
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, k , n, 1.0, A, m, b, k, 0.0, ATbTestFull[i], k ); // Ab = A'b
    }
    cout<<endl;

    delete[] inputsTmp;
    delete[] targetsTmp;


    REAL rmseTestBest = 1e10;

    // greedy through all features
    for ( int outerLoop=bestFeaturesCnt;outerLoop<nFeatures;outerLoop++ )
    {
        time_t t0 = time ( 0 );
        int progress = nFeatures/10 + 1;

        int bestInd = -1;
        REAL bestTrainErr = 1e10;

        REAL* prediction = new REAL[nExamples*nClass*nDomain];
        REAL* ATA = new REAL[ ( bestFeaturesCnt+1 ) * ( bestFeaturesCnt+1 ) ];
        REAL *ATb = new REAL[ ( bestFeaturesCnt+1 ) *nClass*nDomain];
        REAL* AbTrans = new REAL[ ( bestFeaturesCnt+1 ) *nClass*nDomain];

        cout<<"nF:"<<bestFeaturesCnt<<" "<<flush;

        // search over all features
        for ( int innerLoop=0;innerLoop<nFeatures;innerLoop++ )
        {
            if ( innerLoop%progress==0 )
                cout<<"."<<flush;

            if ( usedFeatures[innerLoop] == true )
                continue;

            // select current feature
            bestFeatures[bestFeaturesCnt] = innerLoop;
            bestFeaturesCnt++;

            // cross-fold train prediction
            for ( int f=0;f<nCross;f++ )
            {
                int nPred = trainSplits[f+1] - trainSplits[f];
                int nBest = bestFeaturesCnt;
                REAL* predPtr = prediction + trainSplits[f] * nClass * nDomain;
                REAL* featPtr = inputs + trainSplits[f] * nFeatures;
                predictWithRidgeRegression ( bestFeatures, nBest, nFeatures, nClass*nDomain, ATATrainFull[f], ATbTrainFull[f], featPtr, predPtr, nPred, ATA, ATb, AbTrans );
            }


            // calc train error
            double rmseTrain = 0.0;
            int rmseTrainCnt = 0;
            for ( int i=0;i<testBegin;i++ )
            {
                for ( int j=0;j<nClass*nDomain;j++ )
                {
                    double err = prediction[j + i*nClass*nDomain] - targets[j + i*nClass*nDomain];
                    rmseTrain += err*err;
                    rmseTrainCnt++;
                }
            }
            rmseTrain = sqrt ( rmseTrain/ ( double ) rmseTrainCnt );
            /*AUC auc;
            REAL* tmp = new REAL[trainSize*nClass*nDomain];
            int* tmpL = new int[trainSize*nDomain];
            for(int i=0;i<trainSize;i++)
            {
                for(int j=0;j<nClass*nDomain;j++)
                    tmp[i+j*trainSize] = prediction[j+i*nClass*nDomain];
                for(int j=0;j<nDomain;j++)
                    tmpL[j+i*nDomain] = labels[j+i*nDomain];
            }
            rmseTrain = auc.getAUC(tmp, tmpL, nClass, nDomain, trainSize);
            delete[] tmp;
            delete[] tmpL;
            */
            if ( rmseTrain < bestTrainErr )
            {
                bestTrainErr = rmseTrain;
                bestInd = innerLoop;
            }

            // deselect current feature
            bestFeaturesCnt--;

        }

        // add best found feature
        usedFeatures[bestInd] = true;
        bestFeatures[bestFeaturesCnt] = bestInd;
        bestFeaturesCnt++;

        // cross-fold test prediction
        for ( int f=0;f<nCross;f++ )
        {
            int nPred = testSplits[f+1] - testSplits[f];
            int nBest = bestFeaturesCnt;
            REAL* predPtr = prediction + testSplits[f] * nClass * nDomain;
            REAL* featPtr = inputs + testSplits[f] * nFeatures;
            predictWithRidgeRegression ( bestFeatures, nBest, nFeatures, nClass*nDomain, ATATrainFull[f], ATbTrainFull[f], featPtr, predPtr, nPred, ATA, ATb, AbTrans );
        }

        // calc test error
        double rmseTest = 0.0;
        int rmseTestCnt = 0;
        for ( int i=testBegin;i<nExamples;i++ )
        {
            for ( int j=0;j<nClass*nDomain;j++ )
            {
                double err = prediction[j + i*nClass*nDomain] - targets[j + i*nClass*nDomain];
                rmseTest += err*err;
                rmseTestCnt++;
            }
        }
        rmseTest = sqrt ( rmseTest/ ( double ) rmseTestCnt );
        /*AUC auc;
        REAL* tmp = new REAL[testSize*nClass*nDomain];
        int* tmpL = new int[testSize*nDomain];
        for(int i=0;i<testSize;i++)
        {
            for(int j=0;j<nClass*nDomain;j++)
                tmp[i+j*testSize] = prediction[j+i*nClass*nDomain + trainSize*nClass*nDomain];
            for(int j=0;j<nDomain;j++)
                tmpL[j+i*nDomain] = labels[j+i*nDomain+trainSize*nDomain];
        }
        rmseTest = auc.getAUC(tmp, tmpL, nClass, nDomain, testSize);
        delete[] tmp;
        delete[] tmpL;*/

        delete[] ATA;
        delete[] ATb;
        delete[] AbTrans;
        delete[] prediction;

        cout<<" err:"<<bestTrainErr<<" (testErr:"<<rmseTest<<")"<<" (bestInd:"<<bestInd<<") "<<time ( 0 )-t0<<"[s]";

        if ( rmseTest < rmseTestBest )
        {
            rmseTestBest = rmseTest;
            cout<<" [best test]";

            // write best features (not the first one!)
            // open the feature file and write best (test min.) found features
            fstream featureIDFile;
            featureIDFile.open ( FEATURE_TXT_FILE,ios::out );
            for ( int i=1;i<nFeatures;i++ )
                if ( bestFeatures[i] != -1 )
                    featureIDFile<<bestFeatures[i]<<endl;
            featureIDFile.close();
        }
        else if ( bestFeaturesCnt >= minFeatures )
            return;
        if ( bestFeaturesCnt >= maxFeatures )
            return;

        cout<<endl;
    }



}

void InputFeatureSelector::backwardElimination ( bool* selection, REAL* inputs, int nFeatures, int nExamples, int* labels, REAL* targets, int nClass, int nDomain )
{
    assert ( false );
}

void InputFeatureSelector::selectFeatures ( bool* selection, REAL* inputs, int nFeatures, int nExamples, int* labels, REAL* targets, int nClass, int nDomain )
{
    REAL* inputsWithConst = new REAL[ ( nFeatures+1 ) *nExamples];
    for ( int i=0;i<nExamples;i++ )
    {
        for ( int j=0;j<nFeatures;j++ )
            inputsWithConst[i* ( nFeatures+1 ) +j] = inputs[i*nFeatures+j];
        inputsWithConst[i* ( nFeatures+1 ) +nFeatures] = 1.0;
    }

    forwardSelection ( selection, inputsWithConst, nFeatures+1, nExamples, labels, targets, nClass, nDomain );
    delete[] inputsWithConst;
    return;
}