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

---

mainPrecomputeKernelMatrix.cpp

Go to the documentation of this file.

#include "Basis/CmdOptions.h"
#include "Link/Mpi/MpiFuncs.h"

#include "Core/Matrix/DistributedAccess.h"
#include "Core/Array/Pattern/PatMPixOp.h"
#include "Core/Array/PixSum.h"
#include "Core/Array/Set.h"
#include "Core/Array/Add.h"
#include "Core/Array/DivVal.h"
#include "Core/Array/MulVal.h"
#include "Core/Array/Exp.h"
#include "Core/Array/WriteRaw.h"
#include "Core/Array/MakeFromValue.h"
#include "Core/Array/PrintData.h"
#include "Core/Feature/FeatureTable.h"
#include "Core/Feature/FeatureDefinition.h"
#include "Core/Feature/FeatureTable.h"
#include "Core/Table/Write.h"
#include "Core/Table/Equals.h"
#include "Core/Table/Select.h"
#include "Core/Vector/Apply.h"
#include "Core/VideoSet/MakeVideoSet.h"
#include "Core/ImageSet/MakeImageSet.h"
#include "Core/Training/KernelFunctions.h"
#include "Util/PropertySet.h"
#include "Core/Database/MakeRawDataSet.h"

namespace Impala
{
namespace Application
{

using namespace Core;
using Core::Database::RawDataSet;

bool CheckParameteres(CmdOptions& options,
                      RawDataSet* dataset, RawDataSet* dataset2,
                      std::vector<Feature::FeatureDefinition>& featureDefs,
                      std::vector<double>& weights,
                      String resultname)
{
    /*
    The inputFeatures can no longer be given through the arguments only through
    an option.
    */
    ILOG_VAR(Application.mainPrecomputeKernelMatrix.CheckParameters);

    /* add input features through the INI file */
    String inputFeatures = options.GetString("inputFeatures");
    if(!inputFeatures.empty())
    {
        ILOG_DEBUG("imput features: "<< inputFeatures);
        Util::StringParser p(inputFeatures);
        while(!p.TheEnd())
        {
            double weight = p.GetDouble();
            ILOG_DEBUG("w: "<<weight);
            if(weight == 0.0)
            {
                ILOG_ERROR_ONCE("\"" << inputFeatures << "\" contains a zero weight");
                return false;
            }
            String name = p.GetString(' ', false);
            if(name.empty())
            {
                ILOG_ERROR_ONCE("While parsing --inputFeatures: unbalanced length, empty name");
                return false;
            }
            ILOG_DEBUG("n: "<< name);
            weights.push_back(weight);
            Feature::FeatureDefinition fdef(name);
            name = dataset->GetFilePathFeatureIndex(fdef, "", false, false);
            if(name.empty())
            {
                ILOG_ERROR_ONCE("could not open \"" << fdef.AsString() << "\"");
                return false;
            }
            featureDefs.push_back(fdef);
        }
    }
    
    if(Link::Mpi::MyId() == 0)
    {
        String filename;
        if (dataset2)
            filename = dataset2->GetFilePathPrecomputedKernels
                (resultname+".input.txt", dataset->GetSetName(),
                 true, false);
        else
            filename = dataset->GetFilePathPrecomputedKernels
                (resultname+".input.txt", "", true, false);
            
        Util::Database* db = (dataset2) ? dataset2->GetDatabase()
                                        : dataset->GetDatabase();
        Util::IOBuffer* buf = db->GetIOBuffer(filename, false, false, "tmp");
        if(buf)
        {
            for(int i=0 ; i<weights.size() ; ++i)
            {
                String s = MakeString(weights[i]) + " " + featureDefs[i].AsString();
                buf->Puts(s);
            }
            delete buf;
        }
        else
            ILOG_ERROR_ONCE("could not open input.txt");
    }

    return true;
}


Feature::FeatureTable* OpenFeatureTable(Feature::FeatureDefinition& featureDef,
                                        RawDataSet* dataset)
{
    /* if MPI
    node 0 opens IOBuffer and Broadcast
    every node loads table from IOBuffer
    every node only selects it's own part of the table
    */
    ILOG_VAR(Application.mainPrecomputeKernelMatrix.OpenFeatureTable);
    Feature::FeatureTable* f = 0;
    std::string filename = 
        dataset->GetFilePathFeatureIndex(featureDef, "", false, false);
    ILOG_INFO_ONCE("opening " << featureDef.AsString());
#ifndef MPI_USED
    f = Feature::FeatureTable::MakeFromFile(featureDef, filename,
                                            dataset->GetDatabase());
    //f->SetSize(20);
    //ILOG_WARNING("***debug statment should be removed!***");
#else

    Util::IOBuffer* buf;
    int myId = Link::Mpi::MyId();
    if(myId == 0)
    {
        // use the dataset-relative database (!)
        Util::Database* db = dataset->GetDatabase();
        buf = db->GetIOBuffer(filename, true, true, "");
        if(buf == 0 || !buf->Valid())
        {
            ILOG_ERROR("node0 couldn't open file");
        }
    }
    else
    {
        buf = new Util::IOBuffer();
    }
    Broadcast(buf);
    if (buf && buf->Valid())
    {
        f = new Feature::FeatureTable(featureDef);
        Read(f, buf);
        delete buf;
    }
    else
        ILOG_ERROR(myId << ": could not read file from buffer");
    Link::Mpi::Barrier();
#endif
    ILOG_INFO_ONCE("size = " << f->Size() << "; length = " << f->GetFeatureVectorLength());

    return f;
}

void GetPartialTask(int& partcount, int& row, int& column)
{
    ILOG_VAR(Application.mainPrecomputeKernelMatrix.GetPartialTask);
    if(Link::Mpi::NrProcs() == 1)
    {
        partcount = 1;
        row = 0;
        column = 0;
        return;
    }

    int cpuCount = Link::Mpi::NrProcs();
    partcount = sqrt((double)cpuCount);
    int myId = Link::Mpi::MyId();
    column = myId % partcount;
    row = myId / partcount;
    if(row >= partcount)
    {
        partcount = 0;
        row = 0;
        column = 0;
    }
}

void CheckQuids(Feature::FeatureTable* f, RawDataSet* set, RawDataSet* set2,
                String resultname, int part, int partcount)
{
    ILOG_VAR(Application.mainPrecomputeKernelMatrix.CheckQuids);
    bool equal = true;
    bool fileExists = true;

    String basename = resultname + MakeString(part) + "of" + 
                      MakeString(partcount) + ".quids";
    String filename;
    if (set2)
    {
        // try to open for reading
        filename = set2->GetFilePathPrecomputedKernels
            (basename, set->GetSetName(), false, true);
        if (filename.empty())
        {
            fileExists = false;
            filename = set2->GetFilePathPrecomputedKernels
                (basename, set->GetSetName(), true, false);
        }
    }
    else
    {
        // try to open for reading
        filename = set->GetFilePathPrecomputedKernels(basename, "", false,
                                                      true);
        if (filename.empty())
        {
            fileExists = false;
            filename = set->GetFilePathPrecomputedKernels(basename, "", true,
                                                          false);
        }
    }
    Table::QuidTable* quids = f->GetQuidTable();
    Table::QuidTable* prevquids = new Table::QuidTable;
    if (fileExists)
    {
        Read(prevquids, filename, set->GetDatabase());
        equal = Equals(quids, prevquids);
    }
    else //we could not read so we try to write
    {
        if (!Write(quids, filename, set->GetDatabase(), true))
            ILOG_ERROR("could neither read nor write " << filename);
    }
    delete prevquids;
    delete quids;
    if(!equal)
    {
        ILOG_ERROR("the input tables do not have the same Quids");
        exit(0);
    }
}



Feature::FeatureTable* GetPartial(Feature::FeatureTable* f,
                                  int partnumber, int partcount)
{
    ILOG_VAR(Application.mainPrecomputeKernelMatrix.GetPartial);
    Feature::FeatureTable* part = 
        new Feature::FeatureTable(f->GetFeatureDefinition(),0,
                                  f->GetFeatureVectorLength());
    int s = f->Size();
    Matrix::DistributedAccess::IndexConverter indexConvert(f->Size(), partcount);
    int from = indexConvert.PartToIndex(partnumber);
    int to   = indexConvert.PartToIndex(partnumber + 1);
    //int from = ceil((s*partnumber)/(double)partcount);
    //int to = ceil((s*(partnumber+1))/(double)partcount);
    Table::Select(part, f, from, to, true);
    return part;
}


Matrix::Mat* ComputeMatrix(Feature::FeatureTable *devel,
                           Feature::FeatureTable *test, String resultname,
                           RawDataSet* set, RawDataSet* set2)
{
    ILOG_VAR(Application.mainPrecomputeKernelMatrix.ComputeMatrix);
    int partcount, row, column;
    GetPartialTask(partcount, row, column);
    if(partcount == 0)
        ILOG_ERROR_NODE("partcount == 0: shouldn't reach this point in code");
    Feature::FeatureTable* horizontal;
    Feature::FeatureTable* vertical;
    if(partcount == 1)
    {
        horizontal = devel;
        vertical = test;
    }
    else
    {
        horizontal = GetPartial(devel, column, partcount);
        vertical = GetPartial(test, row, partcount);
        delete devel;
        if(test != devel)
            delete test;
    }
    if(column == 0)
        CheckQuids(vertical, set, set2, resultname + ".rowindices", row, partcount);
    if(row == 0)
        CheckQuids(horizontal, set, set2, resultname + ".columnindices", column, partcount);


    // fill a matrix with chi2 distance
    typedef Feature::FeatureTable::ColumnVectorSet VectorSet;
    VectorSet* vectorsH = horizontal->GetColumn2();
    VectorSet* vectorsV = vertical->GetColumn2();
    vectorsH->SetSize(horizontal->Size());
    vectorsV->SetSize(vertical->Size());
    Matrix::Mat* distanceMatrix = Apply(&Training::Chi2Distance, vectorsH, vectorsV);
    delete horizontal;
    if(horizontal != vertical)
        delete vertical;
    return distanceMatrix;
}


double GetAverage(Matrix::Mat *distanceMatrix)
{
    double average;
    average = Array::PixSum(distanceMatrix);
    average = Link::Mpi::AllReduceSum(average);
    int pixcount = distanceMatrix->W() * distanceMatrix->H();
    pixcount = Link::Mpi::AllReduceSum(pixcount);
    average /= pixcount;
    return average;
}

void
WriteInfoFile(int columns, int rows, int partcount, String filepathname,
              Util::Database* db)
{
    if(Link::Mpi::MyId() != 0)
        return;
    ILOG_VAR(Application.mainPrecomputeKernelMatrix.CheckQuids);

    Util::IOBuffer* buf = db->GetIOBuffer(filepathname, false, false, "tmp");
    Util::PropertySet ps;
    ps.Add("totalrows", rows);
    ps.Add("totalcolumns", columns);
    ps.Add("rowparts", partcount);
    ps.Add("columnparts", partcount);
    ps.Print(buf);
    delete buf;
}

void
WriteAverages(String filepathname, Util::Database* db,
              std::vector<double> averages)
{
    if(Link::Mpi::MyId() == 0)
    {
        Util::IOBuffer* buf = db->GetIOBuffer(filepathname, false, false, "tmp");
        if (buf)
        {
            for(int i=0 ; i<averages.size() ; ++i)
                buf->Write(&averages[i], sizeof(double));
        }
        delete buf;
    }
}

void
LoadAverages(RawDataSet* set, String filepathname, std::vector<double>& averages)
{
        Util::Database* db = set->GetDatabase();
        Util::IOBuffer* buf = db->GetIOBuffer(filepathname, true, false, "");
        if (buf)
        {
            for(int i=0 ; i<averages.size() ; ++i)
                buf->Read(&averages[i], sizeof(double));
        }
        delete buf;
}

// resultname already has the path prepended to it
void
WriteResult(String resultname, Util::Database* db, Matrix::Mat *accumulator)
{
    int partcount, row, column;
    GetPartialTask(partcount, row, column);
    std::ostringstream filename;
    filename << resultname << ".precomputed";
    filename << ".part-R" << row << "-C" << column;
    filename << ".raw";
    WriteRaw(accumulator, filename.str(), db, true);
}

int mainPrecomputeKernelMatrix(CmdOptions& options)
{
    // this name is 26 characters. if you are tracing unexplainable C++ string
    // exceptions, check your log4cpp.properties to see if you are limiting
    // this name to just 25 characters...
    ILOG_VAR(Application.mainPrecomputeKernelMatrix);
    RawDataSet* dataset = Core::Database::MakeRawDataSet(options.GetArg(0));
    if(!dataset)
    {
        ILOG_INFO_ONCE("failed to open dataset " << options.GetArg(0));
        return 1;
    }
    String name2 = options.GetArg(1);
    RawDataSet* dataset2 = 0;
    if((name2 != "0") && (atof(name2) == 0.0))    // it is not a number, or 0
    {
        dataset2 = Core::Database::MakeRawDataSet(name2, true);
    }
    if(dataset2)
        ILOG_INFO_ONCE("2 sets loaded: " << dataset2->GetSetName());

    // for distribution we compute what part of the matrix we must compute
    // in single process case partcount == 1
    // in multi process case it is sqrt(NrProcs()) (or 0 if this node doesn't compute anything)
    int partcount, row, column;
    GetPartialTask(partcount, row, column);
    ILOG_INFO_ONCE("Using a " << partcount << "x" << partcount << "grid");
    int cpuCount = Link::Mpi::NrProcs();
    ILOG_INFO_ONCE("total nodes = " << cpuCount << " unused nodes = " 
                   << cpuCount - partcount*partcount);
    if(partcount == 0)
    {
        ILOG_ERROR("unsupported number of nodes, number of nodes MUST be"
                   << " sqare of a natural number");
        exit(0);
    }
    else
    {
        ILOG_DEBUG_NODE("col = " << column << " row = " << row);
    }

    std::vector<Feature::FeatureDefinition> featureDefs;
    std::vector<double> weights;
    String resultname = options.GetArg(options.GetNrArg()-1);
    if(!CheckParameteres(options, dataset, dataset2, featureDefs, weights, 
                         resultname))
        return 0;
    std::vector<double> averages;
    averages.resize(weights.size());
    if(dataset2)
    {
        String filename = resultname + ".averages.raw";
        filename = dataset->GetFilePathPrecomputedKernels(filename, "", false, false);
        if(filename == "")
        {
            ILOG_ERROR("couldn't read averages");
            exit(0);
        }
        LoadAverages(dataset, filename, averages);
        for (int i = 0; i < averages.size(); i++)
        {
            ILOG_INFO_ONCE("average " << i << " " << averages[i]);
        }
    }

    if(!dataset2)
    {
        String tmp = dataset->GetFilePathPrecomputedKernels(resultname+".info",
                       "", true, true);
        String tmp2 = dataset->GetFilePathPrecomputedKernels(resultname+".averages.raw",
                       "", true, true);
        if(tmp.empty() && tmp2.empty())
        {
            ILOG_INFO("A complete kernel already exists, nothing to do.");
            return 0;
        }
        if(tmp.empty())
        {
            ILOG_ERROR("Incomplete kernel exists! Cleanup files first.");
            return 1;
        }
    }

    Util::Database* db = (dataset2) ? dataset2->GetDatabase()
                                    : dataset->GetDatabase();
    Matrix::Mat* accumulator = 0;
    double totalweight = 0;
    for(int i=0 ; i<weights.size() ; ++i)
    {
        // open feature tables
        double weight = weights[i];
        Feature::FeatureTable* f1 = OpenFeatureTable(featureDefs[i], dataset);
        Feature::FeatureTable* f2 = f1;
        if(dataset2)
            f2 = OpenFeatureTable(featureDefs[i], dataset2);
        if(Link::Mpi::MyId() == 0 && i == 0)
        {
            // create a .info file once
            String filename = resultname+".info";
            if(dataset2)
                filename = dataset2->GetFilePathPrecomputedKernels(filename,
                    dataset->GetSetName(), true, false);
            else
                filename = dataset->GetFilePathPrecomputedKernels(filename,
                    "", true, false);
            if(filename == "")
            {
                ILOG_ERROR("could not save .info file");
                exit(0);
            }
            ILOG_INFO("Saving info in " << filename);
            WriteInfoFile(f1->Size(), f2->Size(), partcount, filename, db);

            // write the total list of devel-quids once
            if (dataset2)
                filename = dataset2->GetFilePathPrecomputedKernels(
                          resultname + ".columns.tab", dataset->GetSetName(),
                          true, false);
            else
                filename = dataset->GetFilePathPrecomputedKernels
                    (resultname + ".columns.tab", "", true, false);
            ILOG_INFO("Saving columns in " << filename);
            if (!filename.empty())
                Write(f1->GetQuidTable(), filename, dataset->GetDatabase(),
                      true);
            else
                ILOG_ERROR("Unable to write columns.tab");

        }
        // compute kernel distances between features
        Matrix::Mat *distanceMatrix = ComputeMatrix(f1, f2, resultname, 
                                                    dataset, dataset2);
        // f1 & f2 are deleted by ComputeMatrix
        
        // if this is the 'learn' kernel matrix we compute the averages and
        // communicate between nodes to make sure everyone has access to them
        // otherwise, the averages are already loaded into the vector
        if(dataset2 == 0)
        {
            double average = GetAverage(distanceMatrix);
            averages[i] = average;
        }

        // accumulate; this is the part inside the 'exp' in the kernel function
        MulVal(distanceMatrix, distanceMatrix, weight);
        DivVal(distanceMatrix, distanceMatrix, -averages[i]);
        if(accumulator == 0)
            Set(accumulator, distanceMatrix);
        else
            Add(accumulator, accumulator, distanceMatrix);
        delete distanceMatrix;

        totalweight += weight;
    }

    ILOG_INFO_ONCE("finalising...");
    // only thing to do is weigh and exp
    DivVal(accumulator, accumulator, totalweight);
    Exp(accumulator, accumulator);

    if(dataset2 == 0)
    {
        if(Link::Mpi::MyId() == 0)
        {
            String filename = resultname + ".averages.raw";
            filename = dataset->GetFilePathPrecomputedKernels(filename, "", true, false);
            WriteAverages(filename, db, averages);
        }
    }
    // write the matrix
    if (dataset2)
    {
        resultname = dataset2->GetFilePathPrecomputedKernels
            (resultname, dataset->GetSetName(), true, false);
    }
    else
    {
        resultname = dataset->GetFilePathPrecomputedKernels
            (resultname, "", true, false);
    }
    WriteResult(resultname, db, accumulator);
    return 0;
}

Matrix::Mat* CreateTestMat(int row, int column, int partCount, int totalSize)
{
    int xbegin = (totalSize*column)/partCount;
    int xend = (totalSize*(column+1))/partCount;
    int ybegin = (totalSize*row)/partCount;
    int yend = (totalSize*(row+1))/partCount;
    Matrix::Mat* m = new Matrix::Mat(xend - xbegin, yend - ybegin, 0, 0);
    for(int y=ybegin ; y<yend ; ++y)
    {
        for(int x=xbegin ; x<xend ; ++x)
            *(m->CPB(x-xbegin,y-ybegin)) = (x+1.) / (y+1.);
    }
    return m;
                           
}

void
CreateTestQuids(String resultname, Util::Database* db, int part, int partCount,
                int totalSize)
{
    ILOG_VAR(Application.CreateTestQuids);
    String filename = resultname + MakeString(part) + "of" + MakeString(partCount) + ".quids";
    int begin = (totalSize*part)/partCount;
    int end = (totalSize*(part+1))/partCount;
    Table::QuidTable* quids = new Table::QuidTable(end-begin);
    for(int i=begin ; i<end ; ++i)
        quids->Add(i);
    // fill table with quids
    if (!Write(quids, filename, db, false))
        ILOG_ERROR("could not write " << filename);
    delete quids;
}

int makeTestMatrix(CmdOptions& options)
{
    ILOG_VAR(Application.makeTestMatrix);
    int partCount = 3;
    int totalSize = 11;
    String resultname = options.GetArg(options.GetNrArg()-1);
    Util::Database* db = &Util::Database::GetInstance();
    WriteInfoFile(totalSize, totalSize, partCount, resultname, db);

    for(int row=0 ; row<partCount ; ++row)
        for(int column=0 ; column<partCount ; ++column)
        {
            Matrix::Mat *testMat = CreateTestMat(row, column, partCount, totalSize); 
            if(column == 0)
                CreateTestQuids(resultname + ".rowindices", db, row, partCount, totalSize);
            if(row == 0)
                CreateTestQuids(resultname + ".columnindices", db, column, partCount, totalSize);
            // write the matrix
            std::ostringstream filename;
            filename << resultname << ".precomputed";
            filename << ".part-R" << row << "-C" << column;
            filename << ".raw";
            WriteRaw(testMat, filename.str(), db, true);
            delete testMat;
        }
    return 1;
}

} // namespace Application
} // namespace Impala

int
main(int argc, char* argv[])
{
    Impala::Link::Mpi::Init(&argc, &argv);
    Impala::CmdOptions& options = Impala::CmdOptions::GetInstance();
    options.Initialise(false, false, true);
    options.AddOption(0, "inputFeatures", "specify the input features and their weights (INI file only, because prun messes up the command line)", "");

    int code = 1;
    if (options.ParseArgs(argc, argv, "<dataset> <result name>", 2))
    {
        if(options.GetArg(0) == "maketest")
            code = Impala::Application::makeTestMatrix(options);
        else
            code = Impala::Application::mainPrecomputeKernelMatrix(options);
    }

    Impala::Link::Mpi::Finalize();
    return code;
}

---