/** @file psRandom.c
*  \brief Random Number Generators
*  \ingroup Math
*
*  This file will hold the functions which allocate, free,
*  and evaluate random number Generators.
*
*  @ingroup Math
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
*  @date $Date: 2008-11-05 11:12:39 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <stdio.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <inttypes.h>

#include "psMemory.h"
#include "psRandom.h"
#include "psScalar.h"
#include "psError.h"
#include "psTrace.h"
#include "psLogMsg.h"
#include "psAssert.h"


psU64 p_psRandomGetSystemSeed(bool log)
{
    FILE*  fd;
    psU64  seedVal = 0;
    time_t timeVal;

    fd = fopen("/dev/urandom","r");
    if(fd == NULL) {
        // Read system clock to get seed
        seedVal = (psU64)time(&timeVal);
    } else {
        // Read urandom to get seed
        if (fread(&seedVal, sizeof(psU64),1,fd)) {;} // ignore return value
        // Close file
        fclose(fd);
    }

    // Send log message of the system seed value used
    if (log) {
        psLogMsg(__func__,PS_LOG_INFO,"System random seed value used = %" PRIx64, seedVal);
    }

    return seedVal;
}

psRandom *psRandomAlloc(psRandomType type,
                        unsigned long seed)
{
    gsl_rng   *r      = NULL;
    psRandom  *myRNG  = NULL;

    switch (type) {
    case PS_RANDOM_TAUS:
        myRNG = (psRandom*)psAlloc(sizeof(psRandom));
        r = gsl_rng_alloc(gsl_rng_taus);
        myRNG->gsl = r;
        if(seed == 0) {
            gsl_rng_set(myRNG->gsl, p_psRandomGetSystemSeed(true));
        } else {
            gsl_rng_set(myRNG->gsl, seed);
        }
        myRNG->type = type;
        break;

    default:
        psError(PS_ERR_UNEXPECTED_NULL, true, _("Unknown Random Number Generator Type"));
        break;
    }

    return(myRNG);
}

void psRandomReset(psRandom *rand,
                   unsigned long seed)
{
    // Check null psRandom
    if(rand==NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                _("Random variable is NULL."));
    } else {
        // Check seed value to see if system seed should be used
        if(seed == 0) {
            gsl_rng_set(rand->gsl,p_psRandomGetSystemSeed(true));
        } else {
            gsl_rng_set(rand->gsl, seed);
        }
    }
}

double psRandomUniform(const psRandom *r)
{
    // Check null psRandom variable
    if(r == NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                _("Random variable is NULL."));
        return(0);
    } else {
        return(gsl_rng_uniform(r->gsl));
    }
}

double psRandomGaussian(const psRandom *r)
{
    // Check null psRandom variable
    if(r == NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                _("Random variable is NULL."));
        return(0);
    } else {
        // XXX: What should sigma be?
        return(gsl_ran_gaussian(r->gsl, 1.0));
    }
}

double p_psRandomGaussian(const psRandom *r, double sigma)
{
    // Check null psRandom variable
    if(r == NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                _("Random variable is NULL."));
        return(0);
    } else {
        return(gsl_ran_gaussian(r->gsl, sigma));
    }
}

double psRandomPoisson(const psRandom *r, double mean)
{
    // Check null psRandom variable
    if(r == NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                _("Random variable is NULL."));
        return(0);
    } else {
        return((psF64) gsl_ran_poisson(r->gsl, mean));
    }
}

