src/random.c

/*
            __
           / /_______________  ____  ___
          / __/ ___/ ___/ __ \/ __ \/ _ \
         / /_(__  ) /__/ /_/ / /_/ /  __/
         \__/____/\___/\____/ .___/\___/
                           /_/
    
    some macros for backward compatibility

    By Michael Stevens

    See license.html for copyright information
*/

#include "../include/tscope.h"
#include "../include/tscope/internal.h"

//------------------------------------------------------------------------
// random seed functions


// get seeds from user
void ts_rseed(int userseed[3])
{
    _rndseed[0] = userseed[0];
    _rndseed[1] = userseed[1];
    _rndseed[2] = userseed[2];
    if (_debug_flag > DEBUG3)
        fprintf(stderr, "User defined random seeds: %d, %d and %d.\n",
                _rndseed[0], _rndseed[1], _rndseed[2]);
    if (!_random_flag) {
        atexit(_random_exit);
        _random_flag = 1;
    }
}

// read seeds from file
void ts_rseedfile(char *file)
{
    FILE *seedfile;
    if ((seedfile = fopen(file, "r")) == NULL)
        ts_fatal("ts_rseedfile: error reading seed file.");
    else {
        fscanf(seedfile, "%d", &_rndseed[0]);
        fscanf(seedfile, "%d", &_rndseed[1]);
        fscanf(seedfile, "%d", &_rndseed[2]);
        fclose(seedfile);
        if (_debug_flag > DEBUG3)
            fprintf(stderr, "Random seeds from file %s: %d, %d and %d.\n",
                    file, _rndseed[0], _rndseed[1], _rndseed[2]);
    }
    if (!_random_flag) {
        atexit(_random_exit);
        _random_flag = 1;
    }
}

//------------------------------------------------------------------------
// basic randomisation function
// all other random functions call this one

// double float [0:1] from uniform distribution
double ts_runif()
{
    double ran;

    // do we have a random seed yet?
    if (!_random_flag)
        _random_init();

    // generate 3 new random seeds with % and /
    _rndseed[0] = 171 * (_rndseed[0] % 177) - 2 * (_rndseed[0] / 177);
    _rndseed[1] = 172 * (_rndseed[1] % 176) - 35 * (_rndseed[1] / 176);
    _rndseed[2] = 170 * (_rndseed[2] % 178) - 63 * (_rndseed[2] / 178);

    // _rndseed must be > 0
    if (_rndseed[0] < 0)
        _rndseed[0] = _rndseed[0] + 30269;
    if (_rndseed[1] < 0)
        _rndseed[1] = _rndseed[1] + 30307;
    if (_rndseed[2] < 0)
        _rndseed[2] = _rndseed[2] + 30323;

    // generate random number with the random seeds
    ran = ((float) _rndseed[0] / 30269
           + (float) _rndseed[1] / 30307 + (float) _rndseed[2] / 30323);

    // throw away integer part, rescale
    ran = ran - (int) ran;
    if (ran <= 0)
        ran = 1e-30;
    if (ran >= 1)
        ran = 0.9999999999;

    return ran;
}

//------------------------------------------------------------------------
// derived random functions: other distributions

// double float from exponential distribution
double ts_rexp()
{
    //take sample from uniform distribution and take the log
    return -log(ts_runif());
}

// double float from normal distribution
double ts_rnorm(double mean, double sd)
{
    double u1, u2, a, z;

    //take sample from uniform distribution and convert
    u1 = ts_runif();

    if (u1 > 0.5)
        u2 = 1 - u1;
    else
        u2 = u1;

    if (u2 < 1e-20)
        z = 10;
    else {
        a = sqrt(-2 * log(u2));
        z = a - ((7.45551 * a + 450.636) * a + 1271.059)
            / (((a + 110.4212) * a + 750.365) * a + 500.756);
    }

    if (u1 > 0.5)
        z = -z;

    if (sd < 0)
        sd = -sd;

    return mean + (z * sd);
}

//------------------------------------------------------------------------
// derived random functions: random integers / integer lists

// integer from uniform distribution [0:nmax-1]
int ts_rint(int nmax)
{
    return (int) (ts_runif() * nmax);
}

// random list without replacement
int ts_rlist(int nmax, int freq, int *list)
{
    int i;
    int rnd, tmp;

    // fill the list
    for (i = 0; i < nmax * freq; i++)
        list[i] = i / freq;

    // rearrange randomly
    for (i = (nmax * freq) - 1; i > 0; i--) {
        rnd = ts_rint(i + 1);
        tmp = list[i];
        list[i] = list[rnd];
        list[rnd] = tmp;
    }
    return 1;
}

// controlled randomisation without replacement
// every element follows every other element the same number of times
int ts_rslist(int nmax, int freq, int *list)
{
    int i, j, tmp;
    int digit_counter[nmax];
    int prev_counter[nmax][nmax];
    int total = (nmax * nmax * freq) + 1;
    int possible = 0;

    // reset counters
    for (i = 0; i < nmax; i++) {
        digit_counter[i] = 0;
        for (j = 0; j < nmax; j++)
            prev_counter[i][j] = 0;
    }

    list[0] = ts_rint(nmax);
    digit_counter[list[0]]++;

    // randomize the list
    for (i = 1; i < total - 1; i++) {

        // draw a possible candidate
        tmp = ts_rint(nmax);

        //every number : max n_max*n_comb times
        if (digit_counter[tmp] >= nmax * freq) {
            i--;
            continue;
        }
        //combinations : max n_comb times
        if (prev_counter[tmp][list[i - 1]] >= freq) {
            i--;
            continue;
        }
        //if all test are ok -> use the number
        list[i] = tmp;
        digit_counter[tmp]++;
        prev_counter[tmp][list[i - 1]]++;

        //test whether it will be possible to draw another number
        if (i < total - 2)      //no test for the last number
        {
            possible = 0;
            for (j = 0; j < nmax; j++) {
                if (j == tmp)
                    continue;

                if (digit_counter[j] < nmax * freq)
                    if (prev_counter[j][tmp] < freq)
                        possible = 1;
            }
        }
        if (possible == 0)
            break;
    }
    list[total - 1] = list[0];

    //restart if it is impossible to draw another number
    if (possible == 0)
        ts_rslist(nmax, freq, list);

    return 1;
}
top