text:
/* randist/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2010 James Theiler, Brian Gough
 *
 * 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cdf.h>

#define N 100000

/* Convient test dimension for multivariant distributions */
#define MULTI_DIM 10

void testMoments (double (*f) (void), const char *name,
                  double a, double b, double p);
void testPDF (double (*f) (void), double (*pdf) (double), const char *name);
void testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
                      const char *name);

void test_shuffle (void);
void test_choose (void);
double test_beta (void);
double test_beta_pdf (double x);
double test_bernoulli (void);
double test_bernoulli_pdf (unsigned int n);

double test_binomial (void);
double test_binomial_pdf (unsigned int n);
double test_binomial_large (void);
double test_binomial_large_pdf (unsigned int n);
double test_binomial_huge (void);
double test_binomial_huge_pdf (unsigned int n);
double test_binomial_max (void);
double test_binomial_max_pdf (unsigned int n);
double test_binomial0 (void);
double test_binomial0_pdf (unsigned int n);
double test_binomial1 (void);
double test_binomial1_pdf (unsigned int n);

double test_binomial_knuth (void);
double test_binomial_knuth_pdf (unsigned int n);
double test_binomial_large_knuth (void);
double test_binomial_large_knuth_pdf (unsigned int n);
double test_binomial_huge_knuth (void);
double test_binomial_huge_knuth_pdf (unsigned int n);

double test_cauchy (void);
double test_cauchy_pdf (double x);
double test_chisq (void);
double test_chisq_pdf (double x);
double test_chisqnu2 (void);
double test_chisqnu2_pdf (double x);
double test_dirichlet (void);
double test_dirichlet_pdf (double x);
double test_dirichlet_small (void);
double test_dirichlet_small_pdf (double x);
void test_dirichlet_moments (void);
double test_discrete1 (void);
double test_discrete1_pdf (unsigned int n);
double test_discrete2 (void);
double test_discrete2_pdf (unsigned int n);
double test_discrete3 (void);
double test_discrete3_pdf (unsigned int n);
double test_erlang (void);
double test_erlang_pdf (double x);
double test_exponential (void);
double test_exponential_pdf (double x);
double test_exppow0 (void);
double test_exppow0_pdf (double x);
double test_exppow1 (void);
double test_exppow1_pdf (double x);
double test_exppow1a (void);
double test_exppow1a_pdf (double x);
double test_exppow2 (void);
double test_exppow2_pdf (double x);
double test_exppow2a (void);
double test_exppow2a_pdf (double x);
double test_exppow2b (void);
double test_exppow2b_pdf (double x);
double test_fdist (void);
double test_fdist_pdf (double x);
double test_fdist_large (void);
double test_fdist_large_pdf (double x);
double test_flat (void);
double test_flat_pdf (double x);
double test_gamma (void);
double test_gamma_pdf (double x);
double test_gamma1 (void);
double test_gamma1_pdf (double x);
double test_gamma_int (void);
double test_gamma_int_pdf (double x);
double test_gamma_large (void);
double test_gamma_large_pdf (double x);
double test_gamma_vlarge (void);
double test_gamma_vlarge_pdf (double x);
double test_gamma_small (void);
double test_gamma_small_pdf (double x);
double test_gamma_mt (void);
double test_gamma_mt_pdf (double x);
double test_gamma_mt1 (void);
double test_gamma_mt1_pdf (double x);
double test_gamma_mt_int (void);
double test_gamma_mt_int_pdf (double x);
double test_gamma_mt_large (void);
double test_gamma_mt_large_pdf (double x);
double test_gamma_mt_small (void);
double test_gamma_mt_small_pdf (double x);
double test_gamma_knuth_vlarge (void);
double test_gamma_knuth_vlarge_pdf (double x);
double test_gaussian (void);
double test_gaussian_pdf (double x);
double test_gaussian_ratio_method (void);
double test_gaussian_ratio_method_pdf (double x);
double test_gaussian_ziggurat (void);
double test_gaussian_ziggurat_pdf (double x);
double test_gaussian_tail (void);
double test_gaussian_tail_pdf (double x);
double test_gaussian_tail1 (void);
double test_gaussian_tail1_pdf (double x);
double test_gaussian_tail2 (void);
double test_gaussian_tail2_pdf (double x);
double test_ugaussian (void);
double test_ugaussian_pdf (double x);
double test_ugaussian_ratio_method (void);
double test_ugaussian_ratio_method_pdf (double x);
double test_ugaussian_tail (void);
double test_ugaussian_tail_pdf (double x);
double test_bivariate_gaussian1 (void);
double test_bivariate_gaussian1_pdf (double x);
double test_bivariate_gaussian2 (void);
double test_bivariate_gaussian2_pdf (double x);
double test_bivariate_gaussian3 (void);
double test_bivariate_gaussian3_pdf (double x);
double test_bivariate_gaussian4 (void);
double test_bivariate_gaussian4_pdf (double x);
void test_multivariate_gaussian_log_pdf (void);
void test_multivariate_gaussian_pdf (void);
void test_multivariate_gaussian (void);
double test_gumbel1 (void);
double test_gumbel1_pdf (double x);
double test_gumbel2 (void);
double test_gumbel2_pdf (double x);
double test_geometric (void);
double test_geometric_pdf (unsigned int x);
double test_geometric1 (void);
double test_geometric1_pdf (unsigned int x);
double test_hypergeometric1 (void);
double test_hypergeometric1_pdf (unsigned int x);
double test_hypergeometric2 (void);
double test_hypergeometric2_pdf (unsigned int x);
double test_hypergeometric3 (void);
double test_hypergeometric3_pdf (unsigned int x);
double test_hypergeometric4 (void);
double test_hypergeometric4_pdf (unsigned int x);
double test_hypergeometric5 (void);
double test_hypergeometric5_pdf (unsigned int x);
double test_hypergeometric6 (void);
double test_hypergeometric6_pdf (unsigned int x);
double test_landau (void);
double test_landau_pdf (double x);
double test_levy1 (void);
double test_levy1_pdf (double x);
double test_levy2 (void);
double test_levy2_pdf (double x);
double test_levy1a (void);
double test_levy1a_pdf (double x);
double test_levy2a (void);
double test_levy2a_pdf (double x);
double test_levy_skew1 (void);
double test_levy_skew1_pdf (double x);
double test_levy_skew2 (void);
double test_levy_skew2_pdf (double x);
double test_levy_skew1a (void);
double test_levy_skew1a_pdf (double x);
double test_levy_skew2a (void);
double test_levy_skew2a_pdf (double x);
double test_levy_skew1b (void);
double test_levy_skew1b_pdf (double x);
double test_levy_skew2b (void);
double test_levy_skew2b_pdf (double x);
double test_logistic (void);
double test_logistic_pdf (double x);
double test_lognormal (void);
double test_lognormal_pdf (double x);
double test_logarithmic (void);
double test_logarithmic_pdf (unsigned int n);
double test_multinomial (void);
double test_multinomial_pdf (unsigned int n);
double test_multinomial_large (void);
double test_multinomial_large_pdf (unsigned int n);
void test_multinomial_moments (void);
double test_negative_binomial (void);
double test_negative_binomial_pdf (unsigned int n);
double test_pascal (void);
double test_pascal_pdf (unsigned int n);
double test_pareto (void);
double test_pareto_pdf (double x);
double test_poisson (void);
double test_poisson_pdf (unsigned int x);
double test_poisson_large (void);
double test_poisson_large_pdf (unsigned int x);
double test_dir2d (void);
double test_dir2d_pdf (double x);
double test_dir2d_trig_method (void);
double test_dir2d_trig_method_pdf (double x);
double test_dir3dxy (void);
double test_dir3dxy_pdf (double x);
double test_dir3dyz (void);
double test_dir3dyz_pdf (double x);
double test_dir3dzx (void);
double test_dir3dzx_pdf (double x);
double test_rayleigh (void);
double test_rayleigh_pdf (double x);
double test_rayleigh_tail (void);
double test_rayleigh_tail_pdf (double x);
double test_tdist1 (void);
double test_tdist1_pdf (double x);
double test_tdist2 (void);
double test_tdist2_pdf (double x);
double test_laplace (void);
double test_laplace_pdf (double x);
double test_weibull (void);
double test_weibull_pdf (double x);
double test_weibull1 (void);
double test_weibull1_pdf (double x);

gsl_rng *r_global;

static gsl_ran_discrete_t *g1 = NULL;
static gsl_ran_discrete_t *g2 = NULL;
static gsl_ran_discrete_t *g3 = NULL;

int
main (void)
{
  gsl_ieee_env_setup ();

  gsl_rng_env_setup ();
  r_global = gsl_rng_alloc (gsl_rng_default);

#define FUNC(x)  test_ ## x,                     "test gsl_ran_" #x
#define FUNC2(x) test_ ## x, test_ ## x ## _pdf, "test gsl_ran_" #x

  test_shuffle ();
  test_choose ();

  testMoments (FUNC (ugaussian), 0.0, 100.0, 0.5);
  testMoments (FUNC (ugaussian), -1.0, 1.0, 0.6826895);
  testMoments (FUNC (ugaussian), 3.0, 3.5, 0.0011172689);
  testMoments (FUNC (ugaussian_tail), 3.0, 3.5, 0.0011172689 / 0.0013498981);
  testMoments (FUNC (exponential), 0.0, 1.0, 1 - exp (-0.5));
  testMoments (FUNC (cauchy), 0.0, 10000.0, 0.5);

  testMoments (FUNC (discrete1), -0.5, 0.5, 0.59);
  testMoments (FUNC (discrete1), 0.5, 1.5, 0.40);
  testMoments (FUNC (discrete1), 1.5, 3.5, 0.01);

  testMoments (FUNC (discrete2), -0.5,  0.5, 1.0/45.0 );
  testMoments (FUNC (discrete2),  8.5,  9.5, 0 );

  testMoments (FUNC (discrete3), -0.5, 0.5, 0.05 );
  testMoments (FUNC (discrete3),  0.5, 1.5, 0.05 );
  testMoments (FUNC (discrete3), -0.5, 9.5, 0.5 );

  test_dirichlet_moments ();
  test_multinomial_moments ();

  testPDF (FUNC2 (beta));
  testPDF (FUNC2 (cauchy));
  testPDF (FUNC2 (chisq));
  testPDF (FUNC2 (chisqnu2));
  testPDF (FUNC2 (dirichlet));
  testPDF (FUNC2 (dirichlet_small));
  testPDF (FUNC2 (erlang));
  testPDF (FUNC2 (exponential));

  testPDF (FUNC2 (exppow0));
  testPDF (FUNC2 (exppow1));
  testPDF (FUNC2 (exppow1a));
  testPDF (FUNC2 (exppow2));
  testPDF (FUNC2 (exppow2a));
  testPDF (FUNC2 (exppow2b));

  testPDF (FUNC2 (fdist));
  testPDF (FUNC2 (fdist_large));
  testPDF (FUNC2 (flat));
  testPDF (FUNC2 (gamma));
  testPDF (FUNC2 (gamma1));
  testPDF (FUNC2 (gamma_int));
  testPDF (FUNC2 (gamma_large));
  testPDF (FUNC2 (gamma_vlarge));
  testPDF (FUNC2 (gamma_knuth_vlarge));
  testPDF (FUNC2 (gamma_small));
  testPDF (FUNC2 (gamma_mt));
  testPDF (FUNC2 (gamma_mt1));
  testPDF (FUNC2 (gamma_mt_int));
  testPDF (FUNC2 (gamma_mt_large));
  testPDF (FUNC2 (gamma_mt_small));
  testPDF (FUNC2 (gaussian));
  testPDF (FUNC2 (gaussian_ratio_method));
  testPDF (FUNC2 (gaussian_ziggurat));
  testPDF (FUNC2 (ugaussian));
  testPDF (FUNC2 (ugaussian_ratio_method));
  testPDF (FUNC2 (gaussian_tail));
  testPDF (FUNC2 (gaussian_tail1));
  testPDF (FUNC2 (gaussian_tail2));
  testPDF (FUNC2 (ugaussian_tail));

  testPDF (FUNC2 (bivariate_gaussian1));
  testPDF (FUNC2 (bivariate_gaussian2));
  testPDF (FUNC2 (bivariate_gaussian3));
  testPDF (FUNC2 (bivariate_gaussian4));

  test_multivariate_gaussian_log_pdf ();
  test_multivariate_gaussian_pdf ();
  test_multivariate_gaussian ();

  testPDF (FUNC2 (gumbel1));
  testPDF (FUNC2 (gumbel2));
  testPDF (FUNC2 (landau));
  testPDF (FUNC2 (levy1));
  testPDF (FUNC2 (levy2));
  testPDF (FUNC2 (levy1a));
  testPDF (FUNC2 (levy2a));
  testPDF (FUNC2 (levy_skew1));
  testPDF (FUNC2 (levy_skew2));
  testPDF (FUNC2 (levy_skew1a));
  testPDF (FUNC2 (levy_skew2a));
  testPDF (FUNC2 (levy_skew1b));
  testPDF (FUNC2 (levy_skew2b));
  testPDF (FUNC2 (logistic));
  testPDF (FUNC2 (lognormal));
  testPDF (FUNC2 (pareto));
  testPDF (FUNC2 (rayleigh));
  testPDF (FUNC2 (rayleigh_tail));
  testPDF (FUNC2 (tdist1));
  testPDF (FUNC2 (tdist2));
  testPDF (FUNC2 (laplace));
  testPDF (FUNC2 (weibull));
  testPDF (FUNC2 (weibull1));

  testPDF (FUNC2 (dir2d));
  testPDF (FUNC2 (dir2d_trig_method));
  testPDF (FUNC2 (dir3dxy));
  testPDF (FUNC2 (dir3dyz));
  testPDF (FUNC2 (dir3dzx));

  testDiscretePDF (FUNC2 (discrete1));
  testDiscretePDF (FUNC2 (discrete2));
  testDiscretePDF (FUNC2 (discrete3));
  testDiscretePDF (FUNC2 (poisson));
  testDiscretePDF (FUNC2 (poisson_large));
  testDiscretePDF (FUNC2 (bernoulli));
  testDiscretePDF (FUNC2 (binomial));
  testDiscretePDF (FUNC2 (binomial0));
  testDiscretePDF (FUNC2 (binomial1));
  testDiscretePDF (FUNC2 (binomial_knuth));
  testDiscretePDF (FUNC2 (binomial_large));
  testDiscretePDF (FUNC2 (binomial_large_knuth));
  testDiscretePDF (FUNC2 (binomial_huge));
  testDiscretePDF (FUNC2 (binomial_huge_knuth));
  testDiscretePDF (FUNC2 (binomial_max));
  testDiscretePDF (FUNC2 (geometric));
  testDiscretePDF (FUNC2 (geometric1));
  testDiscretePDF (FUNC2 (hypergeometric1));
  testDiscretePDF (FUNC2 (hypergeometric2));
  testDiscretePDF (FUNC2 (hypergeometric3));
  testDiscretePDF (FUNC2 (hypergeometric4));
  testDiscretePDF (FUNC2 (hypergeometric5));
  testDiscretePDF (FUNC2 (hypergeometric6));
  testDiscretePDF (FUNC2 (logarithmic));
  testDiscretePDF (FUNC2 (multinomial));
  testDiscretePDF (FUNC2 (multinomial_large));
  testDiscretePDF (FUNC2 (negative_binomial));
  testDiscretePDF (FUNC2 (pascal));

  gsl_rng_free (r_global);
  gsl_ran_discrete_free (g1);
  gsl_ran_discrete_free (g2);
  gsl_ran_discrete_free (g3);

  exit (gsl_test_summary ());
}

void
test_shuffle (void)
{
  double count[10][10];
  int x[10] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
  int i, j, status = 0;

  for (i = 0; i < 10; i++)
    {
      for (j = 0; j < 10; j++)
        {
          count[i][j] = 0;
        }
    }

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < 10; j++)
        x[j] = j;

      gsl_ran_shuffle (r_global, x, 10, sizeof (int));

      for (j = 0; j < 10; j++)
        count[x[j]][j]++;
    }

  for (i = 0; i < 10; i++)
    {
      for (j = 0; j < 10; j++)
        {
          double expected = N / 10.0;
          double d = fabs (count[i][j] - expected);
          double sigma = d / sqrt (expected);
          if (sigma > 5 && d > 1)
            {
              status = 1;
              gsl_test (status,
                        "gsl_ran_shuffle %d,%d (%g observed vs %g expected)",
                        i, j, count[i][j] / N, 0.1);
            }
        }
    }

  gsl_test (status, "gsl_ran_shuffle on {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}");

}

void
test_choose (void)
{
  double count[10];
  int x[10] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
  int y[3] = { 0, 1, 2 };
  int i, j, status = 0;

  for (i = 0; i < 10; i++)
    {
      count[i] = 0;
    }

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < 10; j++)
        x[j] = j;

      gsl_ran_choose (r_global, y, 3, x, 10, sizeof (int));

      for (j = 0; j < 3; j++)
        count[y[j]]++;
    }

  for (i = 0; i < 10; i++)
    {
      double expected = 3.0 * N / 10.0;
      double d = fabs (count[i] - expected);
      double sigma = d / sqrt (expected);
      if (sigma > 5 && d > 1)
        {
          status = 1;
          gsl_test (status,
                    "gsl_ran_choose %d (%g observed vs %g expected)",
                    i, count[i] / N, 0.1);
        }
    }

  gsl_test (status, "gsl_ran_choose (3) on {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}");

}

void
testMoments (double (*f) (void), const char *name,
             double a, double b, double p)
{
  int i;
  double count = 0, expected, sigma;
  int status;

  for (i = 0; i < N; i++)
    {
      double r = f ();
      if (r < b && r > a)
        count++;
    }

  expected = p * N;
  sigma = (expected > 0) ? fabs (count - expected) / sqrt (expected) : fabs(count - expected);

  status = (sigma > 3);

  gsl_test (status, "%s [%g,%g] (%g observed vs %g expected)",
            name, a, b, count / N, p);
}

#define BINS 100

typedef double pdf_func(double);

/* Keep track of invalid values during integration */
static int pdf_errors = 0;
static double pdf_errval = 0.0;

double
wrapper_function (double x, void *params)
{
  pdf_func * pdf = (pdf_func *)params;
  double P = pdf(x);
  if (!gsl_finite(P)) {
    pdf_errors++;
    pdf_errval = P;
    P = 0; /* skip invalid value now, but return pdf_errval at the end */
  }
  return P;
}

double
integrate (pdf_func * pdf, double a, double b)
{
  double result, abserr;
  size_t n = 1000;
  gsl_function f;
  gsl_integration_workspace * w = gsl_integration_workspace_alloc (n);
  f.function = &wrapper_function;
  f.params = (void *)pdf;
  pdf_errors = 0;
  gsl_integration_qags (&f, a, b, 1e-16, 1e-4, n, w, &result, &abserr);
  gsl_integration_workspace_free (w);
  if (pdf_errors) return pdf_errval;
  return result;
}

void
testPDF (double (*f) (void), double (*pdf) (double), const char *name)
{
  double count[BINS], edge[BINS], p[BINS];
  double a = -5.0, b = +5.0;
  double dx = (b - a) / BINS;
  double bin;
  double total = 0, mean;
  int i, j, status = 0, status_i = 0, attempts = 0;
  long int n0 = 0, n = N;

  for (i = 0; i < BINS; i++)
    {
      /* Compute the integral of p(x) from x to x+dx */

      double x = a + i * dx;

      if (fabs (x) < 1e-10)     /* hit the origin exactly */
        x = 0.0;

      p[i]  = integrate (pdf, x, x+dx);
    }

  for (i = 0; i < BINS; i++)
    {
      count[i] = 0;
      edge[i] = 0;
    }

 trial:
  attempts++;

  for (i = n0; i < n; i++)
    {
      double r = f ();
      total += r;

      if (r < b && r > a)
        {
          double u = (r - a) / dx;
          double f = modf(u, &bin);
          j = (int)bin;

          if (f == 0)
            edge[j]++;
          else
            count[j]++;
        }
    }

  /* Sort out where the hits on the edges should go */

  for (i = 0; i < BINS; i++)
    {
      /* If the bin above is empty, its lower edge hits belong in the
         lower bin */

      if (i + 1 < BINS && count[i+1] == 0) {
        count[i] += edge[i+1];
        edge[i+1] = 0;
      }

      count[i] += edge[i];
      edge[i] = 0;
    }

  mean = (total / n);

  status = !gsl_finite(mean);
  if (status) {
    gsl_test (status, "%s, finite mean, observed %g", name, mean);
    return;
  }

  for (i = 0; i < BINS; i++)
    {
      double x = a + i * dx;
      double d = fabs (count[i] - n * p[i]);
      if (!gsl_finite(p[i]))
        {
          status_i = 1;
        }
      else if (p[i] != 0)
        {
          double s = d / sqrt (n * p[i]);
          status_i = (s > 5) && (d > 2);
        }
      else
        {
          status_i = (count[i] != 0);
        }

      /* Extend the sample if there is an outlier on the first attempt
         to avoid spurious failures when running large numbers of tests. */
      if (status_i && attempts < 50)
        {
          n0 = n;
          n = 2.0*n;
          goto trial;
        }

      status |= status_i;
      if (status_i)
        gsl_test (status_i, "%s [%g,%g) (%g/%d=%g observed vs %g expected)",
                  name, x, x + dx, count[i], n, count[i] / n, p[i]);
    }

  if (status == 0)
    gsl_test (status, "%s, sampling against pdf over range [%g,%g) ",
              name, a, b);
}

void
testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
                 const char *name)
{
  double count[BINS], p[BINS];
  unsigned int i;
  int status = 0, status_i = 0;

  for (i = 0; i < BINS; i++)
    count[i] = 0;

  for (i = 0; i < N; i++)
    {
      int r = (int) (f ());
      if (r >= 0 && r < BINS)
        count[r]++;
    }

  for (i = 0; i < BINS; i++)
    p[i] = pdf (i);

  for (i = 0; i < BINS; i++)
    {
      double d = fabs (count[i] - N * p[i]);
      if (p[i] != 0)
        {
          double s = d / sqrt (N * p[i]);
          status_i = (s > 5) && (d > 1);
        }
      else
        {
          status_i = (count[i] != 0);
        }
      status |= status_i;
      if (status_i)
        gsl_test (status_i, "%s i=%d (%g observed vs %g expected)",
                  name, i, count[i] / N, p[i]);
    }

  if (status == 0)
    gsl_test (status, "%s, sampling against pdf over range [%d,%d) ",
              name, 0, BINS);
}

double
test_beta (void)
{
  return gsl_ran_beta (r_global, 2.0, 3.0);
}

double
test_beta_pdf (double x)
{
  return gsl_ran_beta_pdf (x, 2.0, 3.0);
}

double
test_bernoulli (void)
{
  return gsl_ran_bernoulli (r_global, 0.3);
}

double
test_bernoulli_pdf (unsigned int n)
{
  return gsl_ran_bernoulli_pdf (n, 0.3);
}

double
test_binomial (void)
{
  return gsl_ran_binomial (r_global, 0.3, 5);
}

double
test_binomial_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 5);
}

double
test_binomial0 (void)
{
  return gsl_ran_binomial (r_global, 0, 8);
}

double
test_binomial0_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0, 8);
}

double
test_binomial1 (void)
{
  return gsl_ran_binomial (r_global, 1, 8);
}

double
test_binomial1_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 1, 8);
}

double
test_binomial_knuth (void)
{
  return gsl_ran_binomial_knuth (r_global, 0.3, 5);
}

double
test_binomial_knuth_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 5);
}

double
test_binomial_large (void)
{
  return gsl_ran_binomial (r_global, 0.3, 55);
}

double
test_binomial_large_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 55);
}

double
test_binomial_large_knuth (void)
{
  return gsl_ran_binomial_knuth (r_global, 0.3, 55);
}

double
test_binomial_large_knuth_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 55);
}

double
test_binomial_huge (void)
{
  return gsl_ran_binomial (r_global, 0.3, 5500);
}

double
test_binomial_huge_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 5500);
}

double
test_binomial_huge_knuth (void)
{
  return gsl_ran_binomial_knuth (r_global, 0.3, 5500);
}

double
test_binomial_huge_knuth_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 5500);
}

double
test_binomial_max (void)
{
  return gsl_ran_binomial (r_global, 1e-8, 1<<31);
}

double
test_binomial_max_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 1e-8, 1<<31);
}

double
test_cauchy (void)
{
  return gsl_ran_cauchy (r_global, 2.0);
}

double
test_cauchy_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 2.0);
}

double
test_chisq (void)
{
  return gsl_ran_chisq (r_global, 13.0);
}

double
test_chisq_pdf (double x)
{
  return gsl_ran_chisq_pdf (x, 13.0);
}

double
test_chisqnu2 (void)
{
  return gsl_ran_chisq (r_global, 2.0);
}

double
test_chisqnu2_pdf (double x)
{
  return gsl_ran_chisq_pdf (x, 2.0);
}

double
test_dir2d (void)
{
  double x = 0, y = 0, theta;
  gsl_ran_dir_2d (r_global, &x, &y);
  theta = atan2 (x, y);
  return theta;
}

double
test_dir2d_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI);
    }
  else
    {
      return 0;
    }
}

double
test_dir2d_trig_method (void)
{
  double x = 0, y = 0, theta;
  gsl_ran_dir_2d_trig_method (r_global, &x, &y);
  theta = atan2 (x, y);
  return theta;
}

double
test_dir2d_trig_method_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI);
    }
  else
    {
      return 0;
    }
}

double
test_dir3dxy (void)
{
  double x = 0, y = 0, z = 0, theta;
  gsl_ran_dir_3d (r_global, &x, &y, &z);
  theta = atan2 (x, y);
  return theta;
}

double
test_dir3dxy_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI);
    }
  else
    {
      return 0;
    }
}

double
test_dir3dyz (void)
{
  double x = 0, y = 0, z = 0, theta;
  gsl_ran_dir_3d (r_global, &x, &y, &z);
  theta = atan2 (y, z);
  return theta;
}

double
test_dir3dyz_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI);
    }
  else
    {
      return 0;
    }
}

double
test_dir3dzx (void)
{
  double x = 0, y = 0, z = 0, theta;
  gsl_ran_dir_3d (r_global, &x, &y, &z);
  theta = atan2 (z, x);
  return theta;
}

double
test_dir3dzx_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI);
    }
  else
    {
      return 0;
    }
}

double
test_dirichlet (void)
{
  /* This is a bit of a lame test, since when K=2, the Dirichlet distribution
     becomes a beta distribution */
  size_t K = 2;
  double alpha[2] = { 2.5, 5.0 };
  double theta[2] = { 0.0, 0.0 };

  gsl_ran_dirichlet (r_global, K, alpha, theta);

  return theta[0];
}

double
test_dirichlet_pdf (double x)
{
  size_t K = 2;
  double alpha[2] = { 2.5, 5.0 };
  double theta[2];

  if (x <= 0.0 || x >= 1.0)
    return 0.0;                 /* Out of range */

  theta[0] = x;
  theta[1] = 1.0 - x;

  return gsl_ran_dirichlet_pdf (K, alpha, theta);
}

double
test_dirichlet_small (void)
{
  size_t K = 2;
  double alpha[2] = { 2.5e-3, 5.0e-3};
  double theta[2] = { 0.0, 0.0 };

  gsl_ran_dirichlet (r_global, K, alpha, theta);

  return theta[0];
}

double
test_dirichlet_small_pdf (double x)
{
  size_t K = 2;
  double alpha[2] = { 2.5e-3, 5.0e-3 };
  double theta[2];

  if (x <= 0.0 || x >= 1.0)
    return 0.0;                 /* Out of range */

  theta[0] = x;
  theta[1] = 1.0 - x;

  return gsl_ran_dirichlet_pdf (K, alpha, theta);
}

/* Check that the observed means of the Dirichlet variables are
   within reasonable statistical errors of their correct values. */

#define DIRICHLET_K 10

void
test_dirichlet_moments (void)
{
  double alpha[DIRICHLET_K];
  double theta[DIRICHLET_K];
  double theta_sum[DIRICHLET_K];

  double alpha_sum = 0.0;
  double mean, obs_mean, sd, sigma;
  int status, k, n;

  for (k = 0; k < DIRICHLET_K; k++)
    {
      alpha[k] = gsl_ran_exponential (r_global, 0.1);
      alpha_sum += alpha[k];
      theta_sum[k] = 0.0;
    }

  for (n = 0; n < N; n++)
    {
      gsl_ran_dirichlet (r_global, DIRICHLET_K, alpha, theta);
      for (k = 0; k < DIRICHLET_K; k++)
        theta_sum[k] += theta[k];
    }

  for (k = 0; k < DIRICHLET_K; k++)
    {
      mean = alpha[k] / alpha_sum;
      sd =
        sqrt ((alpha[k] * (1. - alpha[k] / alpha_sum)) /
              (alpha_sum * (alpha_sum + 1.)));
      obs_mean = theta_sum[k] / N;
      sigma = sqrt ((double) N) * fabs (mean - obs_mean) / sd;

      status = (sigma > 3.0);

      gsl_test (status,
                "test gsl_ran_dirichlet: mean (%g observed vs %g expected)",
                obs_mean, mean);
    }
}

/* Check that the observed means of the multinomial variables are
   within reasonable statistical errors of their correct values. */

void
test_multinomial_moments (void)
{
  const unsigned int sum_n = 100;

  const double p[MULTI_DIM] ={ 0.2, 0.20, 0.17, 0.14, 0.12,
                               0.07, 0.05, 0.02, 0.02, 0.01 };

  unsigned int  x[MULTI_DIM];
  double x_sum[MULTI_DIM];

  double mean, obs_mean, sd, sigma;
  int status, k, n;

  for (k = 0; k < MULTI_DIM; k++)
    x_sum[k] =0.0;

  for (n = 0; n < N; n++)
    {
      gsl_ran_multinomial (r_global, MULTI_DIM, sum_n, p, x);
      for (k = 0; k < MULTI_DIM; k++)
        x_sum[k] += x[k];
    }

  for (k = 0; k < MULTI_DIM; k++)
    {
      mean = p[k] * sum_n;
      sd = p[k] * (1.-p[k]) * sum_n;

      obs_mean = x_sum[k] / N;
      sigma = sqrt ((double) N) * fabs (mean - obs_mean) / sd;

      status = (sigma > 3.0);

      gsl_test (status,
                "test gsl_ran_multinomial: mean (%g observed vs %g expected)",
                obs_mean, mean);
    }
}

double
test_discrete1 (void)
{
  static double P[3] = { 0.59, 0.4, 0.01 };
  if (g1 == NULL)
    {
      g1 = gsl_ran_discrete_preproc (3, P);
    }
  return gsl_ran_discrete (r_global, g1);
}

double
test_discrete1_pdf (unsigned int n)
{
  return gsl_ran_discrete_pdf ((size_t) n, g1);
}

double
test_discrete2 (void)
{
  static double P[10] = { 1, 9, 3, 4, 5, 8, 6, 7, 2, 0 };
  if (g2 == NULL)
    {
      g2 = gsl_ran_discrete_preproc (10, P);
    }
  return gsl_ran_discrete (r_global, g2);
}

double
test_discrete2_pdf (unsigned int n)
{
  return gsl_ran_discrete_pdf ((size_t) n, g2);
}
double
test_discrete3 (void)
{
  static double P[20];
  if (g3 == NULL)
    { int i;
      for (i=0; i<20; ++i) P[i]=1.0/20;
      g3 = gsl_ran_discrete_preproc (20, P);
    }
  return gsl_ran_discrete (r_global, g3);
}

double
test_discrete3_pdf (unsigned int n)
{
  return gsl_ran_discrete_pdf ((size_t) n, g3);
}

double
test_erlang (void)
{
  return gsl_ran_erlang (r_global, 3.0, 4.0);
}

double
test_erlang_pdf (double x)
{
  return gsl_ran_erlang_pdf (x, 3.0, 4.0);
}

double
test_exponential (void)
{
  return gsl_ran_exponential (r_global, 2.0);
}

double
test_exponential_pdf (double x)
{
  return gsl_ran_exponential_pdf (x, 2.0);
}

double
test_exppow0 (void)
{
  return gsl_ran_exppow (r_global, 3.7, 0.3);
}

double
test_exppow0_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 0.3);
}

double
test_exppow1 (void)
{
  return gsl_ran_exppow (r_global, 3.7, 1.0);
}

double
test_exppow1_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 1.0);
}

double
test_exppow1a (void)
{
  return gsl_ran_exppow (r_global, 3.7, 1.9);
}

double
test_exppow1a_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 1.9);
}

double
test_exppow2 (void)
{
  return gsl_ran_exppow (r_global, 3.7, 2.0);
}

double
test_exppow2_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 2.0);
}

double
test_exppow2a (void)
{
  return gsl_ran_exppow (r_global, 3.7, 3.5);
}

double
test_exppow2a_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 3.5);
}

double
test_exppow2b (void)
{
  return gsl_ran_exppow (r_global, 3.7, 7.5);
}

double
test_exppow2b_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 7.5);
}

double
test_fdist (void)
{
  return gsl_ran_fdist (r_global, 3.0, 4.0);
}

double
test_fdist_pdf (double x)
{
  return gsl_ran_fdist_pdf (x, 3.0, 4.0);
}

/* Test case for bug #28500: overflow in gsl_ran_fdist_pdf */

double
test_fdist_large (void)
{
  return gsl_ran_fdist (r_global, 8.0, 249.0);
}

double
test_fdist_large_pdf (double x)
{
  return gsl_ran_fdist_pdf (x, 8.0, 249.0);
}

double
test_flat (void)
{
  return gsl_ran_flat (r_global, 3.0, 4.0);
}

double
test_flat_pdf (double x)
{
  return gsl_ran_flat_pdf (x, 3.0, 4.0);
}

double
test_gamma (void)
{
  return gsl_ran_gamma (r_global, 2.5, 2.17);
}

double
test_gamma_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 2.5, 2.17);
}

double
test_gamma1 (void)
{
  return gsl_ran_gamma (r_global, 1.0, 2.17);
}

double
test_gamma1_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 1.0, 2.17);
}

double
test_gamma_int (void)
{
  return gsl_ran_gamma (r_global, 10.0, 2.17);
}

double
test_gamma_int_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 10.0, 2.17);
}

double
test_gamma_large (void)
{
  return gsl_ran_gamma (r_global, 20.0, 2.17);
}

double
test_gamma_large_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 20.0, 2.17);
}

double
test_gamma_small (void)
{
  return gsl_ran_gamma (r_global, 0.92, 2.17);
}

double
test_gamma_small_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 0.92, 2.17);
}

double
test_gamma_vlarge (void)
{
  /* Scale the distribution to get it into the range [-5,5] */
  double c = 2.71828181565;
  double b = 6.32899304917e-10;
  double d = 1e4;
  return (gsl_ran_gamma (r_global, 4294967296.0, b) - c) * d;
}

double
test_gamma_vlarge_pdf (double x)
{
  double c = 2.71828181565;
  double b = 6.32899304917e-10;
  double d = 1e4;
  return gsl_ran_gamma_pdf ((x / d) + c, 4294967296.0, b) / d;
}

double
test_gamma_mt (void)
{
  return gsl_ran_gamma_mt (r_global, 2.5, 2.17);
}

double
test_gamma_mt_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 2.5, 2.17);
}

double
test_gamma_mt1 (void)
{
  return gsl_ran_gamma_mt (r_global, 1.0, 2.17);
}

double
test_gamma_mt1_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 1.0, 2.17);
}

double
test_gamma_mt_int (void)
{
  return gsl_ran_gamma_mt (r_global, 10.0, 2.17);
}

double
test_gamma_mt_int_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 10.0, 2.17);
}

double
test_gamma_mt_large (void)
{
  return gsl_ran_gamma_mt (r_global, 20.0, 2.17);
}

double
test_gamma_mt_large_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 20.0, 2.17);
}

double
test_gamma_mt_small (void)
{
  return gsl_ran_gamma_mt (r_global, 0.92, 2.17);
}

double
test_gamma_mt_small_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 0.92, 2.17);
}

double
test_gamma_knuth_vlarge (void)
{
  /* Scale the distribution to get it into the range [-5,5] */
  double c = 2.71828181565;
  double b = 6.32899304917e-10;
  double d = 1e4;
  return (gsl_ran_gamma_knuth (r_global, 4294967296.0, b) - c) * d;
}

double
test_gamma_knuth_vlarge_pdf (double x)
{
  double c = 2.71828181565;
  double b = 6.32899304917e-10;
  double d = 1e4;
  return gsl_ran_gamma_pdf ((x / d) + c, 4294967296.0, b) / d;
}

double
test_gaussian (void)
{
  return gsl_ran_gaussian (r_global, 3.0);
}

double
test_gaussian_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, 3.0);
}

double
test_gaussian_ratio_method (void)
{
  return gsl_ran_gaussian_ratio_method (r_global, 3.0);
}

double
test_gaussian_ratio_method_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, 3.0);
}

double
test_gaussian_ziggurat (void)
{
  return gsl_ran_gaussian_ziggurat (r_global, 3.12);
}

double
test_gaussian_ziggurat_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, 3.12);
}

double
test_gaussian_tail (void)
{
  return gsl_ran_gaussian_tail (r_global, 1.7, 0.25);
}

double
test_gaussian_tail_pdf (double x)
{
  return gsl_ran_gaussian_tail_pdf (x, 1.7, 0.25);
}

double
test_gaussian_tail1 (void)
{
  return gsl_ran_gaussian_tail (r_global, -1.7, 5.0);
}

double
test_gaussian_tail1_pdf (double x)
{
  return gsl_ran_gaussian_tail_pdf (x, -1.7, 5.0);
}

double
test_gaussian_tail2 (void)
{
  return gsl_ran_gaussian_tail (r_global, 0.1, 2.0);
}

double
test_gaussian_tail2_pdf (double x)
{
  return gsl_ran_gaussian_tail_pdf (x, 0.1, 2.0);
}

double
test_ugaussian (void)
{
  return gsl_ran_ugaussian (r_global);
}

double
test_ugaussian_pdf (double x)
{
  return gsl_ran_ugaussian_pdf (x);
}

double
test_ugaussian_ratio_method (void)
{
  return gsl_ran_ugaussian_ratio_method (r_global);
}

double
test_ugaussian_ratio_method_pdf (double x)
{
  return gsl_ran_ugaussian_pdf (x);
}

double
test_ugaussian_tail (void)
{
  return gsl_ran_ugaussian_tail (r_global, 3.0);
}

double
test_ugaussian_tail_pdf (double x)
{
  return gsl_ran_ugaussian_tail_pdf (x, 3.0);
}

double
test_bivariate_gaussian1 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
  return x;
}

double
test_bivariate_gaussian1_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, 3.0);
}

double
test_bivariate_gaussian2 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
  return y;
}

double
test_bivariate_gaussian2_pdf (double y)
{
  int i, n = 10;
  double sum = 0;
  double a = -10, b = 10, dx = (b - a) / n;
  for (i = 0; i < n; i++)
    {
      double x = a + i * dx;
      sum += gsl_ran_bivariate_gaussian_pdf (x, y, 3.0, 2.0, 0.3) * dx;
    }
  return sum;
}

double
test_bivariate_gaussian3 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
  return x + y;
}

double
test_bivariate_gaussian3_pdf (double x)
{
  double sx = 3.0, sy = 2.0, r = 0.3;
  double su = (sx + r * sy);
  double sv = sy * sqrt (1 - r * r);
  double sigma = sqrt (su * su + sv * sv);

  return gsl_ran_gaussian_pdf (x, sigma);
}

double
test_bivariate_gaussian4 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, -0.5, &x, &y);
  return x + y;
}

double
test_bivariate_gaussian4_pdf (double x)
{
  double sx = 3.0, sy = 2.0, r = -0.5;
  double su = (sx + r * sy);
  double sv = sy * sqrt (1 - r * r);
  double sigma = sqrt (su * su + sv * sv);

  return gsl_ran_gaussian_pdf (x, sigma);
}

/* Examples from R (GPL): http://www.r-project.org/
 * library(mvtnorm); packageVersion("mvtnorm") # 1.0.5
 * mu <- c(1, 2)
 * Sigma <- matrix(c(4,2, 2,3), ncol=2)
 * x <- c(0, 0)
 * sprintf("%.15f", dmvnorm(x=x, mean=mu, sigma=Sigma, log=TRUE)) # -3.565097837249263
 */
void
test_multivariate_gaussian_log_pdf (void)
{
  size_t d = 2;
  const double exp_res = -3.565097837249263;
  double obs_res;
  gsl_vector * mu = gsl_vector_calloc(d);
  gsl_matrix * Sigma = gsl_matrix_calloc(d, d);
  gsl_matrix * L = gsl_matrix_calloc(d, d);
  gsl_vector * x = gsl_vector_calloc(d);
  gsl_vector * work = gsl_vector_calloc(d);

  gsl_vector_set(mu, 0, 1);
  gsl_vector_set(mu, 1, 2);

  gsl_matrix_set(Sigma, 0, 0, 4);
  gsl_matrix_set(Sigma, 1, 1, 3);
  gsl_matrix_set(Sigma, 0, 1, 2);
  gsl_matrix_set(Sigma, 1, 0, 2);

  gsl_matrix_memcpy(L, Sigma);
  gsl_linalg_cholesky_decomp1(L);

  gsl_ran_multivariate_gaussian_log_pdf(x, mu, L, &obs_res, work);
  gsl_test_rel(obs_res, exp_res, 1.0e-10, "gsl_ran_multivariate_gaussian_log_pdf");

  gsl_vector_free(mu);
  gsl_matrix_free(Sigma);
  gsl_matrix_free(L);
  gsl_vector_free(x);
  gsl_vector_free(work);
}

/* Examples from R (GPL): http://www.r-project.org/
 * library(mvtnorm); packageVersion("mvtnorm") # 1.0.5
 * mu <- c(1, 2)
 * Sigma <- matrix(c(4,2, 2,3), ncol=2)
 * x <- c(0, 0)
 * sprintf("%.15f", dmvnorm(x=x, mean=mu, sigma=Sigma, log=FALSE)) # 0.028294217120391
 */
void
test_multivariate_gaussian_pdf (void)
{
  size_t d = 2;
  const double exp_res = 0.028294217120391;
  double obs_res = 0;
  gsl_vector * mu = gsl_vector_calloc(d);
  gsl_matrix * Sigma = gsl_matrix_calloc(d, d);
  gsl_matrix * L = gsl_matrix_calloc(d, d);
  gsl_vector * x = gsl_vector_calloc(d);
  gsl_vector * work = gsl_vector_calloc(d);

  gsl_vector_set(mu, 0, 1);
  gsl_vector_set(mu, 1, 2);

  gsl_matrix_set(Sigma, 0, 0, 4);
  gsl_matrix_set(Sigma, 1, 1, 3);
  gsl_matrix_set(Sigma, 0, 1, 2);
  gsl_matrix_set(Sigma, 1, 0, 2);

  gsl_matrix_memcpy(L, Sigma);
  gsl_linalg_cholesky_decomp1(L);

  gsl_ran_multivariate_gaussian_pdf(x, mu, L, &obs_res, work);
  gsl_test_rel(obs_res, exp_res, 1.0e-10, "gsl_ran_multivariate_gaussian_pdf");

  gsl_vector_free(mu);
  gsl_matrix_free(Sigma);
  gsl_matrix_free(L);
  gsl_vector_free(x);
  gsl_vector_free(work);
}

/* Draw N random vectors according to a given MVN(mu,Sigma). Then, check that
 * one can't reject the null hypothesis that the sample mean is equal to
 * the true mean, using Hotelling's test statistic at 95% confidence level.
 * Details in "Applied Multivariate Statistical Analysis" by  Johnson & Wichern
 * (2001), section 5, page 212.
 */
void
test_multivariate_gaussian (void)
{
  size_t d = 2, i = 0;
  int status = 0;
  double T2 = 0, threshold = 0, alpha = 0.05, pvalue = 0;
  gsl_vector * mu = gsl_vector_calloc(d);
  gsl_matrix * Sigma = gsl_matrix_calloc(d, d);
  gsl_matrix * L = gsl_matrix_calloc(d, d);
  gsl_vector * sample = gsl_vector_calloc(d);
  gsl_matrix * samples = gsl_matrix_calloc(N, d);
  gsl_vector * mu_hat = gsl_vector_calloc(d);
  gsl_matrix * Sigma_hat = gsl_matrix_calloc(d, d);
  gsl_vector * mu_hat_ctr = gsl_vector_calloc(d);
  gsl_matrix * Sigma_hat_inv = gsl_matrix_calloc(d, d);
  gsl_vector * tmp = gsl_vector_calloc(d);

  /* set the true values of parameters mu and Sigma */
  gsl_vector_set(mu, 0, 1);
  gsl_vector_set(mu, 1, 2);
  gsl_matrix_set(Sigma, 0, 0, 4);
  gsl_matrix_set(Sigma, 1, 1, 3);
  gsl_matrix_set(Sigma, 0, 1, 2);
  gsl_matrix_set(Sigma, 1, 0, 2);

  /* draw N random vectors */
  gsl_matrix_memcpy(L, Sigma);
  gsl_linalg_cholesky_decomp1(L);
  for (i = 0; i < N; ++i) {
    gsl_ran_multivariate_gaussian(r_global, mu, L, sample);
    gsl_matrix_set_row(samples, i, sample);
  }

  /* compute the maximum-likelihood estimates */
  gsl_ran_multivariate_gaussian_mean (samples, mu_hat);
  gsl_ran_multivariate_gaussian_vcov (samples, Sigma_hat);

  /* compute Hotelling's test statistic:
     T^2 = n (hat{mu} - mu)' hat{Sigma}^-1 (hat{mu} - mu) */
  gsl_vector_memcpy(mu_hat_ctr, mu_hat);
  gsl_vector_sub(mu_hat_ctr, mu);
  gsl_matrix_memcpy(Sigma_hat_inv, Sigma_hat);
  gsl_linalg_cholesky_decomp1(Sigma_hat_inv);
  gsl_linalg_cholesky_invert(Sigma_hat_inv);
  gsl_blas_dgemv(CblasNoTrans, 1, Sigma_hat_inv, mu_hat_ctr, 0, tmp);
  gsl_blas_ddot(mu_hat_ctr, tmp, &T2);
  T2 *= N;

  /* test if the null hypothesis (hat{mu} = mu) can be rejected
   at the alpha level*/
  threshold = (N-1) * d / (double)(N-d) * gsl_cdf_fdist_Pinv(1-alpha, d, N-d);
  status = (T2 > threshold);
  gsl_test(status,
           "test gsl_ran_multivariate_gaussian: T2 %f < %f",
           T2, threshold);
  pvalue = gsl_cdf_fdist_Q(T2, d, N-d);
  status = (pvalue < alpha);
  gsl_test(status,
           "test gsl_ran_multivariate_gaussian: p value %f > %f",
           pvalue, alpha);

  gsl_vector_free(mu);
  gsl_matrix_free(Sigma);
  gsl_matrix_free(L);
  gsl_vector_free(sample);
  gsl_matrix_free(samples);
  gsl_vector_free(mu_hat);
  gsl_matrix_free(Sigma_hat);
  gsl_vector_free(mu_hat_ctr);
  gsl_matrix_free(Sigma_hat_inv);
  gsl_vector_free(tmp);
}

double
test_geometric (void)
{
  return gsl_ran_geometric (r_global, 0.5);
}

double
test_geometric_pdf (unsigned int n)
{
  return gsl_ran_geometric_pdf (n, 0.5);
}

double
test_geometric1 (void)
{
  return gsl_ran_geometric (r_global, 1.0);
}

double
test_geometric1_pdf (unsigned int n)
{
  return gsl_ran_geometric_pdf (n, 1.0);
}

double
test_hypergeometric1 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 4);
}

double
test_hypergeometric1_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 4);
}

double
test_hypergeometric2 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 11);
}

double
test_hypergeometric2_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 11);
}

double
test_hypergeometric3 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 1);
}

double
test_hypergeometric3_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 1);
}

double
test_hypergeometric4 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 20);
}

double
test_hypergeometric4_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 20);
}

double
test_hypergeometric5 (void)
{
  return gsl_ran_hypergeometric (r_global, 2, 7, 5);
}

double
test_hypergeometric5_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 2, 7, 5);
}

double
test_hypergeometric6 (void)
{
  return gsl_ran_hypergeometric (r_global, 2, 10, 3);
}

double
test_hypergeometric6_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 2, 10, 3);
}

double
test_gumbel1 (void)
{
  return gsl_ran_gumbel1 (r_global, 3.12, 4.56);
}

double
test_gumbel1_pdf (double x)
{
  return gsl_ran_gumbel1_pdf (x, 3.12, 4.56);
}

double
test_gumbel2 (void)
{
  return gsl_ran_gumbel2 (r_global, 3.12, 4.56);
}

double
test_gumbel2_pdf (double x)
{
  return gsl_ran_gumbel2_pdf (x, 3.12, 4.56);
}

double
test_landau (void)
{
  return gsl_ran_landau (r_global);
}

double
test_landau_pdf (double x)
{
  return gsl_ran_landau_pdf (x);
}

double
test_levy1 (void)
{
  return gsl_ran_levy (r_global, 5.0, 1.0);
}

double
test_levy1_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy2 (void)
{
  return gsl_ran_levy (r_global, 5.0, 2.0);
}

double
test_levy2_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
}

double
test_levy1a (void)
{
  return gsl_ran_levy (r_global, 5.0, 1.01);
}

double
test_levy1a_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy2a (void)
{
  return gsl_ran_levy (r_global, 5.0, 1.99);
}

double
test_levy2a_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
}

double
test_levy_skew1 (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.0, 0.0);
}

double
test_levy_skew1_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy_skew2 (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 2.0, 0.0);
}

double
test_levy_skew2_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
}

double
test_levy_skew1a (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.01, 0.0);
}

double
test_levy_skew1a_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy_skew2a (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.99, 0.0);
}

double
test_levy_skew2a_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
}

double
test_levy_skew1b (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.01, 0.001);
}

double
test_levy_skew1b_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy_skew2b (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.99, 0.001);
}

double
test_levy_skew2b_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
}

double
test_logistic (void)
{
  return gsl_ran_logistic (r_global, 3.1);
}

double
test_logistic_pdf (double x)
{
  return gsl_ran_logistic_pdf (x, 3.1);
}

double
test_logarithmic (void)
{
  return gsl_ran_logarithmic (r_global, 0.4);
}

double
test_logarithmic_pdf (unsigned int n)
{
  return gsl_ran_logarithmic_pdf (n, 0.4);
}

double
test_lognormal (void)
{
  return gsl_ran_lognormal (r_global, 2.7, 1.3);
}

double
test_lognormal_pdf (double x)
{
  return gsl_ran_lognormal_pdf (x, 2.7, 1.3);
}

double
test_multinomial (void)
{
  const size_t K = 3;
  const unsigned int sum_n = BINS;
  unsigned int n[3];
  /* Test use of weights instead of probabilities. */
  const double p[] = { 2., 7., 1.};

  gsl_ran_multinomial ( r_global, K, sum_n, p, n);

  return n[0];
}

double
test_multinomial_pdf (unsigned int n_0)
{
  /* The margional distribution of just 1 variate  is binomial. */
  size_t K = 2;
  /* Test use of weights instead of probabilities */
  double p[] = { 0.4, 1.6};
  const unsigned int sum_n = BINS;
  unsigned int n[2];

  n[0] = n_0;
  n[1] =sum_n - n_0;

  return gsl_ran_multinomial_pdf (K, p, n);
}

double
test_multinomial_large (void)
{
  const unsigned int sum_n = BINS;
  unsigned int n[MULTI_DIM];
  const double p[MULTI_DIM] = { 0.2, 0.20, 0.17, 0.14, 0.12,
                                0.07, 0.05, 0.04, 0.01, 0.00  };

  gsl_ran_multinomial ( r_global, MULTI_DIM, sum_n, p, n);

  return n[0];
}

double
test_multinomial_large_pdf (unsigned int n_0)
{
  return test_multinomial_pdf(n_0);
}

double
test_negative_binomial (void)
{
  return gsl_ran_negative_binomial (r_global, 0.3, 20.0);
}

double
test_negative_binomial_pdf (unsigned int n)
{
  return gsl_ran_negative_binomial_pdf (n, 0.3, 20.0);
}

double
test_pascal (void)
{
  return gsl_ran_pascal (r_global, 0.8, 3);
}

double
test_pascal_pdf (unsigned int n)
{
  return gsl_ran_pascal_pdf (n, 0.8, 3);
}

double
test_pareto (void)
{
  return gsl_ran_pareto (r_global, 1.9, 2.75);
}

double
test_pareto_pdf (double x)
{
  return gsl_ran_pareto_pdf (x, 1.9, 2.75);
}

double
test_rayleigh (void)
{
  return gsl_ran_rayleigh (r_global, 1.9);
}

double
test_rayleigh_pdf (double x)
{
  return gsl_ran_rayleigh_pdf (x, 1.9);
}

double
test_rayleigh_tail (void)
{
  return gsl_ran_rayleigh_tail (r_global, 2.7, 1.9);
}

double
test_rayleigh_tail_pdf (double x)
{
  return gsl_ran_rayleigh_tail_pdf (x, 2.7, 1.9);
}

double
test_poisson (void)
{
  return gsl_ran_poisson (r_global, 5.0);
}

double
test_poisson_pdf (unsigned int n)
{
  return gsl_ran_poisson_pdf (n, 5.0);
}

double
test_poisson_large (void)
{
  return gsl_ran_poisson (r_global, 30.0);
}

double
test_poisson_large_pdf (unsigned int n)
{
  return gsl_ran_poisson_pdf (n, 30.0);
}

double
test_tdist1 (void)
{
  return gsl_ran_tdist (r_global, 1.75);
}

double
test_tdist1_pdf (double x)
{
  return gsl_ran_tdist_pdf (x, 1.75);
}

double
test_tdist2 (void)
{
  return gsl_ran_tdist (r_global, 12.75);
}

double
test_tdist2_pdf (double x)
{
  return gsl_ran_tdist_pdf (x, 12.75);
}

double
test_laplace (void)
{
  return gsl_ran_laplace (r_global, 2.75);
}

double
test_laplace_pdf (double x)
{
  return gsl_ran_laplace_pdf (x, 2.75);
}

double
test_weibull (void)
{
  return gsl_ran_weibull (r_global, 3.14, 2.75);
}

double
test_weibull_pdf (double x)
{
  return gsl_ran_weibull_pdf (x, 3.14, 2.75);
}

double
test_weibull1 (void)
{
  return gsl_ran_weibull (r_global, 2.97, 1.0);
}

double
test_weibull1_pdf (double x)
{
  return gsl_ran_weibull_pdf (x, 2.97, 1.0);
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef __MATH_FUNCTINS_H__
#define __MATH_FUNCTINS_H__

#include <stdint.h>
#include <cmath>  // for std::fabs and std::signbit
#include <cblas.h>
#include <cudnn.h>
#include <cublas_v2.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <curand.h>
#include <driver_types.h>  // cuda driver types
#include <algorithm>

#include <glog/logging.h>
#define PERMUTELAYER_ORDERNUM 4
#define BLOCK 512
//
// CUDA macros
//

// CUDA: various checks for different function calls.
#define CUDA_CHECK(condition) \
  /* Code block avoids redefinition of cudaError_t error */ \
  do { \
    cudaError_t error = condition; \
    CHECK_EQ(error, cudaSuccess) << " " << cudaGetErrorString(error); \
  } while (0)

#define CUBLAS_CHECK(condition) \
  do { \
    cublasStatus_t status = condition; \
    CHECK_EQ(status, CUBLAS_STATUS_SUCCESS) << " " \
      << cublasGetErrorString(status); \
  } while (0)

#define CURAND_CHECK(condition) \
  do { \
    curandStatus_t status = condition; \
    CHECK_EQ(status, CURAND_STATUS_SUCCESS) << " " \
      << curandGetErrorString(status); \
  } while (0)

// CUDA: grid stride looping
#define CUDA_KERNEL_LOOP(i, n) \
  for (int i = blockIdx.x * blockDim.x + threadIdx.x; \
       i < (n); \
       i += blockDim.x * gridDim.x)

// CUDA: check for error after kernel execution and exit loudly if there is one.
#define CUDA_POST_KERNEL_CHECK CUDA_CHECK(cudaPeekAtLastError())

// CUDA: library error reporting.
const char* cublasGetErrorString(cublasStatus_t error);
const char* curandGetErrorString(curandStatus_t error);

// CUDA: use 512 threads per block
const int TENSORRT_CUDA_NUM_THREADS = 256;

// CUDA: number of blocks for threads.
inline int TENSORRT_GET_BLOCKS(const int N) {
  return (N + TENSORRT_CUDA_NUM_THREADS - 1) / TENSORRT_CUDA_NUM_THREADS;
}

/*
 * function: X[i] = alpha,initialize X with constant alpha
 *
 */
template <typename Dtype>
void tensorrt_gpu_set(const int N, const Dtype alpha, Dtype *X);

/*
 * function: y[index] = pow(a[index], alpha)
 *@params n: the dims of matrix a
 *@params a: matrix
 *@params y: vector
 */
template <typename Dtype>
void tensorrt_gpu_powx(const int n, const Dtype* a, const Dtype alpha, Dtype* y);

/*
 *function:y = alpha*A*x + beta*y;
 *@params handle: handle
 *@params TransA: transpose flag
 *@params M: the rows of A
 *@params N: the cols of A
 *@params alpha: the coefficient of A*x
 *@params A: matrix [M x N]
 *@params x: vector x
 *@params beta: the coefficient of y
 *@params y: vector y
 */
template <typename Dtype>
void tensorrt_gpu_gemv(cublasHandle_t handle,const CBLAS_TRANSPOSE TransA, const int M, const int N,
    const Dtype alpha, const Dtype* A, const Dtype* x, const Dtype beta,
    Dtype* y);

template <typename Dtype>
void tensorrt_gpu_divbsx(const int nthreads, const Dtype* A,
    const Dtype* v, const int rows, const int cols, const CBLAS_TRANSPOSE trans,
    Dtype* B);

template <typename Dtype>
void tensorrt_gpu_mulbsx(const int nthreads, const Dtype* A,
    const Dtype* v, const int rows, const int cols, const CBLAS_TRANSPOSE trans,
    Dtype* B);
cudaError_t tensorrt_gpu_permute(const int nthreads,float* const  bottom_data,const bool forward,
	const int* permute_order,const int* old_steps,const int* new_steps,const int num_axes,float* const top_data,cudaStream_t stream);

cudaError_t SoftmaxLayer(const float *bottom_data, int count, int channels, int outer_num_, int inner_num_, float *scale_data, float *top_data, cudaStream_t stream);

cudaError_t ConcatLayer(int nthreads, const float *bottom_data, bool kForward, int num_concats_, int concat_input_size_, int top_concat_axis, int bottom_concat_axis, int offset_concat_axis, float *top_data, cudaStream_t stream);

//cudaError_t cudaSoftmax(int n, int channels,  float* x, float*y, cudaStream_t stream);

//virtual void Forward_gpu(const vector<Blob<Dtype>*>& bottom,const vector<Blob<Dtype>*>& top);
cudaError_t cudaSoftmax_caffe(int count,int channels,float* x,float* y, cudaStream_t stream);

cudaError_t cudaDetectionOutput_caffe( int bottom0_count,
                                       int bottom1_count,
                                       float* loc_data,
                                       float* bottom1,
                                       float* prior_data,
                                       float* bottom3,
                                       float* bottom4,
                                       float* y,
                                       cudaStream_t stream);

#endif

ext: h


text:
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram2d.h>

int
main (void)
{
  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_histogram2d * h = gsl_histogram2d_alloc (10, 10);

  gsl_histogram2d_set_ranges_uniform (h,
                                      0.0, 1.0,
                                      0.0, 1.0);

  gsl_histogram2d_accumulate (h, 0.3, 0.3, 1);
  gsl_histogram2d_accumulate (h, 0.8, 0.1, 5);
  gsl_histogram2d_accumulate (h, 0.7, 0.9, 0.5);

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  {
    int i;
    gsl_histogram2d_pdf * p
      = gsl_histogram2d_pdf_alloc (h->nx, h->ny);

    gsl_histogram2d_pdf_init (p, h);

    for (i = 0; i < 1000; i++) {
      double x, y;
      double u = gsl_rng_uniform (r);
      double v = gsl_rng_uniform (r);

      gsl_histogram2d_pdf_sample (p, u, v, &x, &y);

      printf ("%g %g\n", x, y);
    }

    gsl_histogram2d_pdf_free (p);
  }

  gsl_histogram2d_free (h);
  gsl_rng_free (r);

  return 0;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified


text:
// Copyright (C) 2021 Intel Corporation
// SPDX-License-Identifier: Apache-2.0

#ifndef _HEBench_ClearText_BE_DataPack_H_7e5fa8c2415240ea93eff148ed73539b
#define _HEBench_ClearText_BE_DataPack_H_7e5fa8c2415240ea93eff148ed73539b

#include <memory>
#include <vector>

#include <gsl/gsl>

#include "clear_error.h"

#include "hebench/api_bridge/types.h"

#include "hebench/api_bridge/cpp/hebench.hpp"

template <class T>
/**
 * @brief Encapsulates a HEBench DataPack.
 */
class ClearDataPack : public hebench::cpp::ITaggedObject
{
private:
    HEBERROR_DECLARE_CLASS_NAME(ClearDataPack)

public:
    explicit ClearDataPack(const hebench::APIBridge::DataPack &data_pack);
    /**
     * @brief Creates a new data pack with specified data sizes.
     * @param sample_count[in] Number of data samples to be contained in this pack.
     * @param p_sample_sizes[in] List of number of items per data sample.
     * @param position[in] Parameter position corresponding to hebench::APIBridge::DataPack::param_position.
     * @details This constructor initializes vector<vector<T>> data structure of the shape
     * @code
     * data = vector<vector<T>>(sample_count)
     * @endcode
     * and
     * @code
     * data[i] = vector<T>(p_sample_sizes[i])
     * @endcode
     */
    explicit ClearDataPack(std::uint64_t sample_count,
                           const std::uint64_t *p_sample_sizes,
                           std::uint64_t position);
    ~ClearDataPack() override {}

    /**
     * @brief Fills pre-allocated HEBench data pack with the data contained in this
     * data pack.
     * @param data_pack[out] HEBench data pack where to store the data.
     * @param b_pad[in] If true, any extra space in \p data_pack will be padded with
     * zeroes. Otherwise, extra space will not be overwritten.
     * @details This method is the inverse of the constructor. Thus, if the constructor
     * performs some type of conversion, it must be reversible, and this method
     * reverses the conversion.
     *
     * Default implementation is a simple copy of the data.
     */
    virtual void fillHEBenchDataPack(hebench::APIBridge::DataPack &data_pack, bool b_pad) const;

    std::uint64_t getSamplesCount() const { return m_data.size(); }
    gsl::span<T> getSample(std::uint64_t sample_index);
    gsl::span<const T> getSample(std::uint64_t sample_index) const;
    /**
     * @brief Copies the data contained in \p sample_data into this data pack's
     * sample specified by \p sample index.
     * @param sample_index[in] Index of data sample in this data pack to overwrite.
     * @param sample_data[in] Data to copy into this data pack sample.
     * @details The size of this data pack sample will be resized to match the
     * size of specified \p sample-data.
     */
    void setSample(std::uint64_t sample_index, const gsl::span<T> &sample_data);
    /**
     * @brief Copies the data contained in \p sample_data into this data pack's
     * sample specified by \p sample index.
     * @param sample_index[in] Index of data sample in this data pack to overwrite.
     * @param sample_data[in] Data to copy into this data pack sample.
     * @param pad_value[in] Value to be used for padding if needed.
     * @details The size of this data pack sample will remain constant. If not
     * enough space, this method will copy as much data as it can from \p sample_data
     * into the specified sample. If this data pack's sample's size is more than
     * the supplied data in \p sample_data, the remaining space will be padded
     * with the specified padding value \p pad_value.
     */
    void setSample(std::uint64_t sample_index, const gsl::span<T> &sample_data, const T &pad_value);

    std::uint64_t position() const { return m_position; }

protected:
    std::vector<std::vector<T>> m_data;

    ClearDataPack(std::uint64_t position = 0) :
        m_position(position) {}
    void allocateBuffers(std::uint64_t sample_count, const std::uint64_t *p_sample_sizes);
    void allocateBuffers(const hebench::APIBridge::DataPack &data_pack);
    /**
     * @brief Performs a simple copy.
     * @param data_pack
     */
    virtual void fillFromHEBenchDataPack(const hebench::APIBridge::DataPack &data_pack);

private:
    std::uint64_t m_position;
};

#include "inl/data_container.inl"

#endif // defined _HEBench_ClearText_BE_DataPack_H_7e5fa8c2415240ea93eff148ed73539b

max_stars_repo_name: jlakness-intel/backend-cpu-helib

max_issues_repo_name: jlakness-intel/backend-cpu-helib

max_forks_repo_name: jlakness-intel/backend-cpu-helib


text:
#pragma once

#include "poll_event.h"
#include "poll_response.h"
#include "poll_target.h"

#include <cstddef>
#include <optional>
#include <string>
#include <vector>

#include <gsl/span>

namespace wrappers::zmq {

class context;

class socket final {

public:
    enum class type {
        pair,
        pub,
        sub,
        req,
        rep,
        dealer,
        router,
        pull,
        push,
        xpub,
        xsub,
        stream,
    };

    explicit socket(context &ctx, type socket_type) noexcept;
    socket(const socket &other) = delete;
    socket &operator=(const socket &other) = delete;
    socket(socket &&other) noexcept;
    socket &operator=(socket &&other) noexcept;
    ~socket() noexcept;

    [[nodiscard]] bool bind(const std::string &endpoint) noexcept;

    [[nodiscard]] bool connect(const std::string &endpoint) noexcept;

    [[nodiscard]] bool blocking_send() noexcept;
    [[nodiscard]] bool blocking_send(gsl::span<std::byte> message) noexcept;

    [[nodiscard]] std::optional<std::vector<std::byte>>
    blocking_receive() noexcept;

    [[nodiscard]] bool
    async_receive(void *data,
                  void (*callback)(void *, gsl::span<std::byte>)) noexcept;

    friend std::optional<std::vector<poll_response>>
    blocking_poll(gsl::span<poll_target> targets) noexcept;

private:
    void *m_socket = nullptr;
};

} // namespace wrappers::zmq

max_stars_repo_name: Longhanks/linkollector-win

max_issues_repo_name: Longhanks/linkollector-win

max_forks_repo_name: Longhanks/linkollector-win


text:
/*
 * MathUtils.h
 *
 * Author:
 *       Oleg Kalashev
 *
 * Copyright (c) 2020 Institute for Nuclear Research, RAS
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#ifndef MATHUTILS_H_
#define MATHUTILS_H_

#include <gsl/gsl_roots.h>
#include <gsl/gsl_integration.h>
#include <iostream>
#include "Utils.h"
#include <limits>

namespace Utils {

template<typename X = double >
class FunctionX
{
public:
	virtual X f(X _x) const = 0;
	virtual X Xmin() const {return -std::numeric_limits<X>::max();}
	virtual X Xmax() const {return std::numeric_limits<X>::max();}
	inline X operator()(X _x) const {return f(_x);};
	virtual ~FunctionX(){};

	virtual FunctionX<X>* Clone() const { NOT_IMPLEMENTED; return 0; }

	inline operator gsl_function() const
	{
		gsl_function result = {GslProxyFunc, (void*)this};
		return result;
	}
	void Print(std::ostream& aOut, int nIntervals, bool aLogScale, double aXmin=-DBL_MAX, double aXmax=DBL_MAX) const
	{
		if(aXmin<Xmin())
			aXmin = Xmin();
		if(aXmax>Xmax())
			aXmax = Xmax();
		ASSERT(aXmin<=aXmax);
		ASSERT(aXmax<DBL_MAX);
		ASSERT(aXmin>-DBL_MAX);
		ASSERT(nIntervals>=1);
		ASSERT((!aLogScale) || aXmin>0);
		double step = aLogScale?pow(aXmax/aXmin,1./nIntervals) : (aXmax-aXmin)/nIntervals;
		double x = aXmin;
		for(int i=0; i<=nIntervals; i++, (x = aLogScale ? x*step : x+step) )
			aOut << x << "\t" << f(x) << "\n";
	}
private:
	static double GslProxyFunc (double x, void * params)
	{
		const FunctionX<X>* f = (const FunctionX<X>*)params;
		return f->f(x);
	}
};

template<typename X = double >
class Function2X
{
public:
	virtual X f(X x, X y) const = 0;
	virtual X MinArg(int aArgNo) const {return std::numeric_limits<X>::min();}
	virtual X MaxArg(int aArgNo) const  {return std::numeric_limits<X>::max();}

	virtual inline X operator()(X _x, X _y) const {return f(_x,_y);};
	virtual ~Function2X(){};
};

template<typename X = double >
class IFunctionCallHandlerX
{
public:
	virtual void OnCall(const FunctionX<X>& aFunc, X aX, X aY) const = 0;
};

template<typename X = double >
class DebugFunctionX : public FunctionX<X>
{
public:
	DebugFunctionX(const FunctionX<X>& aOrigFunc, const IFunctionCallHandlerX<X>& aDebugger):fOrigFunc(aOrigFunc),fDebugger(aDebugger){};
	virtual X f(X _x) const
	{
		X y = fOrigFunc.f(_x);
		fDebugger.OnCall(fOrigFunc, _x, y);
		return y;
	}
	virtual ~DebugFunctionX(){};
	virtual FunctionX<X>* Clone() const { return new DebugFunctionX<X>(fOrigFunc,fDebugger); }
private:
	const FunctionX<X>&				fOrigFunc;
	const IFunctionCallHandlerX<X>&	fDebugger;
};

template<typename X = double >
class FunctionCallLoggerX : public IFunctionCallHandlerX<X>
{
public:
	FunctionCallLoggerX(std::ostream& aOutput) : fOutput(aOutput){}
	virtual void OnCall(const FunctionX<X>& aFunc, X aX, X aY) const
	{
		((std::ostream&)fOutput) << aX << "\t" << aY << "\n";
	}
	virtual ~FunctionCallLoggerX()
	{
		fOutput << std::endl;
	}
private:
	std::ostream& fOutput;
};

typedef FunctionX<double> Function;
typedef Function2X<double> Function2;

template<typename X = double >
class ParamlessFunctionX : public FunctionX<X>
{
public:
	ParamlessFunctionX(X (*aFunction)(X)):fFunction(aFunction){}
	X f(X _x) const { return fFunction(_x); }
private:
	X (*fFunction)(X);
};

typedef ParamlessFunctionX<double> ParamlessFunction;

	template<typename X = double >
	class ISampler {
	public:
		virtual X sample(const FunctionX<X> &f, X aTmin, X aTmax, X aRand,
				 X &aTotRate, X aRelErr, X aAbsErr = 1e300) = 0;
		void UnitTest();
	};

/// Mathematical functions
/// Static methods are thread-safe
/// Non-static methods are not guaranteed to be thread-safe
	///TODO: make implementation using boost odeint and dense output and get rid of nr library
class MathUtils {
public:
	MathUtils();
	~MathUtils();
    static inline double RelDifference(double aVal1, double aVal2){
        double meanNorm = 0.5*fabs(aVal1)+fabs(aVal2);
        return meanNorm==0. ? 0.:fabs(aVal1-aVal2)/meanNorm;
    }

	static double SolveEquation(
			double (*aEquation) (double, void*),
			double x_lo, double x_hi, void* aEquationPars = 0,
			const double relError=1e-3, const int max_iter=100,
			const gsl_root_fsolver_type *T = gsl_root_fsolver_bisection);

	static double SolveEquation(
			gsl_function f,	double x_lo, double x_hi,
			const double relError=1e-3, const int max_iter=100,
			const gsl_root_fsolver_type *T = gsl_root_fsolver_bisection);

	double Integration_qag (
			gsl_function aFunction,
			double aXmin,
			double aXmax,
			double epsabs,
			double epsrel,
			size_t limit,
			int key=GSL_INTEG_GAUSS15);
	template<typename X> bool SampleLogscaleDistribution(const Function& aDistrib, double aRand, X& aOutputX, X& aOutputIntegral, int nStepsS, X xMin, X xMax, double aRelError);
	static bool SampleDistribution(const Function& aDistrib, double aRand, double& aOutputX, double& aOutputIntegral, double xMin, double xMax, double aRelError);
	static bool SampleLogDistribution(const Function& aDistrib, double aRand, double& aOutputX, double& aOutputIntegral, double xMin, double xMax, double aRelError);

	static bool SampleLogDistributionBoost(const Function& aDistrib, double aRand, double& aOutputX, double& aOutputIntegral, double xMin, double xMax, double aRelError);
	static bool SampleLogDistributionNR(const Function& aDistrib, double aRand, double& aOutputX, double& aOutputIntegral, double xMin, double xMax, double aRelError);

	template<typename X> static void RelAccuracy(X& aOutput);
	static void SetLogger(IFunctionCallHandlerX<double>* aLogger) { fLogger = aLogger; }
	static int UnitTest();
private:
	static double GslProxySampleLogscaleDistributionFunc (double x, void * params)
	{
		x=exp(x);
		const Function* f = (const Function*)params;
		return x*f->f(x);
	}
	gsl_integration_workspace* gslQAGintegrator;
	static IFunctionCallHandlerX<double>* fLogger;
};

class GslProxyFunction : public Function
{
public:
	GslProxyFunction(gsl_function aFunction, double aXmin=-DBL_MAX, double aXmax=DBL_MAX) :
		fFunction(aFunction),
		fXmin(aXmin),
		fXmax(aXmax){};
	virtual double f(double _x) const
	{
		return fFunction.function(_x,fFunction.params);
	}
	virtual double Xmin() const {return fXmin;}
	virtual double Xmax() const {return fXmax;}
	Function* Clone() const
	{//There is no way to clone fFunction.params
		ASSERT(0);
		Exception::Throw("GslProxyFunction doesn't support Clone");
		return 0;
	}
private:
	gsl_function fFunction;
	double fXmin;
	double fXmax;
};

template<typename X = double >
class ConstFunctionX : public FunctionX<X>
{
public:
	ConstFunctionX(X aValue, X aXmin=-std::numeric_limits<X>::max(), X aXmax=std::numeric_limits<X>::max()) :
		fValue(aValue),
		fXmin(aXmin),
		fXmax(aXmax){};
	virtual X f(X _x) const
	{
		return fValue;
	}
	virtual X Xmin() const {return fXmin;}
	virtual X Xmax() const {return fXmax;}
	FunctionX<X>* Clone() const
	{
		return new ConstFunctionX<X>(fValue, fXmin, fXmax);
	}
private:
	X fValue;
	X fXmin;
	X fXmax;
};

typedef ConstFunctionX<double> ConstFunction;

} /* namespace Utils */
#endif /* MATHUTILS_H_ */

max_stars_repo_name: alexkorochkin/mcray

max_issues_repo_name: alexkorochkin/mcray

max_forks_repo_name: alexkorochkin/mcray


text:
/* multimin/simplex2.c
 *
 * Copyright (C) 2007, 2008, 2009 Brian Gough
 * Copyright (C) 2002 Tuomo Keskitalo, Ivo Alxneit
 *
 * 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/*
   - Originally written by Tuomo Keskitalo <tuomo.keskitalo@iki.fi>
   - Corrections to nmsimplex_iterate and other functions
     by Ivo Alxneit <ivo.alxneit@psi.ch>
   - Additional help by Brian Gough <bjg@network-theory.co.uk>
   - Optimisations added by Brian Gough <bjg@network-theory.co.uk>
         + use BLAS for frequently-called functions
         + keep track of the center to avoid unnecessary computation
         + compute size as RMS value, allowing linear update on each step
           instead of recomputing from all N+1 vectors.
*/

/* The Simplex method of Nelder and Mead, also known as the polytope
   search alogorithm.  Ref: Nelder, J.A., Mead, R., Computer Journal 7
   (1965) pp. 308-313.

   This implementation uses n+1 corner points in the simplex.
*/

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_matrix_double.h>

typedef struct
{
  gsl_matrix *x1;		/* simplex corner points */
  gsl_vector *y1;		/* function value at corner points */
  gsl_vector *ws1;		/* workspace 1 for algorithm */
  gsl_vector *ws2;		/* workspace 2 for algorithm */
  gsl_vector *center;		/* center of all points */
  gsl_vector *delta;		/* current step */
  gsl_vector *xmc;		/* x - center (workspace) */
  double S2;
  unsigned long count;
}
nmsimplex_state_t;

static int
compute_center (const nmsimplex_state_t * state, gsl_vector * center);
static double
compute_size (nmsimplex_state_t * state, const gsl_vector * center);

static double
try_corner_move (const double coeff,
		 const nmsimplex_state_t * state,
		 size_t corner,
		 gsl_vector * xc, const gsl_multimin_function * f)
{
  /* moves a simplex corner scaled by coeff (negative value represents
     mirroring by the middle point of the "other" corner points)
     and gives new corner in xc and function value at xc as a
     return value
   */

  gsl_matrix *x1 = state->x1;
  const size_t P = x1->size1;
  double newval;

  /* xc = (1-coeff)*(P/(P-1)) * center(all) + ((P*coeff-1)/(P-1))*x_corner */
  {
    double alpha = (1 - coeff) * P / (P - 1.0);
    double beta = (P * coeff - 1.0) / (P - 1.0);
    gsl_vector_const_view row = gsl_matrix_const_row (x1, corner);

    gsl_vector_memcpy (xc, state->center);
    gsl_blas_dscal (alpha, xc);
    gsl_blas_daxpy (beta, &row.vector, xc);
  }

  newval = GSL_MULTIMIN_FN_EVAL (f, xc);

  return newval;
}

static void
update_point (nmsimplex_state_t * state, size_t i,
	      const gsl_vector * x, double val)
{
  gsl_vector_const_view x_orig = gsl_matrix_const_row (state->x1, i);
  const size_t P = state->x1->size1;

  /* Compute delta = x - x_orig */
  gsl_vector_memcpy (state->delta, x);
  gsl_blas_daxpy (-1.0, &x_orig.vector, state->delta);

  /* Compute xmc = x_orig - c */
  gsl_vector_memcpy (state->xmc, &x_orig.vector);
  gsl_blas_daxpy (-1.0, state->center, state->xmc);

  /* Update size: S2' = S2 + (2/P) * (x_orig - c).delta + (P-1)*(delta/P)^2 */
  {
    double d = gsl_blas_dnrm2 (state->delta);
    double xmcd;
    gsl_blas_ddot (state->xmc, state->delta, &xmcd);
    state->S2 += (2.0 / P) * xmcd + ((P - 1.0) / P) * (d * d / P);
  }

  /* Update center:  c' = c + (x - x_orig) / P */

  {
    double alpha = 1.0 / P;
    gsl_blas_daxpy (-alpha, &x_orig.vector, state->center);
    gsl_blas_daxpy (alpha, x, state->center);
  }

  gsl_matrix_set_row (state->x1, i, x);
  gsl_vector_set (state->y1, i, val);
}

static int
contract_by_best (nmsimplex_state_t * state, size_t best,
		  gsl_vector * xc, gsl_multimin_function * f)
{

  /* Function contracts the simplex in respect to best valued
     corner. That is, all corners besides the best corner are moved.
     (This function is rarely called in practice, since it is the last
     choice, hence not optimised - BJG)  */

  /* the xc vector is simply work space here */

  gsl_matrix *x1 = state->x1;
  gsl_vector *y1 = state->y1;

  size_t i, j;
  double newval;

  int status = GSL_SUCCESS;

  for (i = 0; i < x1->size1; i++)
    {
      if (i != best)
	{
	  for (j = 0; j < x1->size2; j++)
	    {
	      newval = 0.5 * (gsl_matrix_get (x1, i, j)
			      + gsl_matrix_get (x1, best, j));
	      gsl_matrix_set (x1, i, j, newval);
	    }

	  /* evaluate function in the new point */

	  gsl_matrix_get_row (xc, x1, i);
	  newval = GSL_MULTIMIN_FN_EVAL (f, xc);
	  gsl_vector_set (y1, i, newval);

	  /* notify caller that we found at least one bad function value.
	     we finish the contraction (and do not abort) to allow the user
	     to handle the situation */

	  if (!gsl_finite (newval))
	    {
	      status = GSL_EBADFUNC;
	    }
	}
    }

  /* We need to update the centre and size as well */
  compute_center (state, state->center);
  compute_size (state, state->center);

  return status;
}

static int
compute_center (const nmsimplex_state_t * state, gsl_vector * center)
{
  /* calculates the center of the simplex and stores in center */

  gsl_matrix *x1 = state->x1;
  const size_t P = x1->size1;
  size_t i;

  gsl_vector_set_zero (center);

  for (i = 0; i < P; i++)
    {
      gsl_vector_const_view row = gsl_matrix_const_row (x1, i);
      gsl_blas_daxpy (1.0, &row.vector, center);
    }

  {
    const double alpha = 1.0 / P;
    gsl_blas_dscal (alpha, center);
  }

  return GSL_SUCCESS;
}

static double
compute_size (nmsimplex_state_t * state, const gsl_vector * center)
{
  /* calculates simplex size as rms sum of length of vectors
     from simplex center to corner points:

     sqrt( sum ( || y - y_middlepoint ||^2 ) / n )
   */

  gsl_vector *s = state->ws1;
  gsl_matrix *x1 = state->x1;
  const size_t P = x1->size1;
  size_t i;

  double ss = 0.0;

  for (i = 0; i < P; i++)
    {
      double t;
      gsl_matrix_get_row (s, x1, i);
      gsl_blas_daxpy (-1.0, center, s);
      t = gsl_blas_dnrm2 (s);
      ss += t * t;
    }

  /* Store squared size in the state */
  state->S2 = (ss / P);

  return sqrt (ss / P);
}

static int
nmsimplex_alloc (void *vstate, size_t n)
{
  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  if (n == 0)
    {
      GSL_ERROR ("invalid number of parameters specified", GSL_EINVAL);
    }

  state->x1 = gsl_matrix_alloc (n + 1, n);

  if (state->x1 == NULL)
    {
      GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
    }

  state->y1 = gsl_vector_alloc (n + 1);

  if (state->y1 == NULL)
    {
      gsl_matrix_free (state->x1);
      GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM);
    }

  state->ws1 = gsl_vector_alloc (n);

  if (state->ws1 == NULL)
    {
      gsl_matrix_free (state->x1);
      gsl_vector_free (state->y1);
      GSL_ERROR ("failed to allocate space for ws1", GSL_ENOMEM);
    }

  state->ws2 = gsl_vector_alloc (n);

  if (state->ws2 == NULL)
    {
      gsl_matrix_free (state->x1);
      gsl_vector_free (state->y1);
      gsl_vector_free (state->ws1);
      GSL_ERROR ("failed to allocate space for ws2", GSL_ENOMEM);
    }

  state->center = gsl_vector_alloc (n);

  if (state->center == NULL)
    {
      gsl_matrix_free (state->x1);
      gsl_vector_free (state->y1);
      gsl_vector_free (state->ws1);
      gsl_vector_free (state->ws2);
      GSL_ERROR ("failed to allocate space for center", GSL_ENOMEM);
    }

  state->delta = gsl_vector_alloc (n);

  if (state->delta == NULL)
    {
      gsl_matrix_free (state->x1);
      gsl_vector_free (state->y1);
      gsl_vector_free (state->ws1);
      gsl_vector_free (state->ws2);
      gsl_vector_free (state->center);
      GSL_ERROR ("failed to allocate space for delta", GSL_ENOMEM);
    }

  state->xmc = gsl_vector_alloc (n);

  if (state->xmc == NULL)
    {
      gsl_matrix_free (state->x1);
      gsl_vector_free (state->y1);
      gsl_vector_free (state->ws1);
      gsl_vector_free (state->ws2);
      gsl_vector_free (state->center);
      gsl_vector_free (state->delta);
      GSL_ERROR ("failed to allocate space for xmc", GSL_ENOMEM);
    }

  state->count = 0;

  return GSL_SUCCESS;
}

static void
nmsimplex_free (void *vstate)
{
  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  gsl_matrix_free (state->x1);
  gsl_vector_free (state->y1);
  gsl_vector_free (state->ws1);
  gsl_vector_free (state->ws2);
  gsl_vector_free (state->center);
  gsl_vector_free (state->delta);
  gsl_vector_free (state->xmc);
}

static int
nmsimplex_set (void *vstate, gsl_multimin_function * f,
	       const gsl_vector * x,
	       double *size, const gsl_vector * step_size)
{
  int status;
  size_t i;
  double val;

  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  gsl_vector *xtemp = state->ws1;

  if (xtemp->size != x->size)
    {
      GSL_ERROR ("incompatible size of x", GSL_EINVAL);
    }

  if (xtemp->size != step_size->size)
    {
      GSL_ERROR ("incompatible size of step_size", GSL_EINVAL);
    }

  /* first point is the original x0 */

  val = GSL_MULTIMIN_FN_EVAL (f, x);

  if (!gsl_finite (val))
    {
      GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
    }

  gsl_matrix_set_row (state->x1, 0, x);
  gsl_vector_set (state->y1, 0, val);

  /* following points are initialized to x0 + step_size */

  for (i = 0; i < x->size; i++)
    {
      status = gsl_vector_memcpy (xtemp, x);

      if (status != 0)
	{
	  GSL_ERROR ("vector memcopy failed", GSL_EFAILED);
	}

      {
	double xi = gsl_vector_get (x, i);
	double si = gsl_vector_get (step_size, i);

	gsl_vector_set (xtemp, i, xi + si);
	val = GSL_MULTIMIN_FN_EVAL (f, xtemp);
      }

      if (!gsl_finite (val))
	{
	  GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
	}

      gsl_matrix_set_row (state->x1, i + 1, xtemp);
      gsl_vector_set (state->y1, i + 1, val);
    }

  compute_center (state, state->center);

  /* Initialize simplex size */
  *size = compute_size (state, state->center);

  state->count++;

  return GSL_SUCCESS;
}

static int
nmsimplex_iterate (void *vstate, gsl_multimin_function * f,
		   gsl_vector * x, double *size, double *fval)
{

  /* Simplex iteration tries to minimize function f value */
  /* Includes corrections from Ivo Alxneit <ivo.alxneit@psi.ch> */

  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  /* xc and xc2 vectors store tried corner point coordinates */

  gsl_vector *xc = state->ws1;
  gsl_vector *xc2 = state->ws2;
  gsl_vector *y1 = state->y1;
  gsl_matrix *x1 = state->x1;

  const size_t n = y1->size;
  size_t i;
  size_t hi, s_hi, lo;
  double dhi, ds_hi, dlo;
  int status;
  double val, val2;

  if (xc->size != x->size)
    {
      GSL_ERROR ("incompatible size of x", GSL_EINVAL);
    }

  /* get index of highest, second highest and lowest point */

  dhi = dlo = gsl_vector_get (y1, 0);
  hi = 0;
  lo = 0;

  ds_hi = gsl_vector_get (y1, 1);
  s_hi = 1;

  for (i = 1; i < n; i++)
    {
      val = (gsl_vector_get (y1, i));
      if (val < dlo)
	{
	  dlo = val;
	  lo = i;
	}
      else if (val > dhi)
	{
	  ds_hi = dhi;
	  s_hi = hi;
	  dhi = val;
	  hi = i;
	}
      else if (val > ds_hi)
	{
	  ds_hi = val;
	  s_hi = i;
	}
    }

  /* try reflecting the highest value point */

  val = try_corner_move (-1.0, state, hi, xc, f);

  if (gsl_finite (val) && val < gsl_vector_get (y1, lo))
    {
      /* reflected point is lowest, try expansion */

      val2 = try_corner_move (-2.0, state, hi, xc2, f);

      if (gsl_finite (val2) && val2 < gsl_vector_get (y1, lo))
	{
	  update_point (state, hi, xc2, val2);
	}
      else
	{
	  update_point (state, hi, xc, val);
	}
    }
  else if (!gsl_finite (val) || val > gsl_vector_get (y1, s_hi))
    {
      /* reflection does not improve things enough, or we got a
         non-finite function value */

      if (gsl_finite (val) && val <= gsl_vector_get (y1, hi))
	{
	  /* if trial point is better than highest point, replace
	     highest point */

	  update_point (state, hi, xc, val);
	}

      /* try one-dimensional contraction */

      val2 = try_corner_move (0.5, state, hi, xc2, f);

      if (gsl_finite (val2) && val2 <= gsl_vector_get (y1, hi))
	{
	  update_point (state, hi, xc2, val2);
	}
      else
	{
	  /* contract the whole simplex about the best point */

	  status = contract_by_best (state, lo, xc, f);

	  if (status != GSL_SUCCESS)
	    {
	      GSL_ERROR ("contraction failed", GSL_EFAILED);
	    }
	}
    }
  else
    {
      /* trial point is better than second highest point.  Replace
         highest point by it */

      update_point (state, hi, xc, val);
    }

  /* return lowest point of simplex as x */

  lo = gsl_vector_min_index (y1);
  gsl_matrix_get_row (x, x1, lo);
  *fval = gsl_vector_get (y1, lo);

  /* Update simplex size */

  {
    double S2 = state->S2;

    if (S2 > 0)
      {
	*size = sqrt (S2);
      }
    else
      {
	/* recompute if accumulated error has made size invalid */
	*size = compute_size (state, state->center);
      }
  }

  return GSL_SUCCESS;
}

static const gsl_multimin_fminimizer_type nmsimplex_type =
{ "nmsimplex2",	/* name */
  sizeof (nmsimplex_state_t),
  &nmsimplex_alloc,
  &nmsimplex_set,
  &nmsimplex_iterate,
  &nmsimplex_free
};

const gsl_multimin_fminimizer_type
  * gsl_multimin_fminimizer_nmsimplex2 = &nmsimplex_type;

static inline double
ran_unif (unsigned long *seed)
{
  unsigned long s = *seed;
  *seed = (s * 69069 + 1) & 0xffffffffUL;
  return (*seed) / 4294967296.0;
}

static int
nmsimplex_set_rand (void *vstate, gsl_multimin_function * f,
		    const gsl_vector * x,
		    double *size, const gsl_vector * step_size)
{
  size_t i, j;
  double val;

  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  gsl_vector *xtemp = state->ws1;

  if (xtemp->size != x->size)
    {
      GSL_ERROR ("incompatible size of x", GSL_EINVAL);
    }

  if (xtemp->size != step_size->size)
    {
      GSL_ERROR ("incompatible size of step_size", GSL_EINVAL);
    }

  /* first point is the original x0 */

  val = GSL_MULTIMIN_FN_EVAL (f, x);

  if (!gsl_finite (val))
    {
      GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
    }

  gsl_matrix_set_row (state->x1, 0, x);
  gsl_vector_set (state->y1, 0, val);

  {
    gsl_matrix_view m =
      gsl_matrix_submatrix (state->x1, 1, 0, x->size, x->size);

    /* generate a random orthornomal basis  */
    unsigned long seed = state->count ^ 0x12345678;

    ran_unif (&seed);		/* warm it up */

    gsl_matrix_set_identity (&m.matrix);

    /* start with random reflections */
    for (i = 0; i < x->size; i++)
      {
        double s = ran_unif (&seed);
        if (s > 0.5) gsl_matrix_set (&m.matrix, i, i, -1.0);
      }

    /* apply random rotations */
    for (i = 0; i < x->size; i++)
      {
	for (j = i + 1; j < x->size; j++)
	  {
	    /* rotate columns i and j by a random angle */
	    double angle = 2.0 * M_PI * ran_unif (&seed);
	    double c = cos (angle), s = sin (angle);
	    gsl_vector_view c_i = gsl_matrix_column (&m.matrix, i);
	    gsl_vector_view c_j = gsl_matrix_column (&m.matrix, j);
	    gsl_blas_drot (&c_i.vector, &c_j.vector, c, s);
	  }
      }

    /* scale the orthonormal basis by the user-supplied step_size in
       each dimension, and use as an offset from the central point x */

    for (i = 0; i < x->size; i++)
      {
	double x_i = gsl_vector_get (x, i);
	double s_i = gsl_vector_get (step_size, i);
	gsl_vector_view c_i = gsl_matrix_column (&m.matrix, i);

	for (j = 0; j < x->size; j++)
	  {
	    double x_ij = gsl_vector_get (&c_i.vector, j);
	    gsl_vector_set (&c_i.vector, j, x_i + s_i * x_ij);
	  }
      }

    /* compute the function values at each offset point */

    for (i = 0; i < x->size; i++)
      {
	gsl_vector_view r_i = gsl_matrix_row (&m.matrix, i);

	val = GSL_MULTIMIN_FN_EVAL (f, &r_i.vector);

	if (!gsl_finite (val))
	  {
	    GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
	  }

	gsl_vector_set (state->y1, i + 1, val);
      }
  }

  compute_center (state, state->center);

  /* Initialize simplex size */
  *size = compute_size (state, state->center);

  state->count++;

  return GSL_SUCCESS;
}

static const gsl_multimin_fminimizer_type nmsimplex2rand_type =
{ "nmsimplex2rand",	/* name */
  sizeof (nmsimplex_state_t),
  &nmsimplex_alloc,
  &nmsimplex_set_rand,
  &nmsimplex_iterate,
  &nmsimplex_free
};

const gsl_multimin_fminimizer_type
  * gsl_multimin_fminimizer_nmsimplex2rand = &nmsimplex2rand_type;

max_stars_repo_name: ielomariala/Hex-Game


text:
#ifndef LIB_INCLUDE_TICK_ARRAY_VECTOR_OPERATIONS_H_
#define LIB_INCLUDE_TICK_ARRAY_VECTOR_OPERATIONS_H_

// License: BSD 3 clause

#include <atomic>
#include <vector>
#include <numeric>
#include <algorithm>
#include <type_traits>

#include "promote.h"
#include "tick/base/defs.h"

#if defined(TICK_USE_MKL)

#include "mkl.h"

#elif defined(TICK_USE_CBLAS)

#if defined(__APPLE__)
#include <Accelerate/Accelerate.h>
// TODO(svp) Disabling this feature until we find
//  a good way to determine if ATLAS is actually available
#else
extern "C" {
#include <cblas.h>
}
#endif  // defined(__APPLE__)

#else

#include "tick/array/vector/ops_unoptimized.h"
namespace tick {
template <typename T>
using vector_operations = detail::vector_operations_unoptimized<T>;
  }

#endif

#if defined(TICK_USE_MKL) || defined(TICK_USE_CBLAS)
#include "tick/array/vector/ops_blas.h"
namespace tick {
template <typename T>
using vector_operations = detail::vector_operations_cblas<T>;
  }
#endif

#include "tick/array/vector/ops_unoptimized_impl.h"

#endif  // LIB_INCLUDE_TICK_ARRAY_VECTOR_OPERATIONS_H_

max_stars_repo_name: sumau/tick

max_issues_repo_name: saurabhdash/tick

max_forks_repo_name: saurabhdash/tick


text:
#ifndef __MATH_LOGSPACE_H__
#define __MATH_LOGSPACE_H__
#include <cmath>
#include <vector>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_gamma.h>
#include "logspace_base.h"
#include "specialfunc.h"

// Given two log vectors, log a_i and log b_i, compute
// log sum (a_i * b_i).
double log_dot_product(const gsl_vector* log_a, const gsl_vector* log_b);

// Given a log vector, log a_i, compute log sum a_i.  Returns the sum.
double log_normalize(gsl_vector* x);

// Compute the log sum over all elements in the vector
double log_sum(const gsl_vector* x);

// Given a log matrix, log a_i, compute log sum a_i.  Returns the sum.
double log_normalize_matrix(gsl_matrix* x);

double log_dirichlet_likelihood(const double sum,
                                const double prior_sum,
                                const std::vector<int>& counts,
                                bool debug = false);

double log_dirichlet_likelihood(const double sum,
                                const double prior_scale,
                                const gsl_vector* prior,
                                const std::vector<int>& counts);

#endif  // __MATH_LOGSPACE_H__

max_stars_repo_name: boomsbloom/dtm-fmri

max_issues_repo_name: boomsbloom/dtm-fmri

max_forks_repo_name: boomsbloom/dtm-fmri