/*
    Copyright (2001) Ascent Software Ltd.
    Contact: Chris McGinlay, Roselynn, Levenwick, Shetland, UK, ZE2 9HX.

    This file is part of the Probability Distribution Namespace (version AUG2001).     
    
    Probability Distribution Namespace (version AUG2001) is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.     
    
    Probability Distribution Namespace (version AUG2001) is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.     
    
    You should have received a copy of the GNU General Public License
    along with Probability Distribution Namespace (version AUG2001); if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/************************************
  REVISION LOG ENTRY
  Revision By: Chris McGinlay
  Revised on 01/08/01 13:20:41
  Comments: adding Binomial class
 ************************************/

/************************************
  CREATION LOG ENTRY
  CREATED BY: Chris McGinlay
  CREATED ON: 17/07/01 21:35:39
  Comments: implementation for random classes
 ************************************/

#pragma warning(disable: 4786)  //For MS VC++6, STL generates C4786 warning during debug - can ignore
#include <cstdlib>  //for random declarations
#include <ctime>    //for seeding generator
#include <cmath>
#include <limits>   //for maximum size of long
#include <iostream> //for test()
#include <map>      //for test()

#include "prob_dist.h" //interface

namespace Prob_Dist = Probability_Distribution;

Prob_Dist::Uniform::Uniform(long lo, long hi, long seed)
{
    //seed generator
    if (!seed)
    {
        time_t thetime;  
        time(&thetime); //get a fairly random seed from system time
        seed = thetime;
    }
    srand(seed);

    if (!lo) lo = default_distrib._lo_bound;
    if (!hi) hi = default_distrib._up_bound;

    if ( lo >= hi ) throw Bad_Distribution();
    _lo_bound = lo;
    _up_bound = hi;
    _width = 1 + _up_bound - _lo_bound;
    return;
}

//return a random number on [lo, hi] with uniform probability
long Prob_Dist::Uniform::draw() const
{
    long result = 0;

    //get number and eliminate slight problem at upper bound
    do
    {
        result = long ( (double(rand()) / RAND_MAX) * _width ); 
    } while ( result== _width);

    return result + _lo_bound;    //on [_lo, _hi]
}

Prob_Dist::Uniform Prob_Dist::Uniform::default_distrib(0,RAND_MAX);  //simulated D6

void Prob_Dist::Uniform::set_default(int lo, int hi)
{
    if ( lo >= hi ) throw Bad_Distribution();
    default_distrib._lo_bound = lo;
    default_distrib._up_bound = hi;
    default_distrib._width = hi-lo+1;
    return;
}

void Prob_Dist::Uniform::test(int n) const
{
    std::map<int, int> bucket;
    for (int i=0; i<n; i++) bucket[draw()]++;
    std::cout << "Uniform Distribution" << '\t' << "On: [" << _lo_bound << ',' << _up_bound << ']' <<'\n';
    std::cout << "Bucket #" << '\t' << "Number of Results in Bucket" << '\n'; 
    for (int j=_lo_bound; j<=_up_bound; j++)
        std::cout << j << "\t\t" << bucket[j] << '\n';

    std::cout << std::endl; //flush
}

//*******************************************
//*******************************************
//*******************************************

const Prob_Dist::Uniform Prob_Dist::Binomial::_u(0,RAND_MAX,0);
long const Prob_Dist::Binomial::_long_max = std::numeric_limits<long>::max();

Prob_Dist::Binomial::Binomial(int n, double p, long seed)
{
    //seed generator
    if (!seed)
    {
        time_t thetime;  
        time(&thetime); //get a fairly random seed from system time
        seed = thetime;
    }
    srand(seed);

    if (!n) n = default_distrib._n;
    if (!p) p = default_distrib._p;

    if (n<0 || p<0.0 || p>1) throw Bad_Distribution();
    _n = n;
    _p = p;

    long f=0;   //temp for probability of i successful trials
    _bucket.resize(n+1);

    //now set up scaled bucket widths
    for (long i=0; i<_bucket.size(); i++)
    {
        //each delta_f is binomial prob density function * _long_max
        //following introduces a rounding (down) error, which has the useful side effect of preventing overflow in summation!

        f += bincoff(n, i) * pow(p,i) * pow(1-p, n-i) * _long_max;     //order of computation is important to prevent overflow
        _bucket.at(i) = f;
    }   

    return;
}

Prob_Dist::Binomial Prob_Dist::Binomial::default_distrib(6,0.5);

long Prob_Dist::Binomial::draw() const
{
    //draw from internal uniform distribution, dumping any result above upper bound of last _bucket
    long result=0;
    do
    {
        result = _u.draw() * (_long_max/RAND_MAX);   //rescale
    } while (result>_bucket.back()) ;

    int i=0;  //loop temp
    //loop until result is found in a bucket interval
    while ( result >= _bucket.at(i) ) i++;
    return long(i);
}

void Prob_Dist::Binomial::set_default(int n, double p)
{
    if (!n) n = default_distrib._n;
    if (!p) p = default_distrib._p;

    if (n<0 || p<0.0 || p>1.0) throw Bad_Distribution();

    default_distrib._n = n;
    default_distrib._p = p;

    long f=0;   //temp for probability of i successful trials
    default_distrib._bucket.resize(n+1);

    //now set up scaled bucket widths
    for (long i=0; i<default_distrib._bucket.size(); i++)
    {
        //each delta_f is binomial prob density function * _long_max
        //following introduces a rounding (down) error, which has the useful side effect of preventing overflow in summation!
        f += bincoff(n, i) * pow(p,i) * pow(1-p, n-i) * _long_max;     //order of computation is important to prevent overflow
        default_distrib._bucket.at(i) = f;
    }   

    return;
}

void Prob_Dist::Binomial::test(int n) const
{
    std::map<int, int> bucket;
    for (int i=0; i<n; i++) bucket[draw()]++;
    std::cout << "Binomial Distribution" << '\t' << "n: " << _n << '\t' << "p: " << _p << '\n';
    std::cout << "Bucket #" << '\t' << "Number of Results in Bucket" << '\n'; 
    for (int j=0; j<=_n; j++)
        std::cout << j << "\t\t" << bucket[j] << '\n';

    std::cout << std::endl; //flush
}

//*******************************************
//*******************************************
//*******************************************
//Poisson

Prob_Dist::Uniform const Prob_Dist::Poisson::_u(6,RAND_MAX,0);
long const Prob_Dist::Poisson::_long_max = std::numeric_limits<long>::max();

Prob_Dist::Poisson::Poisson(double mu, long seed)
{
    //seed generator
    if (!seed)
    {
        time_t thetime;  
        time(&thetime); //get a fairly random seed from system time
        seed = thetime;
    }
    srand(seed);

    if (!mu) mu = default_distrib._mu;

    if (mu<0) throw Bad_Distribution();

    _mu = mu;

    long i=0;   //i successful outcomes (number of)
    long f=0;   //temp for probability of i successful trials
    long width=0;   //another temp

    // _bucket gets bigger until f becomes "small" (dist has infinite tail)
    // note differences in this part in other classes

    //now set up scaled bucket widths

    long max = 0;
    do
    {
        width = double(pow(_mu,i) / factorial(i)) * exp(-_mu) * _long_max;
        if (width <0) throw Bad_Distribution();
        max = width>max ? width : max;
        
        f+=width;        
        _bucket.push_back(f);
        i++;
    } while (width > 0);

    return;
}

Prob_Dist::Poisson Prob_Dist::Poisson::default_distrib(7);

long Prob_Dist::Poisson::draw() const
{
    //draw from internal uniform distribution, dumping any result above upper bound of last _bucket
    long result=0;
    do
    {
        result = _u.draw() * int(_long_max/RAND_MAX);   //rescale
    } while (result>_bucket.back()) ;

    int i=0;  //loop temp
    //loop until result is found in a bucket interval
    while ( result >= _bucket.at(i) ) i++;
    return i;
}

void Prob_Dist::Poisson::set_default(double mu)
{
    if (mu<0) throw Bad_Distribution();

    default_distrib._mu = mu;

    long i=0;   //i successful outcomes (number of)
    long f=0;   //temp for probability of i successful trials
    long width=0;   //another temp

    // _bucket gets bigger until f becomes "small" (dist has infinite tail)
    // note differences in this part in other classes

    //now set up scaled bucket widths

    long max = 0;
    do
    {
        width = double(pow(mu,i) / factorial(i)) * exp(-mu) * _long_max;
        if (width <0) throw Bad_Distribution();
        max = width>max ? width : max;
        
        f+=width;        
        default_distrib._bucket.push_back(f);
        i++;
    } while (width > 0);

    return;
}

void Prob_Dist::Poisson::test(int n) const
{
    std::map<int, int> bucket;
    for (int i=0; i<n; i++) bucket[draw()]++;
    std::cout << "Poisson Distribution" << '\t' << "Mu: " << _mu << '\n';
    std::cout << "Bucket #" << '\t' << "Number of Results in Bucket" << '\n'; 
    int size = bucket.size();
    for (int j=0; j<=size; j++)
        std::cout << j << "\t\t" << bucket[j] << '\n';

    std::cout << std::endl; //flush
}
//*******************************************
//*******************************************
//*******************************************

long double Prob_Dist::factorial(int n)
{
    if (n) return n*factorial(n-1);
    return 1;
}

long double Prob_Dist::factorial(int top, int bot)
{
    if (top!=bot) return top*factorial(top-1, bot);
    return 1;
}

long Prob_Dist::bincoff(int n, int r)
{
    //choose fastest route to computation
    if (n==r || r==0) return 1; //yikes!
    if (r>n-r)
    {
        //r is greater than n-r
        return factorial(n,r)/factorial(n-r);
    }

    //n-r is greater than or equal to r
    return factorial(n, n-r)/factorial(r);
}