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

---

MarkovStationaryFeature.h

Go to the documentation of this file.

#ifndef Impala_Core_Feature_MarkovStationaryFeature_h
#define Impala_Core_Feature_MarkovStationaryFeature_h

#include "Core/Array/Add.h"
#include "Core/Array/MulVal.h"
#include "Core/Vector/VectorTem.h"
#include "Core/Vector/Sum.h"
#include "Core/Matrix/Mat.h"
#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/MatSum.h"
#include "Core/Matrix/MatCopy.h"
#include "Core/Matrix/MatMul.h"
#include "Core/Matrix/GetDiagonal.h"

namespace Impala
{
namespace Core
{
namespace Feature
{


void
Dump(Matrix::Mat* matrix)
{
    using namespace Matrix;
    for(int i = 0; i < MatNrRow(matrix); i++)
    {
        for(int j = 0; j < MatNrCol(matrix); j++)
        {
            std::cout << *MatE(matrix, i, j) << " ";
        }
        std::cout << std::endl;
    }
}


inline Vector::VectorTem<Real64>
MarkovStationaryFeature(Matrix::Mat* cooccuranceMatrix)
{
    using namespace Matrix;
    const int ITERATION_COUNT = 100;

    int n = MatNrRow(cooccuranceMatrix);
    if(MatNrCol(cooccuranceMatrix) != n)
    {
        // it cannot be!
        std::cerr << "cooccuranceMatrix must be square" << std::endl;
        throw "bye";
    }

    /* descriptor storage */
    Vector::VectorTem<Real64> descriptor(2 * n);

    /* normalize the co-ocurrance matrix */
    Mat* p = MatCreate<Mat>(n, n);
    Mat* normalizationFactors = MatSumAxis1(cooccuranceMatrix);
    for(int vj = 0; vj < n; vj++)
    {
        if(*MatE(normalizationFactors, vj, 0) == 0.0)
        {
            *MatE(normalizationFactors, vj, 0) = 1.0;
        }
        for(int i = 0; i < n; i++) 
        {
            *MatE(p, vj, i) = *MatE(cooccuranceMatrix, vj, i) / *MatE(normalizationFactors, vj, 0);
        }
    }
    //Dump(cooccuranceMatrix);
    //Dump(normalizationFactors);
    //Dump(p);
    delete normalizationFactors;

    /* we now have a transition matrix in p */
    Mat* identity = MatCreate<Mat>(n, n);
    SetVal(identity, 0.0);
    // set to identity matrix
    for(int i = 0; i < n; i++) 
    {
        *MatE(identity, i, i) = 1.0;
    }
    Mat* a = MatCopy(identity);
    Mat* previousPower = identity;
    // identity will be deleted in the loop as previousPower
    for(int k = 0; k < ITERATION_COUNT; k++)
    {
        Mat* currentPower = MatMul(previousPower, p);
        // perform matrix addition
        Add(a, a, currentPower);
        delete previousPower;
        previousPower = currentPower;
    }
    //std::cout << "previous:" << std::endl; Dump(previousPower);
    delete previousPower;
    MulVal(a, a, 1.0 / (ITERATION_COUNT + 1));
    //std::cout << "a_n:" << std::endl; Dump(a);
    
    /* construct the stationary disribution by taking the averages of the
       rows of a (all rows of a should be equal upon convergence) */
    Mat* stationarySum = MatSumAxis0(a);
    //Dump(stationarySum);
    for(int i = 0; i < n; i++)
    {
        descriptor[n + i] = *MatE(stationarySum, 0, i) / n;
    }
    delete a;
    delete stationarySum;

    /* construct initial distribution by normalizing the self-transition */
    Vector::VectorTem<Real64> diagonal = GetDiagonal(cooccuranceMatrix);
    Real64 normalizationFactor = Sum(diagonal);
    if(normalizationFactor == 0.0) normalizationFactor = 1.0;
    Vector::VectorTem<Real64> initialDistribution = diagonal * (1.0 / normalizationFactor);
    
    /* combine initial distribution and stationary distribution into a single
       vector */    
    for(int i = 0; i < n; i++)
    {
        descriptor[i] = initialDistribution[i];
    }
    delete p;
    
    return descriptor;
}


inline Vector::VectorTem<Real64>
MarkovStationaryFeatureTimed(Matrix::Mat* m)
{
    ILOG_VAR(Sandbox.koen.MarkovStationaryFeature);
    Timer timer;
    Vector::VectorTem<Real64> temp = MarkovStationaryFeature(m);
    ILOG_INFO("MSF in " << timer.SplitTimeStr());
    return temp;
}


void
testMSF()
{
    ILOG_VAR(Impala.Core.Feature.testMSF);
    using namespace Matrix;
    Mat* a = 0;
    a = MatCreate<Mat>(8, 8);
    
            *MatE(a, 0, 0) = 1;
            *MatE(a, 0, 1) = 2;
            *MatE(a, 0, 2) = 3;
            *MatE(a, 0, 3) = 4;
            *MatE(a, 0, 4) = 5;
            *MatE(a, 1, 0) = 4;
            *MatE(a, 1, 1) = 3;
            *MatE(a, 1, 2) = 2;
            *MatE(a, 1, 3) = 9;
            *MatE(a, 1, 4) = 0;

    a = MatCreate<Mat>(3, 3);
            *MatE(a, 0, 0) = 160;
            *MatE(a, 0, 1) = 30;
            *MatE(a, 0, 2) = 10;
            *MatE(a, 1, 0) = 30;
            *MatE(a, 1, 1) = 340;
            *MatE(a, 1, 2) = 30;
            *MatE(a, 2, 0) = 10;
            *MatE(a, 2, 1) = 30;
            *MatE(a, 2, 2) = 160;

    a = MatCreate<Mat>(3, 3);
            *MatE(a, 0, 0) = 160;
            *MatE(a, 0, 1) = 40;
            *MatE(a, 0, 2) = 0;
            *MatE(a, 1, 0) = 40;
            *MatE(a, 1, 1) = 340;
            *MatE(a, 1, 2) = 40;
            *MatE(a, 2, 0) = 0;
            *MatE(a, 2, 1) = 40;
            *MatE(a, 2, 2) = 160;

    ILOG_INFO("MSF start");
    Vector::VectorTem<Real64> descriptor = MarkovStationaryFeature(a);
    for(int i = 0; i < descriptor.Size(); i++)
    {
        std::cout << descriptor[i] << " ";
    }
    std::cout << std::endl;
    ILOG_INFO("MSF stop");
}

} // namespace Koen
} // namespace Sandbox
} // namespace Impala

#endif

---