#include "misc.h"

/*
int execute_event(event * myEvent)
{
//downcast event object to the proper type
    for(i = 0; i < clock.maxtimestep; i += delta)

}*/


/*
void error (lua_State *L, const char *fmt, ...) {
    va_list argp;
    va_start(argp, fmt);
    vfprintf(stderr, fmt, argp);
    va_end(argp);
    lua_close(L);
    exit(EXIT_FAILURE);
}
*/

// Allows to select the number of threads used
// Static fields allow to initialise on first launch 
#define MULTI_THREAD 4
#if MULTI_THREAD == 0
    double Random(unsigned int seed)
    {
        static mt19937 gen(seed);
        static uniform_real_distribution<double> dist(0.0, 1.0);
        return dist(gen);
    }
#else
    double Random(unsigned int seed)
    {
        static unsigned int nt = MULTI_THREAD;
        static mt19937 gen(seed);
        static uniform_real_distribution<double> dist(0.0, 1.0);
        static RandEngine<mt19937, uniform_real_distribution<double>, 
                          double> dg(gen, dist);
        static ThreadedGen<RandEngine<mt19937, uniform_real_distribution<double>, 
                           double>, double> r1(10000*nt, 10*nt, nt, dg);
        return r1.get_value();
    }
#endif

/*
// TODO: change to cpp 11 random http://en.cppreference.com/w/cpp/numeric/random
// uniforme entre 0 compris et 1 non compris 
double Random() 
{
// Générateur de nombre pseudo-aléatoires. Utilise les seeds S1 et S2 et les
// modifie à chaque appel.
// Reference : L'ECUYER, 1988 in Goudet, J., 1993, The genetics of geographically
// structured populations, School of Biological Sciences University College of 
//North Wales.
    long i;
    double r, z;

    r = S1 / 53668;
    i = (long) r;
    r = (double) i;
    S1 = 40014 * (S1 - r * 53668) - r * 12211;
    if (S1 < 0)
        S1 += 2147483563;
    r = S2 / 52774;
    i = (long) r;
    r = (double) i;
    S2 = 40692 * (S2 - r * 52774) - r * 3791;
    if (S2 < 0)
        S2 += 2147483399;
    z = S1 - S2;
    if (z < 1)
        z += 2147483562;

    return (z * 4.656613e-10);
}
*/
double gasdev()
{
    static int iset = 0;
    static double gset;
    double fac, rsq, v1, v2;

    if (iset == 0)
    {
        do
        {
            v1 = 2.0 * Random() - 1.0;
            v2 = 2.0 * Random() - 1.0;
            rsq = v1 * v1 + v2*v2;
        } while ((rsq >= 1.0) || (rsq == 0.0));
        fac = sqrt(-2.0 * log(rsq) / rsq);
        gset = v1*fac;
        iset = 1;
        return (v2 * fac);
    }
    else
    {
        iset = 0;
        return (gset);
    }
}



int myBinarySearch (vector<double> mySortedVector, double value, int exact)
{
    /* Binary search
       Exact = 1 return false if value not found
       Exact = 2 return the upper bound (false if does not exists)
       Exact = 3 return the lower bound (false if does not exists)
    */
    int low = 0;
    int high = mySortedVector.size();
    int mid = 0;

    while( low < high )
    {
        mid = floor((low+high)/2);
        if(mySortedVector[mid] > value)
        {
            high = mid - 1;
        }
        else if (mySortedVector[mid] < value)
        {
            low = mid + 1;
        }
        else return mid;
    }

    if (exact == 1) return -1; //If equality, fuction should have already ended
    else if (exact == 2)
    {
        if (low >= (int)mySortedVector.size()) return -1; //Value beyond the upper bound of the table
        return low;
    }
    else if (exact == 3)
    {
        if (high < 1) return -1;
        return high;
    }
    else printf("The \"exact\" parameter must take value 1, 2, or 3");

    return -1;
}

//Normalize values between 0 and 1
void set_norm(vector<double> &values)
{
    double mini = values[0];
    double maxi = values[0];
    for(unsigned int i = 0; i<values.size(); i++)
    {
        if(mini > values[i]) mini = values[i];
        if(maxi < values[i]) maxi = values[i];
    }

    for(unsigned int i = 0; i<values.size(); i++)
    {
        values[i] = (values[i]-mini)/(maxi-mini);
    }
}

void exit_error(string error)
{
    cout << "\nERROR:\n" << error << "\nProgram will now exit." << endl;
    exit(0);
}


void fail(string error)
{
    cout << error << endl;
    exit(EXIT_FAILURE);
}

/******************** Move in stat.cpp ****************************************/
//Make inline?
// Straightforward mean calcuation for a vector of double
double calc_mean(vector<double> data)
{
    double sum = std::accumulate(std::begin(data), std::end(data), 0.0);
    return sum / data.size();
}

//Variance of a finite pop
double calc_var(vector<double> data)
{
    double m = calc_mean(data);
    double temp = 0.0;
    std::for_each (std::begin(data), std::end(data), [&](const double d) 
        {
            temp += (d - m) * (d - m);
        });
    return temp / data.size();
}

double calc_std(vector<double> data)
{
    return sqrt(calc_var(data));
}
