D:/~Research/Code/C++/EvolveTraffic/Distribution.cpp

Go to the documentation of this file.

// Distribution.cpp : implementation file
//

#include "stdafx.h"
#include "EvolveTraffic.h"
#include "Distribution.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

// CDistribution

IMPLEMENT_SERIAL(CDistribution, CObject, VERSION_NUMBER)

CDistribution::CDistribution():PI(3.14159265359)
{
        m_DistributionID = DIST_NORMAL;
        m_Location = 10.0;
        m_Scale = 2.0;
        m_Shape = 0.0;
}

CDistribution::CDistribution(double loc, double sc, double sh):PI(3.14159265359)
{
        m_DistributionID = DIST_NORMAL;
        m_Location = loc;
        m_Scale = sc;
        m_Shape = sh;
}

CDistribution::~CDistribution()
{

}

CDistribution::CDistribution(const CDistribution& dist):PI(3.14159265359)
{
        this->setDistributionID( dist.getDistributionID() );
        this->setLocation( dist.getLocation() );
        this->setScale( dist.getScale() );
        this->setShape( dist.getShape() );
}

// Initialise static member of class
MTRand CDistribution::m_RNG;    // CDistribution::m_RNG(500); MAGIC NUMBER - constant seed if needed for debugging

// CDistribution diagnostics

#ifdef _DEBUG
void CDistribution::AssertValid() const
{
        CObject::AssertValid();
}

void CDistribution::Dump(CDumpContext& dc) const
{
        CObject::Dump(dc);
}
#endif //_DEBUG

// CDistribution serialization

void CDistribution::Serialize(CArchive& ar)
{
        if (ar.IsStoring())
        {
                ar      << m_DistributionID
                        << m_Location
                        << m_Scale
                        << m_Shape;
        }
        else
        {
                ar      >> m_DistributionID
                        >> m_Location
                        >> m_Scale
                        >> m_Shape;
        }
}

// CDistribution commands

CDistribution& CDistribution::operator=(const CDistribution &dist)
{
        this->setDistributionID( dist.getDistributionID() );
        this->setLocation( dist.getLocation() );
        this->setScale( dist.getScale() );
        this->setShape( dist.getShape() );

        return *this;
}

WORD CDistribution::getDistributionID() const
{
        return m_DistributionID;
}

double CDistribution::getLocation() const
{
        return m_Location;
}       

double CDistribution::getScale() const
{
        return m_Scale;
}

double CDistribution::getShape() const
{
        return m_Shape;
}

void CDistribution::setDistributionID(WORD id)
{
        m_DistributionID = id;
}

void CDistribution::setLocation(double loc)
{
        m_Location = loc;
}

void CDistribution::setScale(double sc)
{
        m_Scale = sc;
}

void CDistribution::setShape(double sh)
{
        m_Shape = sh;
}

// This function generates a value from the distribution
double CDistribution::Generate()
{
        switch(m_DistributionID)
        {
                case DIST_EXPONENTIAL:
                        return GenerateExponential();
                case DIST_LOGNORMAL:
                        return GenerateLogNormal();
                case DIST_GAMMA:
                        return GenerateGamma();
                case DIST_GUMBEL:
                        return GenerateGumbel();
                case DIST_POISSON:
                        return GeneratePoisson();
                case DIST_GEV:
                        return GenerateGEV();
                case DIST_NORMAL:
                        return GenerateNormal();
                case DIST_CONST:
                        return m_Location;
                default:
                        return m_Location;              // constant
        }
}

// The random deviate generators

double CDistribution::GenerateExponential()
{
        double u01 = m_RNG.rand();
        double val = -log(u01);
        return val;
}

double CDistribution::GenerateLogNormal()
{
        double x[2], u[2];
        u[0] = m_RNG.rand();
        u[1] = m_RNG.rand();
        x[0] = m_Location + m_Scale * sqrt(-2 * log(u[0])) * cos(2 * PI * u[1]);
        x[1] = m_Location + m_Scale * sqrt(-2 * log(u[0])) * sin(2 * PI * u[1]);
        double val = exp( x[ (int)(u[1] + 0.5) ] );

        return val;
}

double CDistribution::GenerateGamma()
{
        double u[2];
        u[0] = m_RNG.rand();
        u[1] = m_RNG.rand();
        u[2] = m_RNG.rand();
        double val = (1/m_Location) 
                                        * (-log(u[2])) 
                                        * pow( u[0], (1/m_Scale) ) 
                                        / pow( u[0], (1/m_Scale) 
                                        + pow( u[1], (1/(1 - m_Scale)) ) );
        
        return val;
}

double CDistribution::GenerateGumbel()
{
        double u01 = m_RNG.rand();
        double val = m_Location - (1 / m_Scale) * log(log(1 / u01));

        return val;
}

double CDistribution::GeneratePoisson()
{
        double x[2], u[2];
        u[0] = m_RNG.rand();
        u[1] = m_RNG.rand();

        double m = m_Location - 0.5;
        double sigma = sqrt(m_Scale);

        x[0] = m + sigma * sqrt(-2 * log(u[0])) * cos(2 * PI * u[1]);
        x[1] = m + sigma * sqrt(-2 * log(u[0])) * sin(2 * PI * u[1]);
        double val = x[ (int)(u[1] + 0.5) ];

        return val;
}

double CDistribution::GenerateGEV()
{
        double u01 = m_RNG.rand();
        double val = m_Location + m_Scale * (1 - pow(-log(u01),m_Shape)) / m_Shape;

        return val;
}

double CDistribution::GenerateNormal()
{
        return m_Location + m_Scale * BoxMuller();
}

double CDistribution::BoxMuller()
{
        // normal random variate generator
        
        double x1, x2, w, y1;
        static double y2;
        static int use_last = 0;

        if (use_last)                   // use value from previous call
        {
                y1 = y2;
                use_last = 0;
        }
        else
        {
                do {
                        x1 = 2.0 * m_RNG.rand() - 1.0;
                        x2 = 2.0 * m_RNG.rand() - 1.0;
                        w = x1 * x1 + x2 * x2;
                } while ( w >= 1.0 );

                w = sqrt( (-2.0 * log( w ) ) / w );
                y1 = x1 * w;
                y2 = x2 * w;
                use_last = 1;
        }
        return( y1 );
}

double CDistribution::GenKeepingCoV(double newLocation)
{
        const double zero = 0.0001;
        double val = m_Location;
        if( (m_Scale >  zero || m_Scale < -zero) && (m_Location >  zero || m_Location < -zero) )
        {
                double CoV = m_Scale / m_Location;
                double OldScale = m_Scale;
                double OldLocation = m_Location;
                
                // set new params & generate new value
                m_Scale = newLocation * CoV;
                m_Location = newLocation;
                val = Generate();
                // restore old params
                m_Scale = OldScale;
                m_Location = OldLocation;
        }
        return val;
}

---