core/random.c

Go to the documentation of this file.

#include <stdlib.h>
#include <math.h>
#include <errno.h>
#include <string.h>
#include <stdarg.h>
#include <time.h>
#include <float.h>
#include "random.h"
#include "find.h"
#include "output.h"
#include "exception.h"

#ifdef WIN32
#define finite _finite
#endif

#define QNAN sqrt(-1) /* quiet NaN used only for low-key failures */

int random_init(void)
{
    /* randomizes the random number generator start value */
    srand((unsigned int)time(NULL));
    return 1;
}

RANDOMTYPE random_type(char *name) 
{
    struct {
        char *name;
        RANDOMTYPE type;
    } *p, map[] = {
        {"degenerate", RT_DEGENERATE},
        {"uniform",RT_UNIFORM},
        {"normal",RT_NORMAL},
        {"bernoulli",RT_BERNOULLI},
        {"sampled",RT_SAMPLED},
        {"pareto",RT_PARETO},
        {"lognormal",RT_LOGNORMAL},
        {"exponential",RT_EXPONENTIAL},
    };
    for (p=map; p<map+sizeof(map)/sizeof(map[0]); p++)
    {
        if (strcmp(p->name,name)==0)
            return p->type;
    }
    return RT_INVALID;
}

/* uniform distribution in range (0,1( */
double randunit(void)
{
    double u = rand()/(RAND_MAX+1.0);
    if (u<1) return u;
    return randunit();
}

double random_degenerate(double a)
{
    /* returns a, i.e., Dirac delta function */
    double aa = fabs(a);
    if (a!=0 && (aa<1e-30 || aa>1e30))
        output_warning("random_degenerate(a=%g): a is outside normal bounds of +/-1e(+/-30)", a);
    return a;
}

double random_uniform(double a, 
                      double b) 
{
    /* uniform distribution in range (a,b( */
    double aa=fabs(a), ab=fabs(b);
    if (a!=0 && (aa<1e-30 || aa>1e30))
        output_warning("random_uniform(a=%g, b=%g): a is outside normal bounds of +/-1e(+/-30)", a, b);
    if (b!=0 && (ab<1e-30 || ab>1e30))
        output_warning("random_uniform(a=%g, b=%g): b is outside normal bounds of +/-1e(+/-30)", a, b);
    if (b<a)
        output_warning("random_uniform(a=%g, b=%g): b is less than a", a, b);
    return randunit()*(b-a)+a;
}

double random_normal(double m, 
                     double s) 
{
    /* normal distribution centered on m, with variance s^2 */
    double r=randunit();
    if (s<0)
        output_warning("random_normal(m=%g, s=%g): s is negative", m, s);
    while (r<=0 || r>1)
        r=randunit();
    return sqrt(-2*log(r)) * sin(2*PI*randunit())*s+m;
}

double random_bernoulli(double p) 
{
    double ap = fabs(p);
    if (ap!=0 && (ap<1e-30 || ap>1e30))
        output_warning("random_bernoulli(p=%g): p is not within the normal bounds of +/-1e(+/-30)", p);
    if (p<0 || p>1)
        output_warning("random_bernoulli(p=%g): p is not between 0 and 1", p);
    /* 1 if rand<=p, 0 if rand>p */
    return (p>=randunit()) ? 1 : 0;
}

double random_sampled(unsigned n, 
                      double *x) 
{
    if (n>0)
    {
        double v = x[(unsigned)(randunit()*n)];
        double av = fabs(v);
        if (v!=0 && (v<1e-30 || v>1e30))
            output_warning("random_sampled(n=%d,...): sampled value is not within normal bounds of +/-1e(+/-30)",n);
        return v;
    }
    else
    {
        throw_exception("random_sampled(n=%d,...): n must be a positive number", n);
        return QNAN; /* never gets here */
    }
}


double random_pareto(double m, 
                     double k) 
{
    double am = fabs(m);
    double ak = fabs(k);
    double r = randunit();
    if (am!=0 && (am<1e-03 || am>1e30))
        output_warning("random_pareto(m=%g, k=%g): m is not within the normal bounds of +/-1e(+/-30)", m, k);
    if (ak>1e30)
        output_warning("random_pareto(m=%g, k=%g): k is not within the normal bounds of +/-1e(+/-30)", m, k);
    if (k<=0)
        throw_exception("random_pareto(m=%g, k=%g): k must be greater than 1", m, k);
    while (r<=0 || r>=1)
        r = randunit();
    return m*pow(r,-1/k);
}

double random_lognormal(double gmu, 
                        double gsigma) 
{
    return exp(random_normal(0,1)*gsigma+gmu);
}

double random_exponential(double lambda) 
{
    double r=randunit();
    if (lambda<=0)
        throw_exception("random_exponential(l=%g): l must be greater than 0", lambda);
    if (lambda<1e-30 || lambda>1e30)
        output_warning("random_exponential(l=%g): l is not within the normal bounds of 1e(+/-30)", lambda);
    while (r<=0 || r>=1)
        r=randunit();
    return -log(r)/lambda;
}

/* internal function that generates a random number */
static double _random_value(RANDOMTYPE type, va_list ptr)
{
    switch (type) {
    case RT_DEGENERATE:/* ... double value */
        return random_degenerate(va_arg(ptr,double));
    case RT_UNIFORM:        /* ... double min, double max */
        {   double min = va_arg(ptr,double);
            double max = va_arg(ptr,double);
            return random_normal(min,max);
        }
    case RT_NORMAL:     /* ... double mean, double stdev */
        {   double mu = va_arg(ptr,double);
            double sigma = va_arg(ptr,double);
            return random_normal(mu,sigma);
        }
    case RT_BERNOULLI:  /* ... double p */
        return random_bernoulli(va_arg(ptr,double));
    case RT_SAMPLED: /* ... unsigned n_samples, double samples[n_samples] */
        {   unsigned n_samples = va_arg(ptr,unsigned);
            double *samples = va_arg(ptr,double*);
            return random_sampled(n_samples,samples);
        }
    case RT_PARETO: /* ... double base, double gamma */
        {   double base = va_arg(ptr,double);
            double gamma = va_arg(ptr,double);
            return random_pareto(base,gamma);
        }
    case RT_LOGNORMAL:  /* ... double gmean, double gsigma */
        {   double gmu = va_arg(ptr,double);
            double gsigma = va_arg(ptr,double);
            return random_lognormal(gmu,gsigma);
        }
    case RT_EXPONENTIAL: /* ... */
        {   double lambda = va_arg(ptr,double);
            return random_exponential(lambda);
        }
    default:
        throw_exception("_random_value(type=%d,...); type is not valid",type);
    }
    return QNAN; /* never gets here */
}

int random_apply(char *group_expression, 
                 char *property, 
                 RANDOMTYPE type, 
                 ...) 
{
    FINDLIST *list = find_objects(FL_GROUP, group_expression);
    OBJECT *obj;
    unsigned count=0;
    va_list ptr;
    va_start(ptr,type);
    for (obj=find_first(list); obj!=NULL; find_next(list,obj))
    {
        /* this is quite slow and should use a class property lookup */
        object_set_double_by_name(obj,property,_random_value(type,ptr));
        count++;
    }
    va_end(ptr);
    return count;
}

double random_value(RANDOMTYPE type, 
                    ...) 
{
    double x;
    va_list ptr;
    va_start(ptr,type);
    x = _random_value(type,ptr);
    va_end(ptr);
    return x;
}

/******************************************************************************/
static double mean(double sample[], unsigned int count)
{
    double sum=0;
    unsigned int i;
    for (i=0; i<count; i++)
        sum += sample[i];
    return sum/count;
}
#ifdef min
#undef min
#endif
static double min(double sample[], unsigned int count)
{
    double min;
    unsigned int i;
    for (i=0; i<count; i++)
    {
        if (i==0) min=sample[0];
        else if (sample[i]<min) min=sample[i];
    }
    return min;
}
#ifdef max
#undef max
#endif
static double max(double sample[], unsigned int count)
{
    double max;
    unsigned int i;
    for (i=0; i<count; i++)
    {
        if (i==0) max=sample[0];
        else if (sample[i]>max) max=sample[i];
    }
    return max;
}
static double stdev(double sample[], unsigned int count)
{
    double sum=0, mean=0;
    unsigned int n = 0, i;
    for (i=0; i<count; i++)
    {
        double delta = sample[i] - mean;
        n++;
        mean += delta/n;
        sum += delta*(sample[i]-mean);
    }
    return sqrt(sum/(n-1));
}
static void sort(double sample[], unsigned int count)
{
}
static void report(char *parameter, double actual, double expected)
{
    if (parameter==NULL)
    {
        output_test("   Parameter       Actual    Expected    Error");
        output_test("---------------- ---------- ---------- ----------");
    }
    else
    {
        if (expected==0)
            output_test("%-16.16s %10.4f %10.4f %9.6f ", parameter,actual,expected,actual-expected);
        else
            output_test("%-16.16s %10.4f %10.4f %9.6f%%", parameter,actual,expected,(actual-expected)/expected*100);
    }
}

int random_test(void)
{
    static double sample[1000000];
    double a, b;
    unsigned int count = sizeof(sample)/sizeof(sample[0]);
    unsigned int i;

    /* Dirac distribution test */
    a = 10*randunit()/2-5;
    output_test("\ndegenerate(x=%g)",a);
    for (i=0; i<count; i++)
    {
        sample[i] = random_degenerate(a);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),a);
    report("Stdev",stdev(sample,count),0);
    report("Min",min(sample,count),a);
    report("Max",max(sample,count),a);

    /* uniform distribution test */
    a = 10*randunit()/2;
    b = 10*randunit()/2 + 5;
    output_test("\nuniform(min=%g, max=%g)",a,b);
    for (i=0; i<count; i++)
    {
        sample[i] = random_uniform(a,b);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),(a+b)/2);
    report("Stdev",stdev(sample,count),sqrt((b-a)*(b-a)/12));
    report("Min",min(sample,count),a);
    report("Max",max(sample,count),b);

    /* Bernoulli distribution test */
    a = randunit();
    output_test("\nBernoulli(prob=%g)",a);
    for (i=0; i<count; i++)
    {
        sample[i] = random_bernoulli(a);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),a);
    report("Stdev",stdev(sample,count),sqrt(a*(1-a)));
    report("Min",min(sample,count),0);
    report("Max",max(sample,count),1);

    /* normal distribution test */
    a = 20*randunit()-5;
    b = 5*randunit();
    output_test("\nnormal(mean=%g, stdev=%g)",a,b);
    for (i=0; i<count; i++)
    {
        sample[i] = random_normal(a,b);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),a);
    report("Stdev",stdev(sample,count),b);
    
    /* exponential distribution test */
    a = 1/randunit()-1;
    output_test("\nexponential(lambda=%g)",a);
    for (i=0; i<count; i++)
    {
        sample[i] = random_exponential(a);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),1/a);
    report("Stdev",stdev(sample,count),1/a);
    report("Min",min(sample,count),0);
    
    /* lognormal distribution test */
    a = 2*randunit()-1;
    b = 2*randunit();
    output_test("\nlognormal(gmean=%g, gstdev=%g)",a,b);
    for (i=0; i<count; i++)
    {
        sample[i] = random_lognormal(a,b);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),exp(a+b*b/2));
    report("Stdev",stdev(sample,count),sqrt((exp(b*b)-1)*exp(2*a+b*b)));
    report("Min",min(sample,count),0);

    /* Pareto distribution test */
    a = 10*randunit();
    b = randunit()*2+2;
    output_test("\nPareto(base=%g, gamma=%g)",a,b);
    for (i=0; i<count; i++)
    {
        sample[i] = random_pareto(a,b);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),b*a/(b-1));
    report("Stdev",stdev(sample,count),sqrt(a*a*b/((b-1)*(b-1)*(b-2))));
    report("Min",min(sample,count),a);

    /* sampled distribution test */
    output_test("\nSampled(count=10, sample=[0..9])");
    for (i=0; i<count; i++)
    {
        double set[10]={0,1,2,3,4,5,6,7,8,9};
        sample[i] = random_sampled(10,set);
        if (!finite(sample[i]))
            output_test("Sample %d is not a finite number!",i--);
    }
    report(NULL,0,0);
    report("Mean",mean(sample,count),4.5);
    //report("Stdev",stdev(sample,count),sqrt(9*9/12));
    report("Stdev",stdev(sample,count),sqrt(99/12)); /* sqrt((b-a+1)^2-1 / 12)*/ /* 2.87 is more accurate and was Mathematica's answer */
    report("Min",min(sample,count),0);
    report("Max",max(sample,count),9);

    output_warning("random_test() not implemented yet");
    return 0;
}

---

GridLAB-D^TM Version 1.0

An open-source project initiated by the US Department of Energy