pd.if.org Git - afpopgen/blob - sim/variate.c

   1 #include <stdint.h>
   2 #include <stdio.h>
   3 #include <unistd.h>
   4
   5 #include <gsl/gsl_sf_bessel.h>
   6 #include <gsl/gsl_rng.h>
   7 #include <gsl/gsl_randist.h>
   8
   9 static gsl_rng *rng;  /* global generator */
  10
  11 void var_seed(uint64_t seed) {
  12         gsl_rng_set(rng, seed);
  13 }
  14
  15 void var_init(uint64_t seed) {
  16         const gsl_rng_type * T;
  17
  18         gsl_rng_env_setup();
  19         T = gsl_rng_default;
  20         rng = gsl_rng_alloc (T);
  21
  22         if (seed == 0) {
  23                 seed = (unsigned long)getpid();
  24         }
  25
  26         var_seed(seed);
  27 }
  28
  29 unsigned int var_rhyper(unsigned int t, unsigned int n1, unsigned int n2) {
  30         return gsl_ran_hypergeometric(rng, n1, n2, t);
  31 }
  32
  33 /* must be two ints in each vector */
  34 unsigned int var_rhyperv(unsigned int t, unsigned int *s, unsigned int *d) {
  35         int x;
  36         x = var_rhyper(t, s[0], s[1]);
  37         d[0] = x; d[1] = t - x;
  38         return x;
  39 }
  40
  41 void var_rbinom(unsigned int t, unsigned int *n, unsigned int *d) {
  42         d[0] = gsl_ran_binomial(rng, (double)n[0]/((double)n[0]+(double)n[1]), t);
  43         d[1] = t - d[0];
  44 }
  45
  46 /* TODO hardcoded to size of arrays == 2 */
  47 void var_rmultinom(unsigned int t, unsigned int *n, unsigned int *d) {
  48         double p[2];
  49         size_t k = 2;
  50
  51         p[0] = n[0]; p[1] = n[1];
  52         gsl_ran_multinomial(rng, k, t, p, d);
  53 }
  54
  55 /* pick one given n weighted probabilities */
  56 int var_pick(unsigned int n, unsigned int *p) {
  57         long sum;
  58         int i;
  59         unsigned long r;
  60
  61         sum = 0;
  62         for (i=0;i<n;i++) {
  63                 sum += p[i];
  64         }
  65         r = gsl_rng_uniform_int(rng, sum);
  66
  67         sum = 0;
  68         for (i=0;i<n;i++) {
  69                 sum += p[i];
  70                 if (sum >= r) return i;
  71         }
  72         return i;
  73 }