/** @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.8 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-06-01 23:51:25 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

#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 "psMemory.h"
#include "psRandom.h"
#include "psScalar.h"
#include "psError.h"
#include "psTrace.h"
#include "psLogMsg.h"
#include "psConstants.h"
#include "psDataManipErrors.h"

psU64 p_psRandomGetSystemSeed()
{
    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
        fread(&seedVal, sizeof(psU64),1,fd);
        // Close file
        fclose(fd);
    }

    // Send log message of the system seed value used
    psLogMsg(__func__,PS_LOG_INFO,"System random seed value used  seed = %llX hex",seedVal);

    return seedVal;
}

psRandom *psRandomAlloc(psRandomType type,
                        psU64 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());
        } else {
            gsl_rng_set(myRNG->gsl,seed);
        }
        myRNG->type = type;
        break;

    default:
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                PS_ERRORTEXT_psRandom_UNKNOWN_RANDFOM_NUMBER_GENERATOR_TYPE);
        break;
    }

    return(myRNG);
}

void psRandomReset(psRandom *rand,
                   psU64 seed)
{
    // Check null psRandom
    if(rand==NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                PS_ERRORTEXT_psRandom_NULL_RANDOM_VAR);
    } else {
        // Check seed value to see if system seed should be used
        if(seed == 0) {
            gsl_rng_set(rand->gsl,p_psRandomGetSystemSeed());
        } else {
            gsl_rng_set(rand->gsl, seed);
        }
    }
}

psF64 psRandomUniform(const psRandom *r)
{
    // Check null psRandom variable
    if(r == NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                true,
                PS_ERRORTEXT_psRandom_NULL_RANDOM_VAR);
        return(0);
    } else {
        return(gsl_rng_uniform(r->gsl));
    }
}

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

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

