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

---

mainPrecomputeKernelMatrix.cpp

Go to the documentation of this file.

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

#include "Persistency/DistributedAccessRepository.h"
#include "Persistency/FeatureTableRepository.h"
#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/VirtualMatrixFactory.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"
#ifdef SSE_USED
#include "Core/Vector/Chi2DistanceSSE.h"
#endif
#include "Core/Vector/Apply.h"
#include "Core/VideoSet/MakeVideoSet.h"
#include "Core/ImageSet/MakeImageSet.h"
#include "Core/Vector/Chi2Distance.h"
#include "Core/Vector/HistogramIntersection.h"
#include "Core/Vector/DotProduct.h"
#include "Util/PropertySet.h"
#include "Core/Database/MakeRawDataSet.h"
#include "Util/TimeStats.h"

namespace Impala
{
namespace Application
{
namespace Precompute
{


using namespace Core;
using namespace Persistency;
using Core::Database::RawDataSet;
using Core::Matrix::DistributedAccess;

bool
CheckFeatures(CmdOptions& options, RawDataSet* thisSet, DistributedAccess* da)
{
    ILOG_VAR(Impala.Application.Precompute.CheckFeatures);

    /* 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_HEADNODE("\"" << inputFeatures <<
                                    "\" contains a zero weight");
                return false;
            }
            String name = p.GetString(' ', false);
            if (name.empty())
            {
                ILOG_ERROR_HEADNODE("inputFeatures: no name for weight");
                return false;
            }
            ILOG_DEBUG("n: "<< name);
            String indexCat = options.GetString("featureIndexCat");
            FeatureLocator loc(thisSet->GetLocator(), false, true, indexCat,
                               name, "");
            if (!FeatureTableRepository().Exists(loc))
            {
                ILOG_ERROR_HEADNODE("Could not open \"" << name << "\"");
                return false;
            }
            da->AddFeature(name, weight, -1);
        }
    }
    return true;
}

Feature::FeatureTable*
OpenFeatureTable(Feature::FeatureDefinition& featureDef, String indexCat,
                 RawDataSet* dataset)
{
    ILOG_VAR(Impala.Application.Precompute.OpenFeatureTable);
    FeatureLocator loc(dataset->GetLocator(), false, true, indexCat,
                       featureDef.AsString(), "");
    loc.SetUseBroadcast(true);
    Feature::FeatureTable* f = FeatureTableRepository().Get(loc);
    ILOG_INFO_HEADNODE(featureDef.AsString() << ", size = " << f->Size() <<
                       "; length = " << f->GetFeatureVectorLength());
    return f;
}

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


Matrix::Mat*
ComputeMatrix(DistributedAccess* da, Feature::FeatureTable* devel,
              Feature::FeatureTable* test, String distanceFunction)
{
    ILOG_VAR(Impala.Application.Precompute.ComputeMatrix);
    Feature::FeatureTable* horizontal;
    Feature::FeatureTable* vertical;
    if (da->GetTotalPartCount() == 1)
    {
        horizontal = devel;
        vertical = test;
    }
    else
    {
        int myPart = Link::Mpi::MyId();
        int column = da->GetColumnPartOfPart(myPart);
        int row = da->GetRowPartOfPart(myPart);
        horizontal = GetPartial(devel, column, da->GetColumnPartCount());
        vertical = GetPartial(test, row, da->GetRowPartCount());
        delete devel;
        if (test != devel)
            delete test;
    }

    // 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 = 0;
    if (distanceFunction == "dot")
    {
        distanceMatrix = Apply(&Vector::DotProduct, vectorsH, vectorsV);
    }
    else if (distanceFunction == "histogramintersection")
    {
        distanceMatrix = Apply(&Vector::HistogramIntersectionDouble,
                               vectorsH, vectorsV);
    }
    else
    {
#ifdef SSE_USED
        distanceMatrix = Apply(&Vector::Chi2DistanceSSE, vectorsH, vectorsV);
#else
        distanceMatrix = Apply(&Vector::Chi2Distance, vectorsH, vectorsV);
#endif
    }
    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;
}

int
Precompute(CmdOptions& options)
{
    ILOG_VAR(Impala.Application.Precompute.Precompute);
    RawDataSet* thisSet = Core::Database::MakeRawDataSet(options.GetArg(0));
    if (!thisSet)
    {
        ILOG_INFO_HEADNODE("Failed to open thisSet " << options.GetArg(0));
        return 1;
    }
    
    String develSetName = options.GetArg(1);
    ILOG_INFO_HEADNODE("DevelSetName = [" << develSetName << "]");
    RawDataSet* develSet = 0;
    if (!develSetName.empty())
    {
        develSet = Core::Database::MakeRawDataSet(develSetName);
        if (!develSet)
        {
            ILOG_INFO_HEADNODE("Failed to open develSet " << develSetName);
            return 1;
        }
    }
    String model = options.GetArg(2);
    String kernel = options.GetArg(3);
    String indexCat = options.GetString("featureIndexCat");
    KernelMatrixLocator daLoc(thisSet->GetLocator(), true, "",
                              develSetName, model, kernel, "");
    daLoc.SetFeatureIndexCat(indexCat);
    if (!options.GetBool("override") &&
        DistributedAccessRepository().Exists(daLoc))
    {
        ILOG_INFO("Kernel already exists.");
        return 0;
    }

    int startFeature = options.GetInt("startFeature");
    DistributedAccess* da = 0;
    if (startFeature == 0)
    {
        da = new DistributedAccess();
        if (!da->Valid(true, false))
            return 1;

        ILOG_INFO_HEADNODE("Using a " << da->GetRowPartCount() << "x" <<
                           da->GetColumnPartCount() << "grid");
        if (!CheckFeatures(options, thisSet, da))
            return 1;
    }
    else
    {
        daLoc.SetIncrement(startFeature);
        daLoc.SetStartNode(0);
        daLoc.SetNodeCount(Link::Mpi::NrProcs());
        da = DistributedAccessRepository().Get(daLoc);
    }


    if (develSet)
    {
        KernelMatrixLocator devLoc(develSet->GetLocator(), true, "", "",
                                   model, kernel, "");
        devLoc.SetFeatureIndexCat(indexCat);
        devLoc.SetDoParts(0);
        devLoc.SetStartNode(0);
        devLoc.SetNodeCount(Link::Mpi::NrProcs());
        DistributedAccess* devDA = DistributedAccessRepository().Get(devLoc);
        for (int i = 0; i < devDA->GetNrFeatures(); i++)
        {
            double average = devDA->GetAverage(i);
            ILOG_INFO_HEADNODE("average " << i << " " << average);
            da->SetAverage(i, average);
        }
        // Could do some more checking here...
        delete devDA;
    }
    else
    {
        if (model == "chi2")
            da->SetHasOwnAverages(true);
    }

    String distanceFunction = options.GetString("distanceFunction");
    int numberFeatures = options.GetInt("numberFeatures");
    if (numberFeatures == -1)
        numberFeatures = da->GetNrFeatures() - startFeature;
    if (startFeature + numberFeatures >= da->GetNrFeatures())
        numberFeatures = da->GetNrFeatures() - startFeature;
    Matrix::Mat* accumulator = 0;
    if (startFeature != 0)
        accumulator = da->StealPart();
    double totalweight = 0;
    for (int i=0 ; i<startFeature ; i++)
        totalweight += da->GetWeight(i);
    Util::TimeStats statsCompute;
    statsCompute.AddGroup("read features");
    statsCompute.AddGroup("compute");
    Util::TimeStats statsOverall;
    statsOverall.AddGroupsFromSub(&statsCompute);
    statsOverall.AddGroup("final exp");
    statsOverall.AddGroup("write");
    statsOverall.MeasureFirst();
    for (int i=startFeature ; i<startFeature+numberFeatures ; i++)
    {
        statsCompute.MeasureFirst();
        // open feature tables
        double weight = da->GetWeight(i);
        Feature::FeatureDefinition fDef(da->GetFeature(i));
        Feature::FeatureTable* f2 = OpenFeatureTable(fDef, indexCat, thisSet);
        Feature::FeatureTable* f1 = f2;
        if (develSet)
            f1 = OpenFeatureTable(fDef, indexCat, develSet);

        if ((Link::Mpi::MyId() == 0) && (i == 0))
        {
            Table::Copy(da->GetColumnQuids(), f1);
            da->SetColumns(f1->Size());
            Table::Copy(da->GetRowQuids(), f2);
            da->SetRows(f2->Size());
        }

        statsCompute.MeasureNext();
        // compute kernel distances between features
        Matrix::Mat *distanceMatrix = ComputeMatrix(da, f1, f2,
                                                    distanceFunction);
        // 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
        // Note : hik kernel does not use averages
        if (da->GetHasOwnAverages())
        {
            if (distanceFunction == "dot")
            {
                da->SetAverage(i, 0.0);
            }
            else
            {
                double average = GetAverage(distanceMatrix);
                da->SetAverage(i, average);
            }
        }

        // accumulate; this is the part inside the 'exp' in the kernel function
        MulVal(distanceMatrix, distanceMatrix, weight);
        if (distanceFunction == "chi2")
        {
            DivVal(distanceMatrix, distanceMatrix, -da->GetAverage(i));
        }
        if (accumulator == 0)
            Set(accumulator, distanceMatrix);
        else
            Add(accumulator, accumulator, distanceMatrix);
        delete distanceMatrix;

        totalweight += weight;
        statsCompute.MeasureLast();
        ILOG_INFO_HEADNODE("Compute " << statsCompute.AsString());
        ILOG_INFO_HEADNODE("Memory " << Process::MemoryInfo::GetUsageString());
    }
    statsOverall.MeasureFromSub(&statsCompute);

    ILOG_INFO_HEADNODE("finalising...");
    if (startFeature + numberFeatures == da->GetNrFeatures())
    {
        // only thing to do is weigh and exp
        DivVal(accumulator, accumulator, totalweight);
        if (distanceFunction == "chi2")
        {
            Exp(accumulator, accumulator);
        }
        daLoc.SetIncrement(0);
        daLoc.SetWriteReal32(true);
    }
    else
    {
        daLoc.SetIncrement(startFeature + numberFeatures);
    }

    statsOverall.MeasureNext();
    typedef Matrix::VirtualMatrixFactory VirtualMatrixFactory;
    VirtualMatrixFactory& vmf = VirtualMatrixFactory::GetInstance();
    Matrix::VirtualMatrix* vm = vmf.ConstructMemory(accumulator);
    da->AddPart(vm);
    DistributedAccessRepository().Add(daLoc, da);
    statsOverall.MeasureLast();
    ILOG_INFO_HEADNODE("Overall " << statsOverall.AsString());
    delete da;
    return 0;
}

int
ChopMatrix(CmdOptions& options)
{
    ILOG_VAR(Impala.Application.Precompute.ChopMatrix);
    RawDataSet* thisSet = Core::Database::MakeRawDataSet(options.GetArg(0));
    if (!thisSet)
    {
        ILOG_INFO_HEADNODE("failed to open dataset " << options.GetArg(0));
        return 1;
    }
    RawDataSet* develSet = 0;
    String develSetName = options.GetArg(1);
    String model = options.GetArg(2);
    String kernel = options.GetArg(3);
    int nrParts = atol(options.GetArg(5));
    KernelMatrixLocator srcLoc(thisSet->GetLocator(), true, "",
                               develSetName, model, kernel, "");
    srcLoc.SetFeatureIndexCat(options.GetString("featureIndexCat"));

    srcLoc.SetStartNode(0);
    srcLoc.SetNodeCount(Link::Mpi::NrProcs());
    srcLoc.SetDoParts(2);
    DistributedAccess* srcDA = DistributedAccessRepository().Get(srcLoc);
    srcDA->Dump();

    Real64* lineBuf = new Real64[srcDA->NrCol()];
    typedef Matrix::Mat32 Mat32;
    DistributedAccess dstDA(srcDA->NrRow(), srcDA->NrCol(), nrParts, nrParts,
                            0, 1);
    dstDA.CopyQuidsFrom(srcDA);
    dstDA.CopyFeaturesFrom(srcDA);
    KernelMatrixLocator dstLoc = srcLoc;
    dstLoc.SetWalkType("chopped");
    dstLoc.SetDoParts(0);
    DistributedAccessRepository().Add(dstLoc, &dstDA);

    for (int row=0 ; row<nrParts ; row++)
    {
        for (int col=0 ; col<nrParts ; col++)
        {
            int startY = dstDA.GetRowStartOfPart(row);
            int endY = dstDA.GetRowEndOfPart(row);
            int startX = dstDA.GetColumnStartOfPart(col);
            int endX = dstDA.GetColumnEndOfPart(col);
            ILOG_INFO("row=" << row << ", col=" << col << ", startX=" << startX
                      << ", endX=" << endX << ", startY=:" << startY
                      << ", endY=" << endY);
            Mat32* mat = Matrix::MatCreate<Mat32>(endY - startY, endX - startX);
            for (int y=startY ; y<endY ; y++)
            {
                srcDA->GetRow(y, lineBuf, srcDA->NrCol());
                float* ptr = mat->CPB(0, y-startY);
                for (int x=startX ; x<endX ; x++)
                    ptr[x-startX] = lineBuf[x];
            }
            Array::PrintDataCorners(mat, 2, 2);

            String filename = "PrecomputedKernels/chopped/" + model + "/"
                + kernel + "/" + kernel + ".precomputed.part-R" + MakeString(row)
                + "-C" + MakeString(col) + ".raw" ;
            Persistency::File file = RepositoryInFileSystem::GetInstance().
                GetFile(dstLoc, "", filename, true, false);
            ILOG_INFO("Saving part in " << file.GetPath());
            Array::WriteRaw(mat, file, 1);
            delete mat;
        }
    }

    delete lineBuf;
    delete srcDA;
    return 0;
}

int
MainPrecompute(int argc, char* argv[])
{
    Link::Mpi::Init(&argc, &argv);
    CmdOptions& options = Impala::CmdOptions::GetInstance();
    options.Initialise(false, false, true);
    options.AddOption(0, "startFeature", "idx", "0");
    options.AddOption(0, "numberFeatures", "nr", "-1");
    options.AddOption(0, "featureIndexCat", "name", "");
    options.AddOption(0, "distanceFunction", "dot|chi2|histogramintersection", "chi2");

    if (! options.ParseArgs(argc, argv, "<dataset> <develset> <model> <kernel>", 4))
    {
        Link::Mpi::Finalize();
        return 1;
    }
        
    ILOG_VAR(Impala.Application.Precompute.Main);

    Timer timer;
    if ((options.GetNrArg() == 6) && (options.GetArg(4) == "chop"))
        ChopMatrix(options);
    else
        Precompute(options);

    ILOG_INFO_HEADNODE("Total execution time = " << timer.SplitTimeStr());
    int nrOfErrors = ILOG_ERROR_COUNT;
    nrOfErrors = Link::Mpi::ReduceSum(nrOfErrors);
    ILOG_INFO_HEADNODE("Root: total nr error = " << nrOfErrors);

    Link::Mpi::Finalize();
    return nrOfErrors;
}

} // namespace Precompute
} // namespace Application
} // namespace Impala

int
main(int argc, char* argv[])
{
    return Impala::Application::Precompute::MainPrecompute(argc, argv);
}

---