PolynomialRegression.cpp

#include "PolynomialRegression.h"

extern StreamOutput cout;

PolynomialRegression::PolynomialRegression()
{
    cout<<"PolynomialRegression"<<endl;
    // init member vars
    m_x = 0;
    m_reg = 0;
    m_polyMean = 0;
    m_polyStd = 0;
    m_polyOrder = 0;
    m_enableCrossInteractions = 0;
    m_inputDim = 0;
}

PolynomialRegression::~PolynomialRegression()
{
    cout<<"descructor PolynomialRegression"<<endl;
    for ( int i=0;i<m_nCross+1;i++ )
    {
        if ( m_x )
        {
            if ( m_x[i] )
                delete[] m_x[i];
            m_x[i] = 0;
        }
    }
    if ( m_x )
        delete[] m_x;
    m_x = 0;
    if ( m_polyMean )
        delete[] m_polyMean;
    m_polyMean = 0;
    if ( m_polyStd )
        delete[] m_polyStd;
    m_polyStd = 0;
}

void PolynomialRegression::readSpecificMaps()
{
    m_polyOrder = m_intMap["polyOrder"];
    m_reg = m_doubleMap["initReg"];
    m_enableCrossInteractions = m_boolMap["enableCrossInteractions"];
}

void PolynomialRegression::modelInit()
{
    // add the tunable parameter
    paramDoubleValues.push_back ( &m_reg );
    paramDoubleNames.push_back ( "reg" );

    m_inputDim = m_polyOrder * m_nFeatures;

    if ( m_enableCrossInteractions )
    {
        m_inputDim = 0;

        // cross-interactions
        for ( int i=0;i<m_polyOrder*m_nFeatures;i++ )
            for ( int j=0;j<m_polyOrder*m_nFeatures+1;j++ )
                if ( j >= i )
                    m_inputDim++;
    }

    cout<<"Input dimension:"<<m_inputDim<<endl;

    // alloc mem for weights
    if ( m_x == 0 )
    {
        m_x = new REAL*[m_nCross+1];
        for ( int i=0;i<m_nCross+1;i++ )
            m_x[i] = new REAL[ ( m_inputDim + 1 ) * m_nClass * m_nDomain];

        // new mean/std
        m_polyMean = new REAL[m_inputDim];
        m_polyStd = new REAL[m_inputDim];
    }
}

void PolynomialRegression::predictAllOutputs ( REAL* rawInputs, REAL* outputs, uint nSamples, uint crossRun )
{
    REAL x,y;

    for ( int i=0;i<nSamples;i++ )
    {
        for ( int j=0;j<m_nClass*m_nDomain;j++ )
        {
            if ( m_enableCrossInteractions )
            {
                REAL sum = 0.0;
                REAL* inputPtr = rawInputs + i*m_nFeatures;
                REAL* xPtr = m_x[crossRun] + j;

                // cross-interactions
                int pos = 0;
                for ( int ii=0;ii<m_polyOrder*m_nFeatures;ii++ )
                    for ( int jj=0;jj<m_polyOrder*m_nFeatures+1;jj++ )
                        if ( jj >= ii )
                        {
                            int index0 = ii%m_nFeatures;
                            int index1 = jj%m_nFeatures;

                            int order0 = ( ii%m_polyOrder ) + 1;
                            int order1 = ( jj%m_polyOrder ) + 1;

                            x = inputPtr[index0];
                            y = inputPtr[index1];

                            x = power ( x,order0 );
                            y = power ( y,order1 );

                            if ( jj < m_polyOrder*m_nFeatures )
                                x = x * y;

                            x = ( x - m_polyMean[pos] ) / m_polyStd[pos];
                            sum += x * xPtr[pos*m_nClass*m_nDomain];
                            pos++;
                        }

                sum += xPtr[pos];
                if ( pos != m_inputDim )
                    assert ( false );

                outputs[i*m_nClass*m_nDomain + j] = sum;
            }
            else  // no feature interaction
            {
                REAL sum = 1.0 * m_x[crossRun][m_polyOrder*m_nFeatures*m_nClass*m_nDomain + j];
                for ( int order=0;order<m_polyOrder;order++ )
                {
                    REAL* inputPtr = rawInputs + i*m_nFeatures;
                    REAL* meanPtr = m_polyMean + order*m_nFeatures;
                    REAL* stdPtr = m_polyStd + order*m_nFeatures;
                    REAL* xPtr = m_x[crossRun] + order*m_nClass*m_nDomain*m_nFeatures + j;
                    for ( int k=0;k<m_nFeatures;k++ )
                    {
                        x = power ( inputPtr[k], order+1 );
                        x = ( x - meanPtr[k] ) / stdPtr[k];
                        sum += x * xPtr[k*m_nClass*m_nDomain];
                    }
                }
                outputs[i*m_nClass*m_nDomain + j] = sum;
            }
        }
    }
}

void PolynomialRegression::modelUpdate ( REAL* input, REAL* target, uint nSamples, uint crossRun )
{
    REAL x,y;
    REAL* crossTrain = 0;

    if ( m_enableCrossInteractions )
    {
        double minStd = 1e10, maxStd = -1e10, minMean = 1e10, maxMean = -1e10;

        // cross-interactions
        int pos = 0;
        for ( int ii=0;ii<m_polyOrder*m_nFeatures;ii++ )
            for ( int jj=0;jj<m_polyOrder*m_nFeatures+1;jj++ )
                if ( jj >= ii )
                {
                    int index0 = ii%m_nFeatures;
                    int index1 = jj%m_nFeatures;

                    int order0 = ( ii%m_polyOrder ) + 1;
                    int order1 = ( jj%m_polyOrder ) + 1;

                    // mean
                    double mean = 0.0, val;
                    for ( int j=0;j<nSamples;j++ )
                    {
                        x = input[j*m_nFeatures + index0];
                        y = input[j*m_nFeatures + index1];

                        x = power ( x,order0 );
                        y = power ( y,order1 );

                        if ( jj < m_polyOrder*m_nFeatures )
                            x = x * y;

                        mean += x;
                    }
                    mean /= ( double ) nSamples;

                    // standard deviation
                    double std = 0.0;
                    for ( int j=0;j<nSamples;j++ )
                    {
                        x = input[j*m_nFeatures + index0];
                        y = input[j*m_nFeatures + index1];

                        x = power ( x,order0 );
                        y = power ( y,order1 );

                        if ( jj < m_polyOrder*m_nFeatures )
                            x = x * y;

                        std += ( mean - x ) * ( mean - x );
                    }
                    std = sqrt ( std/ ( double ) ( nSamples-1 ) );
                    if ( std < m_standardDeviationMin )
                        std = m_standardDeviationMin;

                    minStd = minStd > std? std : minStd;
                    maxStd = maxStd < std? std : maxStd;
                    minMean = minMean > mean? mean : minMean;
                    maxMean = maxMean < mean? mean : maxMean;

                    // save them
                    m_polyMean[pos] = mean;
                    m_polyStd[pos] = std;
                    pos++;
                }

        //cout<<"Min|Max mean: "<<minMean<<"|"<<maxMean<<"   Min|Max std: "<<minStd<<"|"<<maxStd<<endl<<endl;


        if ( pos != m_inputDim )
            assert ( false );

        // apply to trainset
        crossTrain = new REAL[nSamples* ( m_inputDim+1 ) ];
        for ( int i=0;i<nSamples;i++ )
        {
            REAL* inputPtr = input + i*m_nFeatures;
            REAL* trainPtr = crossTrain + i* ( m_inputDim + 1 );

            // cross-interactions
            pos = 0;
            for ( int ii=0;ii<m_polyOrder*m_nFeatures;ii++ )
                for ( int jj=0;jj<m_polyOrder*m_nFeatures+1;jj++ )
                    if ( jj >= ii )
                    {
                        int index0 = ii%m_nFeatures;
                        int index1 = jj%m_nFeatures;

                        int order0 = ( ii%m_polyOrder ) + 1;
                        int order1 = ( jj%m_polyOrder ) + 1;

                        x = inputPtr[index0];
                        y = inputPtr[index1];

                        x = power ( x,order0 );
                        y = power ( y,order1 );

                        if ( jj < m_polyOrder*m_nFeatures )
                            x = x * y;

                        trainPtr[pos] = ( x - m_polyMean[pos] ) / m_polyStd[pos];
                        pos++;
                    }
            trainPtr[pos] = 1.0;

            if ( pos != m_inputDim )
                assert ( false );
        }

        // solve the linear system
        solver.RidgeRegressionMultisolutionSinglecall ( crossTrain, target, m_x[crossRun], nSamples, m_inputDim + 1, m_nClass*m_nDomain, m_reg, true );
    }
    else  // no feature interaction
    {
        crossTrain = new REAL[nSamples* ( m_polyOrder*m_nFeatures+1 ) ];

        //cout<<endl<<"Calculate new mean/std for all poly orders of input features x^(1.."<<m_polyOrder<<")"<<endl;
        double minStd = 1e10, maxStd = -1e10, minMean = 1e10, maxMean = -1e10;
        for ( int order=0;order<m_polyOrder;order++ )
        {
            for ( int i=0;i<m_nFeatures;i++ )
            {
                // mean
                double mean = 0.0, val;
                for ( int j=0;j<nSamples;j++ )
                {
                    val = input[j*m_nFeatures + i];
                    mean += power ( val, order+1 );
                }
                mean /= ( double ) nSamples;

                // standard deviation
                double std = 0.0;
                for ( int j=0;j<nSamples;j++ )
                {
                    val = input[j*m_nFeatures + i];
                    val = power ( val, order+1 );
                    std += ( mean - val ) * ( mean - val );
                }
                std = sqrt ( std/ ( double ) ( nSamples-1 ) );
                if ( std < m_standardDeviationMin )
                    std = m_standardDeviationMin;

                minStd = minStd > std? std : minStd;
                maxStd = maxStd < std? std : maxStd;
                minMean = minMean > mean? mean : minMean;
                maxMean = maxMean < mean? mean : maxMean;

                // save them
                m_polyMean[order*m_nFeatures + i] = mean;
                m_polyStd[order*m_nFeatures + i] = std;
            }
        }
        //cout<<"Min|Max mean: "<<minMean<<"|"<<maxMean<<"   Min|Max std: "<<minStd<<"|"<<maxStd<<endl<<endl;

        // copy train + add a constant input
        for ( int i=0;i<nSamples;i++ )
        {
            for ( int order=0;order<m_polyOrder;order++ )
            {
                int index = i* ( m_polyOrder*m_nFeatures + 1 ) + order * m_nFeatures;
                REAL* inputPtr = input + i*m_nFeatures;
                REAL* meanPtr = m_polyMean + order*m_nFeatures;
                REAL* stdPtr = m_polyStd + order*m_nFeatures;
                REAL* featurePtr = crossTrain + index;
                for ( int k=0;k<m_nFeatures;k++ )
                {
                    x = power ( inputPtr[k], order+1 );
                    x = ( x - meanPtr[k] ) / stdPtr[k];
                    featurePtr[k] = x;
                }
            }
            crossTrain[i* ( m_polyOrder*m_nFeatures + 1 ) + m_polyOrder*m_nFeatures] = 1.0;
        }

        // solve the linear system
        solver.RidgeRegressionMultisolutionSinglecall ( crossTrain, target, m_x[crossRun], nSamples, m_polyOrder * m_nFeatures + 1, m_nClass*m_nDomain, m_reg, true );
    }
    if ( crossTrain )
        delete[] crossTrain;
    crossTrain = 0;
}

REAL PolynomialRegression::power ( REAL x, int e )
{
    REAL tmp = x;
    for ( int i=1;i<e;i++ )
        tmp *= x;
    return tmp;
}

void PolynomialRegression::saveWeights ( int cross )
{
    char buf[1024];
    sprintf ( buf,"%02d",cross );
    string name = m_datasetPath + "/" + m_tempPath + "/" + m_weightFile + "." + buf;
    if ( m_inRetraining )
        cout<<"Save:"<<name<<endl;
    fstream f ( name.c_str(), ios::out );
    f.write ( ( char* ) &m_nTrain, sizeof ( int ) );
    f.write ( ( char* ) &m_nFeatures, sizeof ( int ) );
    f.write ( ( char* ) &m_nClass, sizeof ( int ) );
    f.write ( ( char* ) &m_nDomain, sizeof ( int ) );
    f.write ( ( char* ) &m_inputDim, sizeof ( int ) );
    f.write ( ( char* ) m_x[cross], sizeof ( REAL ) * ( m_inputDim+1 ) *m_nClass*m_nDomain );
    f.write ( ( char* ) m_polyMean, sizeof ( REAL ) *m_inputDim );
    f.write ( ( char* ) m_polyStd, sizeof ( REAL ) *m_inputDim );
    f.write ( ( char* ) &m_maxSwing, sizeof ( double ) );
    f.write ( ( char* ) &m_reg, sizeof ( double ) );
    f.close();
}

void PolynomialRegression::loadWeights ( int cross )
{
    char buf[1024];
    sprintf ( buf,"%02d",cross );
    string name = m_datasetPath + "/" + m_tempPath + "/" + m_weightFile + "." + buf;
    cout<<"Load:"<<name<<endl;
    fstream f ( name.c_str(), ios::in );
    if ( f.is_open() == false )
        assert ( false );
    f.read ( ( char* ) &m_nTrain, sizeof ( int ) );
    f.read ( ( char* ) &m_nFeatures, sizeof ( int ) );
    f.read ( ( char* ) &m_nClass, sizeof ( int ) );
    f.read ( ( char* ) &m_nDomain, sizeof ( int ) );
    f.read ( ( char* ) &m_inputDim, sizeof ( int ) );

    m_x = new REAL*[m_nCross+1];
    for ( int i=0;i<m_nCross+1;i++ )
        m_x[i] = 0;
    m_x[cross] = new REAL[ ( m_inputDim + 1 ) * m_nClass * m_nDomain];
    m_polyMean = new REAL[m_inputDim];
    m_polyStd = new REAL[m_inputDim];

    f.read ( ( char* ) m_x[cross], sizeof ( REAL ) * ( m_inputDim+1 ) *m_nClass*m_nDomain );
    f.read ( ( char* ) m_polyMean, sizeof ( REAL ) *m_inputDim );
    f.read ( ( char* ) m_polyStd, sizeof ( REAL ) *m_inputDim );
    f.read ( ( char* ) &m_maxSwing, sizeof ( double ) );
    f.read ( ( char* ) &m_reg, sizeof ( double ) );
    f.close();
}

void PolynomialRegression::loadMetaWeights ( int cross )
{
    char buf[1024];
    sprintf ( buf,"%02d",cross );
    string name = m_datasetPath + "/" + m_tempPath + "/" + m_weightFile + "." + buf;
    cout<<"LoadMeta:"<<name<<endl;
    fstream f ( name.c_str(), ios::in );
    if ( f.is_open() == false )
        assert ( false );
    f.read ( ( char* ) &m_nTrain, sizeof ( int ) );
    f.read ( ( char* ) &m_nFeatures, sizeof ( int ) );
    f.read ( ( char* ) &m_nClass, sizeof ( int ) );
    f.read ( ( char* ) &m_nDomain, sizeof ( int ) );
    f.read ( ( char* ) &m_inputDim, sizeof ( int ) );

    REAL* tmp = new REAL[ ( m_inputDim+1 ) *m_nClass*m_nDomain];
    f.read ( ( char* ) tmp, sizeof ( REAL ) * ( m_inputDim+1 ) *m_nClass*m_nDomain );
    delete[] tmp;
    tmp = new REAL[m_inputDim];
    f.read ( ( char* ) tmp, sizeof ( REAL ) *m_inputDim );
    f.read ( ( char* ) tmp, sizeof ( REAL ) *m_inputDim );
    delete[] tmp;
    f.read ( ( char* ) &m_maxSwing, sizeof ( double ) );
    f.read ( ( char* ) &m_reg, sizeof ( double ) );
    f.close();
}

string PolynomialRegression::templateGenerator ( int id, string preEffect, int nameID, bool blendStop )
{
    stringstream s;
    s<<"ALGORITHM=PolynomialRegression"<<endl;
    s<<"ID="<<id<<endl;
    s<<"TRAIN_ON_FULLPREDICTOR="<<preEffect<<endl;
    s<<"DISABLE=0"<<endl;
    s<<endl;
    s<<"[int]"<<endl;
    s<<"polyOrder=2"<<endl;
    s<<"maxTuninigEpochs=20"<<endl;
    s<<endl;
    s<<"[double]"<<endl;
    s<<"initMaxSwing=1.0"<<endl;
    s<<"initReg=1e-3"<<endl;
    s<<endl;
    s<<"[bool]"<<endl;
    s<<"enableCrossInteractions=0"<<endl;
    s<<"enableClipping=1"<<endl;
    s<<"enableTuneSwing=0"<<endl;
    s<<endl;
    s<<"minimzeProbe="<< ( !blendStop ) <<endl;
    s<<"minimzeProbeClassificationError=0"<<endl;
    s<<"minimzeBlend="<<blendStop<<endl;
    s<<"minimzeBlendClassificationError=0"<<endl;
    s<<endl;
    s<<"[string]"<<endl;
    s<<"weightFile=PolynomialRegression_"<<nameID<<"_weights.dat"<<endl;
    s<<"fullPrediction=PolynomialRegression_"<<nameID<<".dat"<<endl;

    return s.str();
}