/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * This program is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details.

 * You should have received a copy of the GNU General Public License along with
 * this program. If not, see <https://www.gnu.org/licenses/>.
 */

#include "random.h"
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "vector.h"

static gsl_rng *r;

/* Generates a random number from a normal distribution using the Box-Muller
 * transform. */
double randn(void)
{
    double u1, u2;

    u1 = genrand_real1();
    u2 = genrand_real1();

    return sqrt(-2*log(u1))*cos(2*M_PI*u2);
}

/* Generates a random point in the unit sphere. */
void rand_ball(double *pos, double radius)
{
    pos[0] = (genrand_real2()*2 - 1)*radius;
    pos[1] = (genrand_real2()*2 - 1)*radius;
    pos[2] = (genrand_real2()*2 - 1)*radius;

    while (NORM(pos) >= radius) {
        pos[0] = (genrand_real2()*2 - 1)*radius;
        pos[1] = (genrand_real2()*2 - 1)*radius;
        pos[2] = (genrand_real2()*2 - 1)*radius;
    }
}

/* Generates a random point on the unit sphere. */
void rand_sphere(double *dir)
{
    double u, v, theta, phi;

    u = genrand_real1();
    v = genrand_real1();

    phi = 2*M_PI*u;
    theta = acos(2*v-1);

    dir[0] = sin(theta)*cos(phi);
    dir[1] = sin(theta)*sin(phi);
    dir[2] = cos(theta);
}

/* Randomly choose `k` elements from array `src` with `n` elements with
 * probabilities proportional to `w`. The elements are sampled without
 * replacement and placed into the `dest` array. */
void ran_choose_weighted(int *dest, double *w, size_t k, int *src, size_t n)
{
    int i, j;
    gsl_ran_discrete_t *g;
    size_t index, result;
    int duplicate;

    if (!r)
        r = gsl_rng_alloc(gsl_rng_default);

    g = gsl_ran_discrete_preproc(n, w);

    i = 0;
    while (i < k) {
        index = gsl_ran_discrete(r, g);

        result = src[index];

        duplicate = 0;
        for (j = 0; j < i; j++) {
            if (result == dest[j]) {
                duplicate = 1;
                break;
            }
        }

        if (!duplicate) dest[i++] = result;
    }

    gsl_ran_discrete_free(g);
}