Prob.cpp

#include "Headers.h"
#include "Prob.h"
#include "Gamma.h"

#define INTEGRAL_PRECISION  1

CProb::CProb(
    double      alpha,
    double      xmin,
    double      xmax,
    const char* cacheFile
    )
{
    this->alpha = alpha;
    this->xmin = xmin;
    this->xmax = xmax;

    this->hyp = new Hyp2F1(1, 2-2*this->alpha, 2-this->alpha);

    this->c1 = Gamma(1-this->alpha) / Gamma(2-this->alpha);
    this->c2 = ((1-this->alpha)*(1-this->alpha)/pow((pow(this->xmax,1-this->alpha)-pow(this->xmin,1-this->alpha)), 2));
    this->c3 = (1-this->alpha)*(1-this->alpha)/(pow(this->xmax,1-this->alpha)-pow(this->xmin,1-this->alpha));
    this->c4 = pow(this->xmax,1-this->alpha);

    // Populate cache from the cache file
    ifstream file;
    try
    {
        file.open(cacheFile, ios_base::in | ios_base::binary);

        cout << "\n[i] Reading abandon probability for the cache file \'" << cacheFile << "\'..." << endl;

        file.read((char*)&this->cacheEpiMax, sizeof(this->cacheEpiMax));
        file.read((char*)&this->cacheEciMax, sizeof(this->cacheEciMax));
        file.read((char*)&this->cacheAlpha, sizeof(this->cacheAlpha));
        file.read((char*)&this->cacheXmin, sizeof(this->cacheXmin));
        file.read((char*)&this->cacheXmax, sizeof(this->cacheXmax));
        file.read((char*)&this->cacheEpMin, sizeof(this->cacheEpMin));
        file.read((char*)&this->cacheEpMax, sizeof(this->cacheEpMax));
        file.read((char*)&this->cacheEpDelta, sizeof(this->cacheEpDelta));
        file.read((char*)&this->cacheEcMin, sizeof(this->cacheEcMin));
        file.read((char*)&this->cacheEcMax, sizeof(this->cacheEcMax));
        file.read((char*)&this->cacheEcDelta, sizeof(this->cacheEcDelta));

        cout << "[i] Cache file parameters are: alpha=" << this->cacheAlpha << " xmin=" << this->cacheXmin << " xmax=" << this->cacheXmax <<
            " ep_min=" << this->cacheEpMin << " ep_max=" << this->cacheEpMax << " ep_delta=" << this->cacheEpDelta <<
            " ec_min=" << this->cacheEcMin << " ec_max=" << this->cacheEcMax << " ec_delta=" << this->cacheEcDelta << endl;

        // Check the cache file
        if((this->cacheAlpha == this->alpha) && (this->cacheXmin == this->xmin) && (this->cacheXmax == this->xmax))
        {
            // Allocate the cache space
            this->cache = new double[this->cacheEciMax*this->cacheEpiMax];

            cout << "[i] Cache file okay. Reading cache data...";
            cout.flush();

            file.read((char*)this->cache, sizeof(double)*this->cacheEpiMax*this->cacheEciMax);

            cout << "done." << endl;
        }
        else
        {
            cout << "[!] Cache file mismatch (alpha=" << this->alpha << " xmin=" << xmin << " xmax=" << xmax << ")." << endl;
            this->CacheDefault();
        }

        file.close();
    }
    catch(...)
    {
        cout << "\n[!] Reading abandon probability from the cache file \'" << cacheFile << "\' failed." << endl;
        this->CacheDefault();
    }

    assert(this->cache);
}

CProb::~CProb()
{
    delete this->hyp;
    delete[] this->cache;
}

void CProb::CacheDefault()
{
    this->cacheAlpha = this->alpha;
    this->cacheXmin = cacheXmin;
    this->cacheXmax = cacheXmax;
    this->cacheEpMin = 0.0;
    this->cacheEpMax = this->xmax;
    this->cacheEpDelta = 1.0;
    this->cacheEcMin = 0.0;
    this->cacheEcMax = this->xmax;
    this->cacheEcDelta = 1.0;

    // Calculate the cache size
    this->cacheEpiMax = (unsigned int)((this->cacheEpMax - this->cacheEpMin) / this->cacheEpDelta) + 1;
    this->cacheEciMax = (unsigned int)((this->cacheEcMax - this->cacheEcMin) / this->cacheEcDelta) + 1;

    // Allocate the cache
    this->cache = new double[this->cacheEciMax*this->cacheEpiMax];

    // Set all cache values to default values
    for(unsigned int index = 0; index < this->cacheEciMax*this->cacheEpiMax; index++)
        this->cache[index] = -1.0;

    cout << "[i] Creating default probability cache with parameters: alpha=" << this->cacheAlpha << " xmin=" << this->cacheXmin << " xmax=" << this->cacheXmax <<
            " ep_min=" << this->cacheEpMin << " ep_max=" << this->cacheEpMax << " ep_delta=" << this->cacheEpDelta <<
            " ec_min=" << this->cacheEcMin << " ec_max=" << this->cacheEcMax << " ec_delta=" << this->cacheEcDelta << endl;
}

double CProb::Prob(double ep, double ec)
{
    // Calculate the cache index (do not go outside of cache boundary)
    if(ep > this->cacheEpMax) ep = this->cacheEpMax;
    if(ec > this->cacheEcMax) ec = this->cacheEcMax;

    if(ep < this->cacheEpMin) ep = this->cacheEpMin;
    if(ec < this->cacheEcMin) ec = this->cacheEcMin;

    unsigned int epi = (unsigned int)((ep - this->cacheEpMin)/this->cacheEpDelta);
    unsigned int eci = (unsigned int)((ec - this->cacheEcMin)/this->cacheEcDelta);

    assert(epi < this->cacheEpiMax);
    assert(eci < this->cacheEciMax);

    unsigned int cindex = eci * this->cacheEpiMax + epi;

    /*
     * Part I. If value is in the probability in the cache, return the cached value
     */
    if(this->cache[cindex] >= 0) return this->cache[cindex];

    /*
     * Part II. Else, calculate new value and insert it in the probability cache
     */
    this->mutex.Lock();
    this->cache[cindex] = this->InvCdf0(epi, eci);
    this->mutex.Unlock();

    return this->cache[cindex];
}

double CProb::InvCdf0(double ep, double ec)
{
    double f = 0;
    double z;

    if(ec < this->xmin)
        if(ep < this->xmin)
        {
            // ep < this->xmin, ec < this->xmin
            z = (0 > this->xmin-this->xmax+ec-ep)?0:this->xmin-this->xmax+ec-ep;
            for(; z < ec-ep; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ec,z,ep,ec) - this->Integral(this->xmin-ep-z,z,ep,ec)) * this->c2;
            for(; z <= this->xmax-this->xmin+ec-ep; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ep-z,z,ep,ec) - this->Integral(this->xmin-ec,z,ep,ec)) * this->c2;
        }
        else
        {
            double c = this->c3 / (this->c4-pow(ep,1-this->alpha));

            // ep >= this->xmin, ec < this->xmin
            z = (0 > ec-this->xmax)?0:ec-this->xmax;
            for(; z < ec-ep; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ec,z,ep,ec) - this->Integral(-z,z,ep,ec)) * c;
            for(; z < ec-this->xmin; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ep-z,z,ep,ec) - this->Integral(-z,z,ep,ec)) * c;
            for(; z < this->xmax-this->xmin+ec-ep; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ep-z,z,ep,ec) - this->Integral(this->xmin-ec,z,ep,ec)) * c;
        }
    else
        if(ep < this->xmin)
        {
            double c = this->c3 / (this->c4-pow(ec,1-this->alpha));

            // ep < this->xmin, ec >= this->xmin
            z = (0 > this->xmin-this->xmax+ec-ep)?0:this->xmin-this->xmax+ec-ep;
            for(; z < this->xmin-ep; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ec,z,ep,ec) - this->Integral(this->xmin-ep-z,z,ep,ec)) * c;
            for(; z < ec-ep; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ec,z,ep,ec) - this->Integral(0,z,ep,ec)) * c;
            for(; z < this->xmax-ep; z += INTEGRAL_PRECISION)
                f += fabs(this->Integral(this->xmax-ep-z,z,ep,ec) - this->Integral(0,z,ep,ec)) * c;
        }
        else
            if(ep < ec)
            {
                double c = (1-this->alpha)*(1-this->alpha)/((this->c4-pow(ep,1-this->alpha))*(this->c4-pow(ec,1-this->alpha)));

                // ep >= this->xmin, ec >= this->xmin, ep < ec
                z = 0;
                for(; z < ec-ep; z += INTEGRAL_PRECISION)
                    f += fabs(this->Integral(this->xmax-ec,z,ep,ec) - this->Integral(0,z,ep,ec)) * c;
                for(; z < this->xmax-ep; z += INTEGRAL_PRECISION)
                    f += fabs(this->Integral(this->xmax-ep-z,z,ep,ec) - this->Integral(0,z,ep,ec)) * c;
            }
            else
            {
                double c = (1-this->alpha)*(1-this->alpha)/((this->c4-pow(ep,1-this->alpha))*(this->c4-pow(ec,1-this->alpha)));

                // ep >= this->xmin, ec >= this->xmin, ep >= ec
                z = (0 > ec-this->xmax)?0:ec-this->xmax;
                for(; z < ec-ep; z += INTEGRAL_PRECISION)
                    f += fabs(this->Integral(this->xmax-ec,z,ep,ec) - this->Integral(-z,z,ep,ec)) * c;
                for(; z < 0; z += INTEGRAL_PRECISION)
                    f += fabs(this->Integral(this->xmax-ep-z,z,ep,ec) - this->Integral(-z,z,ep,ec)) * c;
                for(; z < this->xmax-ep; z += INTEGRAL_PRECISION)
                    f += fabs(this->Integral(this->xmax-ep-z,z,ep,ec) - this->Integral(0,z,ep,ec)) * c;
            }
    return f;
}

double CProb::Integral(double x, double z, double ep, double ec)
{
    double b = z + ep;
    double c = ec;

    if(b != c) return this->c1 * pow((x+b)*(x+c), 1-this->alpha) * fabs((*hyp)(-(x+c)/(b-c))) / (b-c);
    else return pow(x+b,1-2*this->alpha)/(1-2*this->alpha);
}