|| || || || || || || ||

---

Fisher.h

Go to the documentation of this file.

#ifndef Impala_Core_Training_Fisher_h
#define Impala_Core_Training_Fisher_h

#include "Core/Training/Classifier.h"
#include "Core/Training/Utils.h"
#include "Core/Table/Sort.h"

#include "Core/Array/Set.h"
#include "Core/Array/SetVal.h"
#include "Core/Array/Array2dTem.h"
#include "Core/Array/Arrays.h"
#include "Core/Array/AddVal.h"
#include "Core/Array/MulVal.h"
#include "Core/Array/EqualsVal.h"
#include "Core/Array/PrintData.h"
#include "Core/Array/PixStat.h"

#include "Core/Matrix/MatTranspose.h"
#include "Core/Matrix/MatSVD.h"
#include "Core/Matrix/MatPInv.h"
#include "Core/Matrix/MatMul.h"
#include "Core/Matrix/MatCopy.h"
#include "Core/Matrix/MatReplicateMatrix.h"

#include "Core/Vector/VectorTem.h"
#include "Core/Vector/VectorSet.h"
#include "Core/Vector/Avg.h"
#include "Core/Vector/AddAssign.h"
#include "Core/Vector/MulAssign.h"

#include <iomanip>


namespace Impala
{
namespace Core
{
namespace Training
{

class Fisher : public Classifier
{
public:
    Fisher()
    {
        mTrainLabels=0;
        mTrainFeatures=0;
        mTestLabels=0;
        mTestFeatures=0;
        mRotation=0;
        mOffset=0;
        mMapped=0;
        mProbabilities=0;
        mCount=0;
    }

    ~Fisher()
    {
        delete mTrainLabels;
        delete mTrainFeatures;
        delete mTestLabels;
        delete mTestFeatures;
        delete mRotation;
        delete mOffset;
        delete mMapped;
        delete mProbabilities;
    }

    virtual void SetTrainSet(Feature::FeatureTable* data,
                             Table::AnnotationTable* annotation,
                             Table::QuidTable* mask)
    {
        ILOG_INFO("SetTrainSet");
        int featureLength = data->GetFeatureVectorLength();
        int sampleLength  = data->Size();
        ILOG_INFO(sampleLength<<" samples of "<<featureLength<<" features");
        
        if(mTrainFeatures)
            delete mTrainFeatures;
        if(mTrainLabels)
            delete mTrainLabels;

        mTrainLabels=0;
        mTrainFeatures=0;

        mTrainLabels = Array::ArrayCreate<Array::Array2dScalarReal64>
                    (1,sampleLength,0,0);
        mTrainFeatures = Array::ArrayCreate<Array::Array2dScalarReal64>
                    (featureLength,sampleLength,0,0);

        for(int i=0;i<sampleLength;i++)
        {
            const Impala::Real64* values = data->GetColumn2()->GetVectorData(i);
            for( int j=0;j<featureLength;j++)
                 mTrainFeatures->SetValue(values[j],j,i);

            mTrainLabels->SetValue(annotation->Get2(i),0,i);
            if((i==0)||(i==sampleLength-1)){
                ILOG_INFO("FeatureLength:"<<featureLength);
                ILOG_INFO(i<<".\t"<<mTrainFeatures->Value(0,i)<<" ... "
                                <<mTrainFeatures->Value(featureLength-1,i)<<" => "
                                <<mTrainLabels->Value(0,i));
            }
        }

        ILOG_INFO("FeatureLength:"<<featureLength);
        ILOG_INFO(0<<".\t"<<mTrainFeatures->Value(0,0)<<" ... "
                          <<mTrainFeatures->Value(featureLength-1,0)<<" => "
                          <<mTrainLabels->Value(0,0));
        ILOG_INFO("End of SetTrainSet");
        //CopyFeatureTable(data,annotation,&mTrainFeatures,&mTrainLabels);
    }

    virtual void SetTestSet(Feature::FeatureTable* data,
                            Table::AnnotationTable* annotation,
                            Table::QuidTable* mask)
    {
        ILOG_INFO("SetTestSet");
        int featureLength = data->GetFeatureVectorLength();
        int sampleLength  = data->Size();
        ILOG_INFO(sampleLength<<" samples of "<<featureLength<<" features");
        
        if(mTestFeatures)
            delete (mTestFeatures);
        if(mTestLabels)
            delete (mTestLabels);

        mTestLabels=0;
        mTestFeatures=0;

        mTestLabels = Array::ArrayCreate<Array::Array2dScalarReal64>
                    (1,sampleLength,0,0);
        mTestFeatures = Array::ArrayCreate<Array::Array2dScalarReal64>
                    (featureLength,sampleLength,0,0);

        for(int i=0;i<sampleLength;i++)
        {
            const Impala::Real64* values = data->GetColumn2()->GetVectorData(i);
            for( int j=0;j<featureLength;j++)
                mTestFeatures->SetValue(values[j],j,i);

            mTestLabels->SetValue(annotation->Get2(i),0,i);
            if((i==0)||(i==sampleLength-1)){
                ILOG_INFO(i<<".\t"<<mTestFeatures->Value(0,i)<<" ... "
                                <<mTestFeatures->Value(featureLength-1,i)<<" => "
                                <<mTestLabels->Value(0,i));
            }
            
        }
        ILOG_INFO("End of SetTestSet");
        //CopyFeatureTable(data,annotation,&mTestFeatures,&mTestLabels);
    }


    void SetTestSet(const Vector::VectorSet<Array::Array2dScalarReal64>* set)
    {
        int featureLength = set->GetVectorLength(0);
        int sampleLength  = set->Size();
        if(mTestFeatures)
            delete mTestFeatures;
        if(mTestLabels)
            delete mTestLabels;
        mTestLabels=0;
        mTestFeatures=0;

        mTestLabels = Array::ArrayCreate<Array::Array2dScalarReal64>
                    (1,sampleLength,0,0);
        mTestFeatures = Array::ArrayCreate<Array::Array2dScalarReal64>
                    (featureLength,sampleLength,0,0);
        SetVal(mTestLabels,0);

        for(int i=0;i<sampleLength;i++){
            const Impala::Real64* values = set->GetVectorData(i);
            for( int j=0;j<featureLength;j++)
                mTestFeatures->SetValue(values[j],i,j);
        }
        
    }

    void Train(Util::PropertySet* properties)
    { 
        int i,j;
        bool debug=false;
        int W = mTrainFeatures->CW();
        int H = mTrainFeatures->CH();

        Matrix::Mat* tmp    = 0;
        Matrix::Mat* tmp2   = 0;
        Matrix::Mat* tmp3   = 0;

        Matrix::Mat* Targets= 0;
        Matrix::Mat* FeatureAverages   = 0;
        Matrix::Mat* CenteredFeatureInverse = 0;
        Array::Array2dScalarReal64* CenteredFeatures   = 0;
        Array::Array2dScalarReal64* Features = 0;
        Array::Array2dScalarReal64* Labels   = 0;

        Vector::VectorSet<Array::Array2dScalarReal64> vs(true, W, H);

        Array::Set(Features,mTrainFeatures);
        Array::Set(Labels,mTrainLabels);

        //Clear the transformation if it was already calculated
        if(mRotation!=0)
            delete mRotation;
        if(mOffset!=0)
            delete mOffset;

        //Get The Average for each feature
        vs.SetStorage(Features);
        Vector::VectorTem<Real64> featureMeans=0;
        Avg(&featureMeans,&vs,H,0);
        
        //Center the data around 0 mean:
        Vector::MulAssign(featureMeans,Impala::Real64(-1));
        for(i=0;i<H;i++){
            Vector::VectorTem<Real64> tmpVec = vs.GetVector(i,true);
            AddAssign(tmpVec,featureMeans);
        }
        Vector::MulAssign(featureMeans,Impala::Real64(-1));

        //Now The data is centered around origin
        //We should add a column with 1's to the end of the feature vectors

        //If the rank of this Centered matrix with the ones column is 
        //less than the initial number of features for each sample,
        //Then we should find the Pseudo Inverse of this B matrice, and multiply
        //it with a two column labels matrix
        CenteredFeatures = Array::ArrayCreate<Array::Array2dScalarReal64>(W+1,H);
        SetVal(CenteredFeatures,1);
        SetPart(CenteredFeatures,Features,0,0,W,H,0,0);


        Targets         = Matrix::MatCreate<Array::Array2dScalarReal64>(H,2);
        FeatureAverages = Matrix::MatCreate<Array::Array2dScalarReal64>(1,W);

        for(int i=0;i<W;i++){
            FeatureAverages->SetValue(featureMeans.Elem(i),i,0);
        }

        SetPart(Targets,Labels,0,0,1,H,0,0);
        for(int i=0;i<H;i++){
            Targets->SetValue(Targets->Value(0,i),1,i);
            Targets->SetValue(-Targets->Value(0,i),0,i);
        }
        CenteredFeatureInverse=Matrix::MatPInv(CenteredFeatures);

            
        tmp = Matrix::MatMul(CenteredFeatureInverse,Targets);
        tmp2 = Matrix::MatCreate<Array::Array2dScalarReal64>(W,2);

        for(int i=0;i<W;i++){
            tmp2->SetValue(tmp->Value(0,i),0,i);
            tmp2->SetValue(tmp->Value(1,i),1,i);
        }
        tmp3 = Matrix::MatMul(FeatureAverages,tmp2);

        mOffset = Matrix::MatCreate<Array::Array2dScalarReal64>(1,2);
        mOffset->SetValue(tmp->Value(0,W)-tmp3->Value(0,0),0,0);
        mOffset->SetValue(tmp->Value(1,W)-tmp3->Value(1,0),1,0);

        mRotation = Matrix::MatCreate<Array::Array2dScalarReal64>(W,2);
        SetPart(mRotation,tmp,0,0,2,W,0,0);

        delete tmp;
        delete tmp2;
        delete tmp3;

        delete Targets;
        delete Labels;
        //delete Features; //will be released by the vector set
        delete FeatureAverages;
        delete CenteredFeatures;;
        delete CenteredFeatureInverse;
        CNormC();
        
    }
    // Classifier normalisation for ML posterior probabilities
    void ApplyMapping(Array::Array2dScalarReal64* feats)
    {
        if(mMapped)
            delete mMapped;
        if(mProbabilities)
            delete mProbabilities;


        int W=feats->CW();
        int H=feats->CH();

        Matrix::Mat* Features = Matrix::MatCreate<Array::Array2dScalarReal64>(H,W+1);

        Array::SetVal(Features,1);
        Array::SetPart(Features,feats,0,0,W,H,0,0);
        
        Matrix::Mat* Rot = Matrix::MatCreate<Array::Array2dScalarReal64>(W+1,2);
        Array::SetPart(Rot,mRotation,0,0,2,W,0,0);

        Rot->SetValue(mOffset->Value(0,0),0,W);
        Rot->SetValue(mOffset->Value(1,0),1,W);
        
        mMapped = Matrix::MatMul(Features,Rot);

        mProbabilities = SigmoidMap(mMapped,1);

        delete Rot;
        delete Features;
    
    }

    void CNormC()
    {
        int W=mTrainFeatures->CW();
        int H=mTrainFeatures->CH();

        ApplyMapping(mTrainFeatures);

        Impala::Real64 scale = 1E-10;
        Impala::Real64 reg = 1E-7;
        Impala::Real64 epsilon = 1E-6;
        Impala::Real64 likelihood= - std::numeric_limits<double>::infinity();
        Impala::Real64 likelihood_new= std::numeric_limits<double>::min() ;
        Impala::Real64 rmin= std::numeric_limits<double>::min() ;

        Impala::Real64 p;
        Array::Array2dScalarReal64* Positive = 0;
        Array::EqualsVal(Positive,mTrainLabels,1);
        p = Array::PixSum(Positive)/H;
        delete Positive;

        Matrix::Mat* pax;
        Matrix::Mat* pbx;
        Matrix::Mat* PosMapped = Matrix::MatCreate<Array::Array2dScalarReal64>(H,1);
        for(int i=0;i<H;i++)
            PosMapped->SetValue(mMapped->Value(mTrainLabels->Value(i,0)==1?1:0,i),0,i);

        Matrix::Mat* Probs = Matrix::MatCreate<Array::Array2dScalarReal64>(H,1);
        for(int i=0;i<H;i++)
            Probs->SetValue( ( mTrainLabels->Value(i,0)==1 ? p : 1-p ),0,i);

        while(abs(likelihood_new - likelihood) > epsilon)
        {
            Matrix::Mat* Mapped = 0;
            Array::MulVal(Mapped,PosMapped,scale);

            likelihood = likelihood_new;

            pax = SigmoidMap(Mapped,1);
            pbx = Matrix::MatCreate<Array::Array2dScalarReal64>(H,1);
            Array::MulVal(pbx,pax,-1);
            Array::AddVal(pbx,pbx,1);

            //Calculate new likelihood
            //
            Matrix::Mat* tmp = Matrix::MatCreate<Array::Array2dScalarReal64>(H,1);
            for(int i=0;i<H;i++)
                tmp->SetValue(Probs->Value(0,i)*log(pax->Value(0,i)+rmin),0,i);
            Impala::Real64 mean = PixSum(tmp)/H;

            likelihood_new = mean - reg*log(scale+1);

            //Calculate the new scale
            //
            for(int i=0;i<H;i++)
                tmp->SetValue(Probs->Value(0,i)*pbx->Value(0,i),0,i);
            
            Matrix::Mat* tmp2 = Matrix::MatCreate<Array::Array2dScalarReal64>(1,H);
            for(int i=0;i<H;i++)
                tmp2->SetValue(tmp->Value(0,i),i,0);

            Matrix::Mat* tmp3=Matrix::MatMul(tmp2,PosMapped);

            Matrix::Mat* tmp4 = Matrix::MatCreate<Array::Array2dScalarReal64>(1,H);
            Matrix::Mat* tmp5 = Matrix::MatCreate<Array::Array2dScalarReal64>(H,1);

            for(int i=0;i<H;i++)
                tmp4->SetValue(Probs->Value(0,i)*PosMapped->Value(0,i)*pax->Value(0,i),i,0);

            for(int i=0;i<H;i++)
                tmp5->SetValue(PosMapped->Value(0,i)*pbx->Value(0,i),0,i);

            Matrix::Mat* tmp6 = Matrix::MatMul(tmp4,tmp5);

            Impala::Real64 a = tmp3->Value(0,0);
            Impala::Real64 b = reg*H/(scale+1);
            Impala::Real64 c = tmp6->Value(0,0) - H* reg/((scale+1)*(scale+1)) + rmin;

            scale = scale + (a - b) / c;

            delete Mapped;
            delete pax;
            delete pbx;
            delete tmp;
            delete tmp2;
            delete tmp3;
            delete tmp4;
            delete tmp5;
            delete tmp6;
            pax=0;
            pbx=0;
        }

        Impala::Real64 var;
        Impala::Real64* o=0;
        
        Array::PixStat(PosMapped,o,o,o,&var);
        
        // Dirty trick to avoid 0 and inf.
        scale = std::min(std::max(scale,0.1/(sqrt(var)+1e-16)),1e16);

        Array::MulVal(mRotation,mRotation,scale);
        Array::MulVal(mOffset,mOffset,scale);

        delete Probs;
        delete PosMapped;
    }

    Matrix::Mat*
    SigmoidMap(Matrix::Mat* Map,Impala::Real64 scale=1)
    {
        Matrix::Mat* res = 0;
        res = Matrix::MatCopy(Map);
        for(int i=0;i<Map->CW();i++)
            for(int j=0;j<Map->CH();j++)
                res->SetValue(1/(1+exp(-res->Value(i,j)/scale)),i,j);

        return res;
    }

    virtual double Predict(const Vector::VectorTem<double>* feature)
    {
        
        Matrix::Mat* tmp = Matrix::MatCreate<Array::Array2dScalarReal64>
                                            (1,feature->Size(),
                                            const_cast<Impala::Real64*>(feature->GetData())
                                            ,true);
        ApplyMapping(tmp);
        delete tmp;
        return mProbabilities->Value(1,0);
    }

    virtual RankingTableType* Rank()
    {
        if(mTestFeatures == 0)
        {
            ILOG_ERROR("[Fisher::Train] no test set");
            return NULL;
        }
        int size=mTestFeatures->CH();
        RankingTableType* result = MakeRankingTable(size);

        ApplyMapping(mTestFeatures);

        for(int i=0 ; i<size ; i++)
        {
            result->Add("", i, mTestLabels->Value(i,0), mProbabilities->Value(1,i));
        }

        return result;
    }

    virtual void
    Predict(Table::SimilarityTableSet::SimTableType* result)
    {

        if (mTestFeatures == 0)
        {
            ILOG_ERROR("[Fisher::Predict] no test set");
            return;
        }
        result->SetSize(0);
        ApplyMapping(mTestFeatures);
        int size=mTestFeatures->CH();
        
        for (int i=0 ; i<size ; i++)
        {
            result->Add(mProbabilities->Value(1,i));
        }
    }

    void LoadModel(const std::string& name, Util::Database* db)
    {
        if (Util::Database::GetInstance().GetDataChannel())
        {
            std::string tmpName = FileNameTmp();
            Util::FileCopyRemoteToLocal(name, tmpName);
            
            Util::Database& db = Util::Database::GetInstance();
            Util::IOBuffer* buf = db.GetIOBuffer(tmpName, true, false, "");
            ReadRaw(mRotation,buf);
            ReadRaw(mOffset,buf);
            delete buf;

            unlink(tmpName.c_str());
        }
        else
        {
            Util::Database& db = Util::Database::GetInstance();
            Util::IOBuffer* buf = db.GetIOBuffer(name, true, false, "");
            ReadRaw(mRotation,buf);
            ReadRaw(mOffset,buf);
            delete buf;
        }
    }
    void SaveModel(const std::string& name, Util::Database* db)
    {
        if (Util::Database::GetInstance().GetDataChannel())
        {
            std::string tmpName = FileNameTmp();
            //if(svm_save_model(name.c_str(), mModel) == -1)
            Util::Database& db = Util::Database::GetInstance();
            Util::IOBuffer* buf = db.GetIOBuffer(tmpName, true, false, "");

            WriteRaw(mRotation,buf,true);
            WriteRaw(mOffset,buf,true);
            delete buf;

            Util::FileCopyLocalToRemote(tmpName, name);
            unlink(tmpName.c_str());
        }
        else
        {
            Util::Database& db = Util::Database::GetInstance();
            Util::IOBuffer* buf = db.GetIOBuffer(name, true, false, "");

            WriteRaw(mRotation,buf,true);
            WriteRaw(mOffset,buf,true);
            delete buf;
        }
    }
private:
    Array::Array2dScalarReal64* mTrainFeatures;
    Array::Array2dScalarReal64* mTrainLabels;
    Array::Array2dScalarReal64* mTestFeatures;
    Array::Array2dScalarReal64* mTestLabels;

    Matrix::Mat* mRotation;
    Matrix::Mat* mOffset;
    Matrix::Mat* mMapped;
    Matrix::Mat* mProbabilities;

    int mSampleLength;
    int mFeatureLength;
    ILOG_VAR_DEC;
    int mCount;
};

ILOG_VAR_INIT(Fisher,Impala.Core.Training);

}//namespace Core
}//namespace Training
}//Impala


#endif

---