text:
/* multifit/gcv.c
 *
 * Copyright (C) 2016 Patrick Alken
 *
 * 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.
 */

/*
 * References:
 *
 * [1] P. C. Hansen, "Discrete Inverse Problems: Insight and Algorithms,"
 * SIAM Press, 2010.
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_min.h>

typedef struct
{
  const gsl_vector * S;
  const gsl_vector * UTy;
  double delta0;
  size_t np;
  gsl_vector * workp;
} gcv_params;

static double gcv_func(double lambda, void * params);

/*
gsl_multifit_linear_gcv_init()
  Initialize Generalized Cross Validation parameters

Inputs: y         - right hand side vector
        reg_param - (output) regularization parameters
        UTy       - (output) U^T y
        delta0    - (output) delta0
        work      - workspace
*/

int
gsl_multifit_linear_gcv_init(const gsl_vector * y,
                             gsl_vector * reg_param,
                             gsl_vector * UTy,
                             double * delta0,
                             gsl_multifit_linear_workspace * work)
{
  const size_t n = y->size;

  if (n != work->n)
    {
      GSL_ERROR("y vector does not match workspace", GSL_EBADLEN);
    }
  else if (UTy->size != work->p)
    {
      GSL_ERROR ("UTy vector does not match workspace", GSL_EBADLEN);
    }
  else
    {
      const size_t p = work->p;

      gsl_matrix_view U = gsl_matrix_submatrix(work->A, 0, 0, n, p);
      gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);

      const double smax = gsl_vector_get(&S.vector, 0);
      const double smin = gsl_vector_get(&S.vector, p - 1);

      double dr; /* residual error from projection */

      double normy = gsl_blas_dnrm2(y);
      double normUTy;

      /* compute projection UTy = U^T y */
      gsl_blas_dgemv (CblasTrans, 1.0, &U.matrix, y, 0.0, UTy);
      normUTy = gsl_blas_dnrm2(UTy);

      /* dr = ||y||^2 - ||U^T y||^2 */
      dr = (normy + normUTy) * (normy - normUTy);

      /* calculate regularization parameters */
      gsl_multifit_linear_lreg(smin, smax, reg_param);

      if (n > p && dr > 0.0)
        *delta0 = dr;
      else
        *delta0 = 0.0;

      return GSL_SUCCESS;
    }
}

/*
gsl_multifit_linear_gcv_curve()
  Calculate Generalized Cross Validation curve for a set
of regularization parameters

Inputs: reg_param - regularization parameters
        UTy       - U^T y vector, size p
        delta0    - delta0
        G         - (output) GCV curve values
        work      - workspace
*/

int
gsl_multifit_linear_gcv_curve(const gsl_vector * reg_param,
                              const gsl_vector * UTy,
                              const double delta0,
                              gsl_vector * G,
                              gsl_multifit_linear_workspace * work)
{
  const size_t n = work->n;
  const size_t p = work->p;
  const size_t N = reg_param->size; /* number of points on GCV curve */

  if (UTy->size != p)
    {
      GSL_ERROR("UTy vector does not match workspace", GSL_EBADLEN);
    }
  else if (G->size != N)
    {
      GSL_ERROR ("size of reg_param and G vectors do not match",
                 GSL_EBADLEN);
    }
  else
    {
      size_t i;

      gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
      gsl_vector_view workp = gsl_matrix_subcolumn(work->QSI, 0, 0, p);

      gcv_params params;

      params.S = &S.vector;
      params.UTy = UTy;
      params.delta0 = delta0;
      params.np = n - p;
      params.workp = &workp.vector;

      for (i = 0; i < N; ++i)
        {
          double lambdai = gsl_vector_get(reg_param, i);
          double Gi = gcv_func(lambdai, &params);

          gsl_vector_set(G, i, Gi);
        }

      return GSL_SUCCESS;
    }
}

/*
gsl_multifit_linear_gcv_min()
  Find regularization parameter which minimizes GCV curve

Inputs: reg_param - regularization parameters
        UTy       - U^T y vector, size p
        G         - GCV curve values
        delta0    - delta0
        lambda    - (output) optimal regularization parameter
        work      - workspace
*/

int
gsl_multifit_linear_gcv_min(const gsl_vector * reg_param,
                            const gsl_vector * UTy,
                            const gsl_vector * G,
                            const double delta0,
                            double * lambda,
                            gsl_multifit_linear_workspace * work)
{
  const size_t n = work->n;
  const size_t p = work->p;
  const size_t npts = reg_param->size; /* number of points on GCV curve */

  if (UTy->size != p)
    {
      GSL_ERROR("UTy vector does not match workspace", GSL_EBADLEN);
    }
  else if (G->size != npts)
    {
      GSL_ERROR ("size of reg_param and G vectors do not match",
                 GSL_EBADLEN);
    }
  else
    {
      int status;
      const size_t max_iter = 500;
      const double tol = 1.0e-4;
      gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
      gsl_vector_view workp = gsl_matrix_subcolumn(work->QSI, 0, 0, p);
      gcv_params params;
      int idxG = (int) gsl_vector_min_index(G);
      double a = gsl_vector_get(reg_param, GSL_MIN(idxG + 1, (int) npts - 1));
      double b = gsl_vector_get(reg_param, GSL_MAX(idxG - 1, 0));
      double m = gsl_vector_get(reg_param, idxG);
      size_t iter = 0;
      gsl_function F;

      /* XXX FIXME */
      gsl_min_fminimizer *min_workspace_p;

      if (idxG == 0 || idxG == ((int)npts - 1))
        {
          /* the minimum is an endpoint of the curve, no need to search */
          *lambda = m;
          return GSL_SUCCESS;
        }

      /* XXX FIXME */
      min_workspace_p = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);

      params.S = &S.vector;
      params.UTy = UTy;
      params.delta0 = delta0;
      params.np = n - p;
      params.workp = &workp.vector;

      F.function = gcv_func;
      F.params = &params;

      gsl_min_fminimizer_set(min_workspace_p, &F, m, a, b);

      do
        {
          iter++;
          status = gsl_min_fminimizer_iterate(min_workspace_p);

          a = gsl_min_fminimizer_x_lower(min_workspace_p);
          b = gsl_min_fminimizer_x_upper(min_workspace_p);

          status = gsl_min_test_interval(a, b, 0.0, tol);
        }
      while (status == GSL_CONTINUE && iter < max_iter);

      if (status == GSL_SUCCESS)
        *lambda = gsl_min_fminimizer_minimum(min_workspace_p);
      else
        status = GSL_EMAXITER;

      gsl_min_fminimizer_free(min_workspace_p);

      return status;
    }
}

/*
gsl_multifit_linear_gcv_calc()
  Calculate GCV function G(lambda) for given lambda

Inputs: reg_param - regularization parameters
        UTy       - U^T y vector, size p
        delta0    - delta0
        G         - (output) GCV curve values
        work      - workspace
*/

double
gsl_multifit_linear_gcv_calc(const double lambda,
                             const gsl_vector * UTy,
                             const double delta0,
                             gsl_multifit_linear_workspace * work)
{
  const size_t n = work->n;
  const size_t p = work->p;

  if (UTy->size != p)
    {
      GSL_ERROR_VAL("UTy vector does not match workspace", GSL_EBADLEN, 0.0);
    }
  else
    {
      gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
      gsl_vector_view workp = gsl_matrix_subcolumn(work->QSI, 0, 0, p);
      gcv_params params;
      double G;

      params.S = &S.vector;
      params.UTy = UTy;
      params.delta0 = delta0;
      params.np = n - p;
      params.workp = &workp.vector;

      G = gcv_func(lambda, &params);

      return G;
    }
}

/*
gsl_multifit_linear_gcv()
  Calculate Generalized Cross Validation curve for a set
of regularization parameters

Inputs: y         - right hand side vector
        reg_param - (output) regularization parameters
        G         - (output) GCV curve values
        lambda    - (output) optimal regularization parameter which
                    minimizes GCV curve
        G_lambda  - (output) G(lambda) value at optimal parameter
        work      - workspace
*/

int
gsl_multifit_linear_gcv(const gsl_vector * y,
                        gsl_vector * reg_param,
                        gsl_vector * G,
                        double * lambda,
                        double * G_lambda,
                        gsl_multifit_linear_workspace * work)
{
  const size_t n = y->size;
  const size_t N = G->size; /* number of points on GCV curve */

  if (n != work->n)
    {
      GSL_ERROR("y vector does not match workspace", GSL_EBADLEN);
    }
  else if (reg_param->size != N)
    {
      GSL_ERROR ("size of reg_param and G vectors do not match",
                 GSL_EBADLEN);
    }
  else
    {
      int status;
      const size_t p = work->p;
      gsl_vector_view UTy = gsl_vector_subvector(work->xt, 0, p);
      double delta0;

      status = gsl_multifit_linear_gcv_init(y, reg_param, &UTy.vector, &delta0, work);
      if (status)
        return status;

      status = gsl_multifit_linear_gcv_curve(reg_param, &UTy.vector, delta0, G, work);
      if (status)
        return status;

      status = gsl_multifit_linear_gcv_min(reg_param, &UTy.vector, G, delta0, lambda, work);
      if (status)
        return status;

      *G_lambda = gsl_multifit_linear_gcv_calc(*lambda, &UTy.vector, delta0, work);

      return GSL_SUCCESS;
    }
}

static double
gcv_func(double lambda, void * params)
{
  gcv_params * par = (gcv_params *) params;
  const gsl_vector *S = par->S;
  const gsl_vector *UTy = par->UTy;
  double delta0 = par->delta0;
  size_t np = par->np;
  gsl_vector *workp = par->workp;
  const size_t p = S->size;
  size_t i;
  double lambda_sq = lambda * lambda;
  double G, d, norm;
  double sumf = 0.0;

  /* compute workp = 1 - filter_factors */
  for (i = 0; i < p; ++i)
    {
      double si = gsl_vector_get(S, i);
      double fi = lambda_sq / (si * si + lambda_sq);
      gsl_vector_set(workp, i, fi);
      sumf += fi;
    }

  d = (double)np + sumf;

  gsl_vector_mul(workp, UTy);
  norm = gsl_blas_dnrm2(workp);

  G = (norm*norm + delta0) / (d * d);

  return G;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef _MATH_FUNCTIONS_H_
#define _MATH_FUNCTIONS_H_

extern "C" {
#include <cblas.h>
}
#include <math.h>
#include <xmmintrin.h>
//#include <omp.h>
#include <cublas_v2.h>
#include "arena.h"

void cpu_rearrange(int batch_size, int input_dim, int num_blocks,
                  const float *T, float *X);
void gpu_rearrange(int batcch_size, int input_dim, int num_blocks,
                  const float *T, float *X);

void svd2x2(const float a[4], float u[4], float s[2], float v[4]);
void matmulcomplex2x2(const float A[8], const float B[8], float C[8]);
void conjugate_transpose2x2(float M[8]);
void conjugate_transpose2x2(const float M[8], float R[8]);
void cpu_complex_modulus_interleaved_relu(int M, int N, const float *X,
                                          const float* b, float* Y);

void cpu_permute_forward(int N,  const float*X, const int* PI,
                         float* Y);

void cpu_permute_backward(int N,  const float*X, const int* PI,
                          float* Y);

void cpu_diag_gemm(int M, int N, float alpha,  const float* X,
                   const float* W, float beta, float *Y );
void gpu_diag_gemm(int M, int N, float alpha, const float* X,
                   const float* W, float beta, float *Y );

void cpu_permute(int M, int N, const float* PI,
                const float* X, float* Y);

void gpu_permute(int M, int N, const float* PI,
                const float* X, float* Y);

float cpu_relative_error(int N,  const float* grad_ptr,
                         const float* numer_grad_ptr);

void gpu_real_svd(int M, int N, float *A,
                 float *U, float *S, float *V);
void gpu_unitary_svd(int M, int N, float *A,
                 float *U, float *S, float *V);
void gpu_random_orthogonal_matrix(int N, float *W);
void gpu_random_unitary_matrix(int N, float *W);

void cpu_initialize_svd2x2(int N, float* W);
void cpu_initialize_complexsvd2x2(int N, float* W);
void cpu_initialize_unitary(int NF, int* factor_sizes,  float *W);
void cpu_initialize_orthogonal(int NF, int* factor_sizes,  float *W);
void gpu_initialize_unitary(int NF, int* factor_sizes,  float *W);
void gpu_initialize_orthogonal(int NF, int* factor_sizes,  float *W);
void cpu_initialize_glorot10(int NF, int* factor_sizes,  float *W);
float cpu_complexdet2x2(const float W[8], float D[2]);
void cpu_gemm(const CBLAS_TRANSPOSE TransA,
              const CBLAS_TRANSPOSE TransB, int M, int N, int K,
              const float alpha, const float* A, const float* B, const float beta,
              float* C);
void gpu_gemm(const CBLAS_TRANSPOSE TransA,
              const CBLAS_TRANSPOSE TransB, int M, int N, int K,
              const float alpha, const float* A, const float* B, const float beta,
              float* C);

void cpu_complex_gemm(const CBLAS_TRANSPOSE TransA,
                      const CBLAS_TRANSPOSE TransB, int M, int N, int K,
                      const float alpha, const float* A, const float* B,
                      const float beta, float* C);

void cpu_real_complex_gemm(const CBLAS_TRANSPOSE TransA,
                           const CBLAS_TRANSPOSE TransB, int M, int N, int K,
                           const float alpha, const float* A, const float* B,
                           const float beta, float* C);

void gpu_complex_gemm(const CBLAS_TRANSPOSE TransA,
                      const CBLAS_TRANSPOSE TransB, int M, int N, int K,
                      const float alpha, const float* A, const float* B,
                      const float beta, float* C);

void gpu_real_complex_gemm(const CBLAS_TRANSPOSE TransA,
                           const CBLAS_TRANSPOSE TransB, int M, int N, int K,
                           const float alpha, const float* A, const float* B,
                           const float beta, float* C);

//Y = X + a * I
void gpu_apx(int N, float a,  const float* X, float* Y);

void gpu_complex_apx(int N, float a,  const float* X, float* Y);

void gpu_geam(const CBLAS_TRANSPOSE TransA, const CBLAS_TRANSPOSE TransB,
              int m, int n, const float *alpha, const float *A,
              const float *beta, const float *B, float *C);

void cpu_gemv(const CBLAS_TRANSPOSE TransA, int M, int N,
              const float alpha, const float* A, const float* x, const float beta,
              float* y);
void gpu_gemv(const CBLAS_TRANSPOSE TransA, int M, int N,
              const float alpha, const float* A, const float* x, const float beta,
              float* y);
void cpu_axpy(int N, const float alpha, const float* X,
              float* Y);
void gpu_axpy(int N, const float alpha, const float* X,
              float* Y);
void gpu_complex_axpy(int N, const float alpha, const float* X,
                      float* Y);
void cpu_addscalar(int N, const float alpha, const float* X,
                   float* Y);
void gpu_addscalar(int N, const float alpha, const float* X,
                   float* Y);
void cpu_axpby(int N, const float alpha, const float* X,
               const float beta, float* Y);
void gpu_axpby(int N, const float alpha, const float* X,
               const float beta, float* Y);
void cpu_copy(int N, const float *X, float *Y);
void gpu_copy(int N, const float *X, float *Y);
void cpu_scal(int N, const float alpha, float *X);
float cpu_max(int N, const float *X);
int cpu_max_index(int N, const float *X);
void cpu_exp(int N, const float *X, float *Y);
void gpu_exp(int N, const float *X, float *Y);
void cpu_dexp(int N, const float *in_data, const float *in_diff, float *out);
void cpu_tanh(int N, const float *X, float *Y);
void gpu_tanh(int N, const float *X, float *Y);

void cpu_sin(int N, const float *X, float *Y);
void gpu_sin(int N, const float *X, float *Y);

void cpu_cos(int N, const float *X, float *Y);
void gpu_cos(int N, const float *X, float *Y);

void cpu_dtanh(int N, const float *in_data, const float *in_diff, float *out);
void gpu_dtanh(int N, const float *in_data, const float *in_diff, float *out);

void cpu_sigmoid(int N, const float *X, float *Y);
void gpu_sigmoid(int N, const float *X, float *Y);
void cpu_dsigmoid(int N, const float *in_data, const float *in_diff, float *out);
void gpu_dsigmoid(int N, const float *in_data, const float *in_diff, float *out);

void cpu_hadamard(int M, int N, float *X);
void gpu_hadamard(int M, int N, float *X);

void cpu_complex_modulus_relu(int M, int N, const float  *X, const float* b, float* Y);
void cpu_complex_modulus_drelu(int M, int N, const float  *X,
                               const float* b, const float* gradY, float* gradX,
                               float* gradb);
void cpu_complex_modulus_tanh(int M, int N, const float  *X, const float* b, float* Y);
void cpu_complex_modulus_dtanh(int M, int N, const float  *X,
                               const float* b, const float* gradY, float* gradX,
                               float* gradb);
/*
    From Efficient backprop, Yann Lecun, Leon Bottou 1998
*/
void cpu_lecun98(int N, int fan_in, float *x);
void gpu_lecun98(int N, int fan_in, float *x);
/*
    From Understanding the diffcuilt of traing neural nets, Xavier Glorot, Yoshua Bengio : AISTATS 2010
*/
void cpu_glorot10(int fan_in, int fan_out, float *x);
void gpu_glorot10(int fan_in, int fan_out, float *x);
/*
    From Delving deep into rectifiers, Kaiming He Xiangyu Zhang Shaoqing Ren Jian Sun:  Arxiv 2015
*/

void cpu_he15(int N, int fan_in, float *x);
void cpu_allzero(int N, float *x);
void cpu_normal(int N, const float mean, const float stddev, float *x);
void cpu_uniform(int N, const float r1, const float r2, float *x);
void gpu_uniform(int N, const float r1, const float r2, float *x);
void gpu_he15(int N, int fan_in, float *x);
void gpu_allzero(int N, float *x);
void gpu_normal(int N, const float mean, const float stddev, float *x);
void cpu_log(int N, const float *X, float *Y);
void gpu_log(int N, const float *X, float *Y);
void cpu_sqr(int N, const float *X, float *Y);
void gpu_sqr(int N, const float *X, float *Y);
void cpu_abs(int N, const float *X, float *Y);
void cpu_clip(int N, const float *X, const float p, float *Y);
void gpu_clip(int N, const float *X, const float p, float *Y);
void cpu_pow(int N, const float *X, const float p, float *Y);
void gpu_pow(int N, const float *X, const float p, float *Y);
void gpu_scal(int N, const float alpha, float *X);
float cpu_dot(int n, const float* x, const float* y);
void cpu_add(int N, const float* a, const float* b, float* y);
void gpu_add(int N, const float* a, const float* b, float* y);
void gpu_mul(int N, const float* a, const float* b, float* y);
void cpu_mul(int N, const float* a, const float* b, float* y);
void gpu_div(int N, const float* a, const float* b, float* y);
void cpu_div(int N, const float* a, const float* b, float* y);
void gpu_fill(int N, const float a, float* y);
void cpu_fill(int N, const float a, float* y);
float cpu_max(int N, const float* y);
float cpu_nrm2(int N, const float* y);
float gpu_nrm2(int N, const float* y);
float cpu_min(int N, const float* y);
float gpu_dot(int n, const float* x, const float* y);

void cpu_chi(int N, int dof, float* x);
void gpu_chi(int N, int dof, float* x);

void cpu_rademacher(int N, float* x);
void gpu_rademacher(int N, float* x);

void cpu_permutation_array(int N, float* x);
void gpu_permutation_array(int N, float* x);

inline void cblas_saxpby(int N, const float alpha, const float* X,
                         int incX, const float beta, float* Y,
                         int incY) {
  cblas_sscal(N, beta, Y, incY);
  cblas_saxpy(N, alpha, X, incX, Y, incY);
}

#endif

ext: h

max_stars_repo_name: cijose/RandomFourierFeatures

max_issues_repo_name: cijose/RandomFourierFeatures

max_forks_repo_name: cijose/RandomFourierFeatures


text:
#ifndef __INC_MULTIMODELSTRUCT__
#define __INC_MULTIMODELSTRUCT__

#include "modelstruct.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

struct modelstruct;

/**
 * \struct multimodelstruct
 * handle a model with t-variate output i.e y_m = {y_1, ..., y_t}^{T}
 *
 * this contains some r <= t modelstructs representing the emulators in the PCA projected
 * space
 */
typedef struct multi_modelstruct{
	int nt; // number of observables in OBSERVED space
	int nr; // number of pca-observables kept in the ROTATED space
	int nparams; // number of dimensions in the PARAMETER SPACE
	int nmodel_points;  // number of points in the design
	int cov_fn_index;
	int regression_order;

	/**
	 * the design (rows:nmodel_points) (cols:nparams)
	 */
	gsl_matrix *xmodel;
	/**
	 * the full training set in the original space
	 * (rows:nmodel_points) (cols:nt)
	 */
	gsl_matrix *training_matrix;
	/**
	 * a t length vector of the the mean values of the cols of training_matrix
	 */
	gsl_vector *training_mean;

	/** array of (r) pca decomposed models
	 *
	 * pointers from this back to xmodel, is this bad?
	 */
	modelstruct** pca_model_array;

	/**
	 * the eigenvalues and vectors from the pca decomp (t xt) , not saved for now...
	 */
	//gsl_vector *pca_eigenvalues;
	//gsl_matrix *pca_eigenvectors;

	gsl_vector *pca_evals_r;	/** just the first r evals */
	gsl_matrix *pca_evecs_r;  /** the first r evecs (t x nr)*/

	/**
	 * an (r x nmodel_points) matrix of the pca-transformed training points, these are
	 * used to init each of the r entries of pca_model_array */
	gsl_matrix *pca_zmatrix; /** z <- (1/sqrt(pca_evals_r)) %*% t(pca_evecs_r) %*% (training_matrix - training_mean) */

} multi_modelstruct;

multi_modelstruct* alloc_multimodelstruct(gsl_matrix *xmodel_in, gsl_matrix *training_matrix_in,
																				 int cov_fn_index, int regression_order, double varfrac);

void gen_pca_decomp(multi_modelstruct *m, double vfrac);
void gen_pca_model_array(multi_modelstruct *m);

void dump_multi_modelstruct(FILE* fptr, multi_modelstruct *m );
multi_modelstruct *load_multi_modelstruct(FILE* fptr);

double vector_elt_sum(gsl_vector* vec, int nstop);

void free_multimodelstruct(multi_modelstruct *m);

#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * 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.
 */

#pragma once

#ifndef mesh_7c5c4ae5_0f3e_4a2d_936e_0447239ac3a7_h
#define mesh_7c5c4ae5_0f3e_4a2d_936e_0447239ac3a7_h

#include <gslib/type.h>
#include <gslib/std.h>
#include <ariel/type.h>

__ariel_begin__

struct mesh_uv_set;
struct mesh_tangent_set;

typedef vector<vec3> mesh_points;
typedef vector<vec3> mesh_normals;
typedef vector<vec3> mesh_tangents;
typedef vector<vec3> mesh_colors;
typedef vector<vec2> mesh_uv_coords;
typedef vector<int32> mesh_faces;
typedef vector<mesh_uv_set> mesh_uv_sets;
typedef vector<mesh_tangent_set> mesh_tangent_sets;

struct mesh_uv_set
{
    string              name;
    mesh_uv_coords      coords;
};

struct mesh_tangent_set
{
    string              name;
    mesh_tangents       tangents;
};

class mesh:
    public res_node
{
    friend class mesh_io;

public:
    mesh();
    virtual ~mesh() {}
    virtual res_type get_type() const override { return res_mesh; }
    bool has_indices() const { return !_faces.empty(); }
    bool has_points() const { return !_points.empty(); }
    bool has_normals() const { return !_normals.empty(); }
    bool has_colors() const { return !_vcolors.empty(); }
    bool has_uvs() const { return !_uvsets.empty(); }
    bool has_tangents() const { return !_tangentsets.empty(); }
    int32 num_uv_sets() const { return (int32)_uvsets.size(); }
    int32 num_tangent_sets() const { return (int32)_tangentsets.size(); }
    int32 num_points() const { return (int32)_points.size(); }
    int32 num_indices() const { return (int32)_faces.size(); }
    const mesh_points& get_const_points() const { return _points; }
    const mesh_faces& get_const_indices() const { return _faces; }
    const mesh_normals& get_const_normals() const { return _normals; }
    const mesh_colors& get_const_colors() const { return _vcolors; }

protected:
    mesh_points         _points;
    mesh_normals        _normals;
    mesh_colors         _vcolors;
    mesh_uv_sets        _uvsets;
    mesh_tangent_sets   _tangentsets;
    mesh_faces          _faces;
    matrix              _transforms;

public:
    /*
     * Caution:
     * If there were splited vertices in the mesh, the normals & tangents could be incorrect,
     * in such situations, you should always trust the data from exporter.
     */
    void calc_normals();
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <math.h>
#include <stdio.h>
#include <iostream>
#include <time.h>

#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include "gsl/gsl_cdf.h"
#include "gsl/gsl_randist.h"
//
#include "mex.h"

max_stars_repo_name: ferjorosa/test-glfm

max_issues_repo_name: ferjorosa/test-glfm

max_forks_repo_name: ferjorosa/test-glfm


text:
/* statistics/test_nist.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Jim Davies, 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.
 */

/* These tests are based on the NIST Statistical Reference Datasets
   See http://www.nist.gov/itl/div898/strd/index.html for more
   information. */

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_test.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_ieee_utils.h>

int
test_nist (void)
{
  size_t i ;

  const size_t nlew = 200 ;

  const double lew[200] = {
    -213, -564,  -35,  -15,  141,  115, -420, -360,  203, -338, -431,  194,
    -220, -513,  154, -125, -559,   92,  -21, -579,  -52,   99, -543, -175,
     162, -457, -346,  204, -300, -474,  164, -107, -572,   -8,   83, -541,
    -224,  180, -420, -374,  201, -236, -531,   83,   27, -564, -112,  131,
    -507, -254,  199, -311, -495,  143,  -46, -579,  -90,  136, -472, -338,
     202, -287, -477,  169, -124, -568,   17,   48, -568, -135,  162, -430,
    -422,  172,  -74, -577,  -13,   92, -534, -243,  194, -355, -465,  156,
     -81, -578,  -64,  139, -449, -384,  193, -198, -538,  110,  -44, -577,
      -6,   66, -552, -164,  161, -460, -344,  205, -281, -504,  134,  -28,
    -576, -118,  156, -437, -381,  200, -220, -540,   83,   11, -568, -160,
     172, -414, -408,  188, -125, -572,  -32,  139, -492, -321,  205, -262,
    -504,  142,  -83, -574,    0,   48, -571, -106,  137, -501, -266,  190,
    -391, -406,  194, -186, -553,   83,  -13, -577,  -49,  103, -515, -280,
     201,  300, -506,  131,  -45, -578,  -80,  138, -462, -361,  201, -211,
    -554,   32,   74, -533, -235,  187, -372, -442,  182, -147, -566,   25,
      68, -535, -244,  194, -351, -463,  174, -125, -570,   15,   72, -550,
    -190,  172, -424, -385,  198, -218, -536,   96 } ;

  const size_t nlottery = 218 ;

  const double lottery[218] = {
    162, 671, 933, 414, 788, 730, 817, 33, 536, 875, 670, 236, 473, 167,
    877, 980, 316, 950, 456, 92, 517, 557, 956, 954, 104, 178, 794, 278,
    147, 773, 437, 435, 502, 610, 582, 780, 689, 562, 964, 791, 28, 97,
    848, 281, 858, 538, 660, 972, 671, 613, 867, 448, 738, 966, 139, 636,
    847, 659, 754, 243, 122, 455, 195, 968, 793, 59, 730, 361, 574, 522,
    97, 762, 431, 158, 429, 414, 22, 629, 788, 999, 187, 215, 810, 782,
    47, 34, 108, 986, 25, 644, 829, 630, 315, 567, 919, 331, 207, 412,
    242, 607, 668, 944, 749, 168, 864, 442, 533, 805, 372, 63, 458, 777,
    416, 340, 436, 140, 919, 350, 510, 572, 905, 900, 85, 389, 473, 758,
    444, 169, 625, 692, 140, 897, 672, 288, 312, 860, 724, 226, 884, 508,
    976, 741, 476, 417, 831, 15, 318, 432, 241, 114, 799, 955, 833, 358,
    935, 146, 630, 830, 440, 642, 356, 373, 271, 715, 367, 393, 190, 669,
    8, 861, 108, 795, 269, 590, 326, 866, 64, 523, 862, 840, 219, 382,
    998, 4, 628, 305, 747, 247, 34, 747, 729, 645, 856, 974, 24, 568, 24,
    694, 608, 480, 410, 729, 947, 293, 53, 930, 223, 203, 677, 227, 62,
    455, 387, 318, 562, 242, 428, 968 } ;

  const size_t nmavro = 50 ;

  const double mavro[50] = {
    2.00180, 2.00170, 2.00180, 2.00190, 2.00180, 2.00170, 2.00150,
    2.00140, 2.00150, 2.00150, 2.00170, 2.00180, 2.00180, 2.00190,
    2.00190, 2.00210, 2.00200, 2.00160, 2.00140, 2.00130, 2.00130,
    2.00150, 2.00150, 2.00160, 2.00150, 2.00140, 2.00130, 2.00140,
    2.00150, 2.00140, 2.00150, 2.00160, 2.00150, 2.00160, 2.00190,
    2.00200, 2.00200, 2.00210, 2.00220, 2.00230, 2.00240, 2.00250,
    2.00270, 2.00260, 2.00260, 2.00260, 2.00270, 2.00260, 2.00250,
    2.00240 } ;

  const size_t nmichelson = 100 ;

  const double michelson [100] = {
    299.85, 299.74, 299.90, 300.07, 299.93, 299.85, 299.95, 299.98,
    299.98, 299.88, 300.00, 299.98, 299.93, 299.65, 299.76, 299.81,
    300.00, 300.00, 299.96, 299.96, 299.96, 299.94, 299.96, 299.94,
    299.88, 299.80, 299.85, 299.88, 299.90, 299.84, 299.83, 299.79,
    299.81, 299.88, 299.88, 299.83, 299.80, 299.79, 299.76, 299.80,
    299.88, 299.88, 299.88, 299.86, 299.72, 299.72, 299.62, 299.86,
    299.97, 299.95, 299.88, 299.91, 299.85, 299.87, 299.84, 299.84,
    299.85, 299.84, 299.84, 299.84, 299.89, 299.81, 299.81, 299.82,
    299.80, 299.77, 299.76, 299.74, 299.75, 299.76, 299.91, 299.92,
    299.89, 299.86, 299.88, 299.72, 299.84, 299.85, 299.85, 299.78,
    299.89, 299.84, 299.78, 299.81, 299.76, 299.81, 299.79, 299.81,
    299.82, 299.85, 299.87, 299.87, 299.81, 299.74, 299.81, 299.94,
    299.95, 299.80, 299.81, 299.87 } ;

  const size_t npidigits = 5000 ;
  const double pidigits [5000] = { 3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8,
    9, 7, 9, 3, 2, 3, 8, 4, 6, 2, 6, 4, 3, 3, 8, 3, 2, 7, 9, 5, 0, 2,
    8, 8, 4, 1, 9, 7, 1, 6, 9, 3, 9, 9, 3, 7, 5, 1, 0, 5, 8, 2, 0, 9,
    7, 4, 9, 4, 4, 5, 9, 2, 3, 0, 7, 8, 1, 6, 4, 0, 6, 2, 8, 6, 2, 0,
    8, 9, 9, 8, 6, 2, 8, 0, 3, 4, 8, 2, 5, 3, 4, 2, 1, 1, 7, 0, 6, 7,
    9, 8, 2, 1, 4, 8, 0, 8, 6, 5, 1, 3, 2, 8, 2, 3, 0, 6, 6, 4, 7, 0,
    9, 3, 8, 4, 4, 6, 0, 9, 5, 5, 0, 5, 8, 2, 2, 3, 1, 7, 2, 5, 3, 5,
    9, 4, 0, 8, 1, 2, 8, 4, 8, 1, 1, 1, 7, 4, 5, 0, 2, 8, 4, 1, 0, 2,
    7, 0, 1, 9, 3, 8, 5, 2, 1, 1, 0, 5, 5, 5, 9, 6, 4, 4, 6, 2, 2, 9,
    4, 8, 9, 5, 4, 9, 3, 0, 3, 8, 1, 9, 6, 4, 4, 2, 8, 8, 1, 0, 9, 7,
    5, 6, 6, 5, 9, 3, 3, 4, 4, 6, 1, 2, 8, 4, 7, 5, 6, 4, 8, 2, 3, 3,
    7, 8, 6, 7, 8, 3, 1, 6, 5, 2, 7, 1, 2, 0, 1, 9, 0, 9, 1, 4, 5, 6,
    4, 8, 5, 6, 6, 9, 2, 3, 4, 6, 0, 3, 4, 8, 6, 1, 0, 4, 5, 4, 3, 2,
    6, 6, 4, 8, 2, 1, 3, 3, 9, 3, 6, 0, 7, 2, 6, 0, 2, 4, 9, 1, 4, 1,
    2, 7, 3, 7, 2, 4, 5, 8, 7, 0, 0, 6, 6, 0, 6, 3, 1, 5, 5, 8, 8, 1,
    7, 4, 8, 8, 1, 5, 2, 0, 9, 2, 0, 9, 6, 2, 8, 2, 9, 2, 5, 4, 0, 9,
    1, 7, 1, 5, 3, 6, 4, 3, 6, 7, 8, 9, 2, 5, 9, 0, 3, 6, 0, 0, 1, 1,
    3, 3, 0, 5, 3, 0, 5, 4, 8, 8, 2, 0, 4, 6, 6, 5, 2, 1, 3, 8, 4, 1,
    4, 6, 9, 5, 1, 9, 4, 1, 5, 1, 1, 6, 0, 9, 4, 3, 3, 0, 5, 7, 2, 7,
    0, 3, 6, 5, 7, 5, 9, 5, 9, 1, 9, 5, 3, 0, 9, 2, 1, 8, 6, 1, 1, 7,
    3, 8, 1, 9, 3, 2, 6, 1, 1, 7, 9, 3, 1, 0, 5, 1, 1, 8, 5, 4, 8, 0,
    7, 4, 4, 6, 2, 3, 7, 9, 9, 6, 2, 7, 4, 9, 5, 6, 7, 3, 5, 1, 8, 8,
    5, 7, 5, 2, 7, 2, 4, 8, 9, 1, 2, 2, 7, 9, 3, 8, 1, 8, 3, 0, 1, 1,
    9, 4, 9, 1, 2, 9, 8, 3, 3, 6, 7, 3, 3, 6, 2, 4, 4, 0, 6, 5, 6, 6,
    4, 3, 0, 8, 6, 0, 2, 1, 3, 9, 4, 9, 4, 6, 3, 9, 5, 2, 2, 4, 7, 3,
    7, 1, 9, 0, 7, 0, 2, 1, 7, 9, 8, 6, 0, 9, 4, 3, 7, 0, 2, 7, 7, 0,
    5, 3, 9, 2, 1, 7, 1, 7, 6, 2, 9, 3, 1, 7, 6, 7, 5, 2, 3, 8, 4, 6,
    7, 4, 8, 1, 8, 4, 6, 7, 6, 6, 9, 4, 0, 5, 1, 3, 2, 0, 0, 0, 5, 6,
    8, 1, 2, 7, 1, 4, 5, 2, 6, 3, 5, 6, 0, 8, 2, 7, 7, 8, 5, 7, 7, 1,
    3, 4, 2, 7, 5, 7, 7, 8, 9, 6, 0, 9, 1, 7, 3, 6, 3, 7, 1, 7, 8, 7,
    2, 1, 4, 6, 8, 4, 4, 0, 9, 0, 1, 2, 2, 4, 9, 5, 3, 4, 3, 0, 1, 4,
    6, 5, 4, 9, 5, 8, 5, 3, 7, 1, 0, 5, 0, 7, 9, 2, 2, 7, 9, 6, 8, 9,
    2, 5, 8, 9, 2, 3, 5, 4, 2, 0, 1, 9, 9, 5, 6, 1, 1, 2, 1, 2, 9, 0,
    2, 1, 9, 6, 0, 8, 6, 4, 0, 3, 4, 4, 1, 8, 1, 5, 9, 8, 1, 3, 6, 2,
    9, 7, 7, 4, 7, 7, 1, 3, 0, 9, 9, 6, 0, 5, 1, 8, 7, 0, 7, 2, 1, 1,
    3, 4, 9, 9, 9, 9, 9, 9, 8, 3, 7, 2, 9, 7, 8, 0, 4, 9, 9, 5, 1, 0,
    5, 9, 7, 3, 1, 7, 3, 2, 8, 1, 6, 0, 9, 6, 3, 1, 8, 5, 9, 5, 0, 2,
    4, 4, 5, 9, 4, 5, 5, 3, 4, 6, 9, 0, 8, 3, 0, 2, 6, 4, 2, 5, 2, 2,
    3, 0, 8, 2, 5, 3, 3, 4, 4, 6, 8, 5, 0, 3, 5, 2, 6, 1, 9, 3, 1, 1,
    8, 8, 1, 7, 1, 0, 1, 0, 0, 0, 3, 1, 3, 7, 8, 3, 8, 7, 5, 2, 8, 8,
    6, 5, 8, 7, 5, 3, 3, 2, 0, 8, 3, 8, 1, 4, 2, 0, 6, 1, 7, 1, 7, 7,
    6, 6, 9, 1, 4, 7, 3, 0, 3, 5, 9, 8, 2, 5, 3, 4, 9, 0, 4, 2, 8, 7,
    5, 5, 4, 6, 8, 7, 3, 1, 1, 5, 9, 5, 6, 2, 8, 6, 3, 8, 8, 2, 3, 5,
    3, 7, 8, 7, 5, 9, 3, 7, 5, 1, 9, 5, 7, 7, 8, 1, 8, 5, 7, 7, 3, 0,
    5, 3, 2, 1, 7, 1, 2, 2, 6, 8, 0, 6, 6, 1, 3, 0, 0, 1, 9, 2, 7, 8,
    7, 6, 6, 1, 1, 1, 9, 5, 9, 0, 9, 2, 1, 6, 4, 2, 0, 1, 9, 8, 9, 3,
    8, 0, 9, 5, 2, 5, 7, 2, 0, 1, 0, 6, 5, 4, 8, 5, 8, 6, 3, 2, 7, 8,
    8, 6, 5, 9, 3, 6, 1, 5, 3, 3, 8, 1, 8, 2, 7, 9, 6, 8, 2, 3, 0, 3,
    0, 1, 9, 5, 2, 0, 3, 5, 3, 0, 1, 8, 5, 2, 9, 6, 8, 9, 9, 5, 7, 7,
    3, 6, 2, 2, 5, 9, 9, 4, 1, 3, 8, 9, 1, 2, 4, 9, 7, 2, 1, 7, 7, 5,
    2, 8, 3, 4, 7, 9, 1, 3, 1, 5, 1, 5, 5, 7, 4, 8, 5, 7, 2, 4, 2, 4,
    5, 4, 1, 5, 0, 6, 9, 5, 9, 5, 0, 8, 2, 9, 5, 3, 3, 1, 1, 6, 8, 6,
    1, 7, 2, 7, 8, 5, 5, 8, 8, 9, 0, 7, 5, 0, 9, 8, 3, 8, 1, 7, 5, 4,
    6, 3, 7, 4, 6, 4, 9, 3, 9, 3, 1, 9, 2, 5, 5, 0, 6, 0, 4, 0, 0, 9,
    2, 7, 7, 0, 1, 6, 7, 1, 1, 3, 9, 0, 0, 9, 8, 4, 8, 8, 2, 4, 0, 1,
    2, 8, 5, 8, 3, 6, 1, 6, 0, 3, 5, 6, 3, 7, 0, 7, 6, 6, 0, 1, 0, 4,
    7, 1, 0, 1, 8, 1, 9, 4, 2, 9, 5, 5, 5, 9, 6, 1, 9, 8, 9, 4, 6, 7,
    6, 7, 8, 3, 7, 4, 4, 9, 4, 4, 8, 2, 5, 5, 3, 7, 9, 7, 7, 4, 7, 2,
    6, 8, 4, 7, 1, 0, 4, 0, 4, 7, 5, 3, 4, 6, 4, 6, 2, 0, 8, 0, 4, 6,
    6, 8, 4, 2, 5, 9, 0, 6, 9, 4, 9, 1, 2, 9, 3, 3, 1, 3, 6, 7, 7, 0,
    2, 8, 9, 8, 9, 1, 5, 2, 1, 0, 4, 7, 5, 2, 1, 6, 2, 0, 5, 6, 9, 6,
    6, 0, 2, 4, 0, 5, 8, 0, 3, 8, 1, 5, 0, 1, 9, 3, 5, 1, 1, 2, 5, 3,
    3, 8, 2, 4, 3, 0, 0, 3, 5, 5, 8, 7, 6, 4, 0, 2, 4, 7, 4, 9, 6, 4,
    7, 3, 2, 6, 3, 9, 1, 4, 1, 9, 9, 2, 7, 2, 6, 0, 4, 2, 6, 9, 9, 2,
    2, 7, 9, 6, 7, 8, 2, 3, 5, 4, 7, 8, 1, 6, 3, 6, 0, 0, 9, 3, 4, 1,
    7, 2, 1, 6, 4, 1, 2, 1, 9, 9, 2, 4, 5, 8, 6, 3, 1, 5, 0, 3, 0, 2,
    8, 6, 1, 8, 2, 9, 7, 4, 5, 5, 5, 7, 0, 6, 7, 4, 9, 8, 3, 8, 5, 0,
    5, 4, 9, 4, 5, 8, 8, 5, 8, 6, 9, 2, 6, 9, 9, 5, 6, 9, 0, 9, 2, 7,
    2, 1, 0, 7, 9, 7, 5, 0, 9, 3, 0, 2, 9, 5, 5, 3, 2, 1, 1, 6, 5, 3,
    4, 4, 9, 8, 7, 2, 0, 2, 7, 5, 5, 9, 6, 0, 2, 3, 6, 4, 8, 0, 6, 6,
    5, 4, 9, 9, 1, 1, 9, 8, 8, 1, 8, 3, 4, 7, 9, 7, 7, 5, 3, 5, 6, 6,
    3, 6, 9, 8, 0, 7, 4, 2, 6, 5, 4, 2, 5, 2, 7, 8, 6, 2, 5, 5, 1, 8,
    1, 8, 4, 1, 7, 5, 7, 4, 6, 7, 2, 8, 9, 0, 9, 7, 7, 7, 7, 2, 7, 9,
    3, 8, 0, 0, 0, 8, 1, 6, 4, 7, 0, 6, 0, 0, 1, 6, 1, 4, 5, 2, 4, 9,
    1, 9, 2, 1, 7, 3, 2, 1, 7, 2, 1, 4, 7, 7, 2, 3, 5, 0, 1, 4, 1, 4,
    4, 1, 9, 7, 3, 5, 6, 8, 5, 4, 8, 1, 6, 1, 3, 6, 1, 1, 5, 7, 3, 5,
    2, 5, 5, 2, 1, 3, 3, 4, 7, 5, 7, 4, 1, 8, 4, 9, 4, 6, 8, 4, 3, 8,
    5, 2, 3, 3, 2, 3, 9, 0, 7, 3, 9, 4, 1, 4, 3, 3, 3, 4, 5, 4, 7, 7,
    6, 2, 4, 1, 6, 8, 6, 2, 5, 1, 8, 9, 8, 3, 5, 6, 9, 4, 8, 5, 5, 6,
    2, 0, 9, 9, 2, 1, 9, 2, 2, 2, 1, 8, 4, 2, 7, 2, 5, 5, 0, 2, 5, 4,
    2, 5, 6, 8, 8, 7, 6, 7, 1, 7, 9, 0, 4, 9, 4, 6, 0, 1, 6, 5, 3, 4,
    6, 6, 8, 0, 4, 9, 8, 8, 6, 2, 7, 2, 3, 2, 7, 9, 1, 7, 8, 6, 0, 8,
    5, 7, 8, 4, 3, 8, 3, 8, 2, 7, 9, 6, 7, 9, 7, 6, 6, 8, 1, 4, 5, 4,
    1, 0, 0, 9, 5, 3, 8, 8, 3, 7, 8, 6, 3, 6, 0, 9, 5, 0, 6, 8, 0, 0,
    6, 4, 2, 2, 5, 1, 2, 5, 2, 0, 5, 1, 1, 7, 3, 9, 2, 9, 8, 4, 8, 9,
    6, 0, 8, 4, 1, 2, 8, 4, 8, 8, 6, 2, 6, 9, 4, 5, 6, 0, 4, 2, 4, 1,
    9, 6, 5, 2, 8, 5, 0, 2, 2, 2, 1, 0, 6, 6, 1, 1, 8, 6, 3, 0, 6, 7,
    4, 4, 2, 7, 8, 6, 2, 2, 0, 3, 9, 1, 9, 4, 9, 4, 5, 0, 4, 7, 1, 2,
    3, 7, 1, 3, 7, 8, 6, 9, 6, 0, 9, 5, 6, 3, 6, 4, 3, 7, 1, 9, 1, 7,
    2, 8, 7, 4, 6, 7, 7, 6, 4, 6, 5, 7, 5, 7, 3, 9, 6, 2, 4, 1, 3, 8,
    9, 0, 8, 6, 5, 8, 3, 2, 6, 4, 5, 9, 9, 5, 8, 1, 3, 3, 9, 0, 4, 7,
    8, 0, 2, 7, 5, 9, 0, 0, 9, 9, 4, 6, 5, 7, 6, 4, 0, 7, 8, 9, 5, 1,
    2, 6, 9, 4, 6, 8, 3, 9, 8, 3, 5, 2, 5, 9, 5, 7, 0, 9, 8, 2, 5, 8,
    2, 2, 6, 2, 0, 5, 2, 2, 4, 8, 9, 4, 0, 7, 7, 2, 6, 7, 1, 9, 4, 7,
    8, 2, 6, 8, 4, 8, 2, 6, 0, 1, 4, 7, 6, 9, 9, 0, 9, 0, 2, 6, 4, 0,
    1, 3, 6, 3, 9, 4, 4, 3, 7, 4, 5, 5, 3, 0, 5, 0, 6, 8, 2, 0, 3, 4,
    9, 6, 2, 5, 2, 4, 5, 1, 7, 4, 9, 3, 9, 9, 6, 5, 1, 4, 3, 1, 4, 2,
    9, 8, 0, 9, 1, 9, 0, 6, 5, 9, 2, 5, 0, 9, 3, 7, 2, 2, 1, 6, 9, 6,
    4, 6, 1, 5, 1, 5, 7, 0, 9, 8, 5, 8, 3, 8, 7, 4, 1, 0, 5, 9, 7, 8,
    8, 5, 9, 5, 9, 7, 7, 2, 9, 7, 5, 4, 9, 8, 9, 3, 0, 1, 6, 1, 7, 5,
    3, 9, 2, 8, 4, 6, 8, 1, 3, 8, 2, 6, 8, 6, 8, 3, 8, 6, 8, 9, 4, 2,
    7, 7, 4, 1, 5, 5, 9, 9, 1, 8, 5, 5, 9, 2, 5, 2, 4, 5, 9, 5, 3, 9,
    5, 9, 4, 3, 1, 0, 4, 9, 9, 7, 2, 5, 2, 4, 6, 8, 0, 8, 4, 5, 9, 8,
    7, 2, 7, 3, 6, 4, 4, 6, 9, 5, 8, 4, 8, 6, 5, 3, 8, 3, 6, 7, 3, 6,
    2, 2, 2, 6, 2, 6, 0, 9, 9, 1, 2, 4, 6, 0, 8, 0, 5, 1, 2, 4, 3, 8,
    8, 4, 3, 9, 0, 4, 5, 1, 2, 4, 4, 1, 3, 6, 5, 4, 9, 7, 6, 2, 7, 8,
    0, 7, 9, 7, 7, 1, 5, 6, 9, 1, 4, 3, 5, 9, 9, 7, 7, 0, 0, 1, 2, 9,
    6, 1, 6, 0, 8, 9, 4, 4, 1, 6, 9, 4, 8, 6, 8, 5, 5, 5, 8, 4, 8, 4,
    0, 6, 3, 5, 3, 4, 2, 2, 0, 7, 2, 2, 2, 5, 8, 2, 8, 4, 8, 8, 6, 4,
    8, 1, 5, 8, 4, 5, 6, 0, 2, 8, 5, 0, 6, 0, 1, 6, 8, 4, 2, 7, 3, 9,
    4, 5, 2, 2, 6, 7, 4, 6, 7, 6, 7, 8, 8, 9, 5, 2, 5, 2, 1, 3, 8, 5,
    2, 2, 5, 4, 9, 9, 5, 4, 6, 6, 6, 7, 2, 7, 8, 2, 3, 9, 8, 6, 4, 5,
    6, 5, 9, 6, 1, 1, 6, 3, 5, 4, 8, 8, 6, 2, 3, 0, 5, 7, 7, 4, 5, 6,
    4, 9, 8, 0, 3, 5, 5, 9, 3, 6, 3, 4, 5, 6, 8, 1, 7, 4, 3, 2, 4, 1,
    1, 2, 5, 1, 5, 0, 7, 6, 0, 6, 9, 4, 7, 9, 4, 5, 1, 0, 9, 6, 5, 9,
    6, 0, 9, 4, 0, 2, 5, 2, 2, 8, 8, 7, 9, 7, 1, 0, 8, 9, 3, 1, 4, 5,
    6, 6, 9, 1, 3, 6, 8, 6, 7, 2, 2, 8, 7, 4, 8, 9, 4, 0, 5, 6, 0, 1,
    0, 1, 5, 0, 3, 3, 0, 8, 6, 1, 7, 9, 2, 8, 6, 8, 0, 9, 2, 0, 8, 7,
    4, 7, 6, 0, 9, 1, 7, 8, 2, 4, 9, 3, 8, 5, 8, 9, 0, 0, 9, 7, 1, 4,
    9, 0, 9, 6, 7, 5, 9, 8, 5, 2, 6, 1, 3, 6, 5, 5, 4, 9, 7, 8, 1, 8,
    9, 3, 1, 2, 9, 7, 8, 4, 8, 2, 1, 6, 8, 2, 9, 9, 8, 9, 4, 8, 7, 2,
    2, 6, 5, 8, 8, 0, 4, 8, 5, 7, 5, 6, 4, 0, 1, 4, 2, 7, 0, 4, 7, 7,
    5, 5, 5, 1, 3, 2, 3, 7, 9, 6, 4, 1, 4, 5, 1, 5, 2, 3, 7, 4, 6, 2,
    3, 4, 3, 6, 4, 5, 4, 2, 8, 5, 8, 4, 4, 4, 7, 9, 5, 2, 6, 5, 8, 6,
    7, 8, 2, 1, 0, 5, 1, 1, 4, 1, 3, 5, 4, 7, 3, 5, 7, 3, 9, 5, 2, 3,
    1, 1, 3, 4, 2, 7, 1, 6, 6, 1, 0, 2, 1, 3, 5, 9, 6, 9, 5, 3, 6, 2,
    3, 1, 4, 4, 2, 9, 5, 2, 4, 8, 4, 9, 3, 7, 1, 8, 7, 1, 1, 0, 1, 4,
    5, 7, 6, 5, 4, 0, 3, 5, 9, 0, 2, 7, 9, 9, 3, 4, 4, 0, 3, 7, 4, 2,
    0, 0, 7, 3, 1, 0, 5, 7, 8, 5, 3, 9, 0, 6, 2, 1, 9, 8, 3, 8, 7, 4,
    4, 7, 8, 0, 8, 4, 7, 8, 4, 8, 9, 6, 8, 3, 3, 2, 1, 4, 4, 5, 7, 1,
    3, 8, 6, 8, 7, 5, 1, 9, 4, 3, 5, 0, 6, 4, 3, 0, 2, 1, 8, 4, 5, 3,
    1, 9, 1, 0, 4, 8, 4, 8, 1, 0, 0, 5, 3, 7, 0, 6, 1, 4, 6, 8, 0, 6,
    7, 4, 9, 1, 9, 2, 7, 8, 1, 9, 1, 1, 9, 7, 9, 3, 9, 9, 5, 2, 0, 6,
    1, 4, 1, 9, 6, 6, 3, 4, 2, 8, 7, 5, 4, 4, 4, 0, 6, 4, 3, 7, 4, 5,
    1, 2, 3, 7, 1, 8, 1, 9, 2, 1, 7, 9, 9, 9, 8, 3, 9, 1, 0, 1, 5, 9,
    1, 9, 5, 6, 1, 8, 1, 4, 6, 7, 5, 1, 4, 2, 6, 9, 1, 2, 3, 9, 7, 4,
    8, 9, 4, 0, 9, 0, 7, 1, 8, 6, 4, 9, 4, 2, 3, 1, 9, 6, 1, 5, 6, 7,
    9, 4, 5, 2, 0, 8, 0, 9, 5, 1, 4, 6, 5, 5, 0, 2, 2, 5, 2, 3, 1, 6,
    0, 3, 8, 8, 1, 9, 3, 0, 1, 4, 2, 0, 9, 3, 7, 6, 2, 1, 3, 7, 8, 5,
    5, 9, 5, 6, 6, 3, 8, 9, 3, 7, 7, 8, 7, 0, 8, 3, 0, 3, 9, 0, 6, 9,
    7, 9, 2, 0, 7, 7, 3, 4, 6, 7, 2, 2, 1, 8, 2, 5, 6, 2, 5, 9, 9, 6,
    6, 1, 5, 0, 1, 4, 2, 1, 5, 0, 3, 0, 6, 8, 0, 3, 8, 4, 4, 7, 7, 3,
    4, 5, 4, 9, 2, 0, 2, 6, 0, 5, 4, 1, 4, 6, 6, 5, 9, 2, 5, 2, 0, 1,
    4, 9, 7, 4, 4, 2, 8, 5, 0, 7, 3, 2, 5, 1, 8, 6, 6, 6, 0, 0, 2, 1,
    3, 2, 4, 3, 4, 0, 8, 8, 1, 9, 0, 7, 1, 0, 4, 8, 6, 3, 3, 1, 7, 3,
    4, 6, 4, 9, 6, 5, 1, 4, 5, 3, 9, 0, 5, 7, 9, 6, 2, 6, 8, 5, 6, 1,
    0, 0, 5, 5, 0, 8, 1, 0, 6, 6, 5, 8, 7, 9, 6, 9, 9, 8, 1, 6, 3, 5,
    7, 4, 7, 3, 6, 3, 8, 4, 0, 5, 2, 5, 7, 1, 4, 5, 9, 1, 0, 2, 8, 9,
    7, 0, 6, 4, 1, 4, 0, 1, 1, 0, 9, 7, 1, 2, 0, 6, 2, 8, 0, 4, 3, 9,
    0, 3, 9, 7, 5, 9, 5, 1, 5, 6, 7, 7, 1, 5, 7, 7, 0, 0, 4, 2, 0, 3,
    3, 7, 8, 6, 9, 9, 3, 6, 0, 0, 7, 2, 3, 0, 5, 5, 8, 7, 6, 3, 1, 7,
    6, 3, 5, 9, 4, 2, 1, 8, 7, 3, 1, 2, 5, 1, 4, 7, 1, 2, 0, 5, 3, 2,
    9, 2, 8, 1, 9, 1, 8, 2, 6, 1, 8, 6, 1, 2, 5, 8, 6, 7, 3, 2, 1, 5,
    7, 9, 1, 9, 8, 4, 1, 4, 8, 4, 8, 8, 2, 9, 1, 6, 4, 4, 7, 0, 6, 0,
    9, 5, 7, 5, 2, 7, 0, 6, 9, 5, 7, 2, 2, 0, 9, 1, 7, 5, 6, 7, 1, 1,
    6, 7, 2, 2, 9, 1, 0, 9, 8, 1, 6, 9, 0, 9, 1, 5, 2, 8, 0, 1, 7, 3,
    5, 0, 6, 7, 1, 2, 7, 4, 8, 5, 8, 3, 2, 2, 2, 8, 7, 1, 8, 3, 5, 2,
    0, 9, 3, 5, 3, 9, 6, 5, 7, 2, 5, 1, 2, 1, 0, 8, 3, 5, 7, 9, 1, 5,
    1, 3, 6, 9, 8, 8, 2, 0, 9, 1, 4, 4, 4, 2, 1, 0, 0, 6, 7, 5, 1, 0,
    3, 3, 4, 6, 7, 1, 1, 0, 3, 1, 4, 1, 2, 6, 7, 1, 1, 1, 3, 6, 9, 9,
    0, 8, 6, 5, 8, 5, 1, 6, 3, 9, 8, 3, 1, 5, 0, 1, 9, 7, 0, 1, 6, 5,
    1, 5, 1, 1, 6, 8, 5, 1, 7, 1, 4, 3, 7, 6, 5, 7, 6, 1, 8, 3, 5, 1,
    5, 5, 6, 5, 0, 8, 8, 4, 9, 0, 9, 9, 8, 9, 8, 5, 9, 9, 8, 2, 3, 8,
    7, 3, 4, 5, 5, 2, 8, 3, 3, 1, 6, 3, 5, 5, 0, 7, 6, 4, 7, 9, 1, 8,
    5, 3, 5, 8, 9, 3, 2, 2, 6, 1, 8, 5, 4, 8, 9, 6, 3, 2, 1, 3, 2, 9,
    3, 3, 0, 8, 9, 8, 5, 7, 0, 6, 4, 2, 0, 4, 6, 7, 5, 2, 5, 9, 0, 7,
    0, 9, 1, 5, 4, 8, 1, 4, 1, 6, 5, 4, 9, 8, 5, 9, 4, 6, 1, 6, 3, 7,
    1, 8, 0, 2, 7, 0, 9, 8, 1, 9, 9, 4, 3, 0, 9, 9, 2, 4, 4, 8, 8, 9,
    5, 7, 5, 7, 1, 2, 8, 2, 8, 9, 0, 5, 9, 2, 3, 2, 3, 3, 2, 6, 0, 9,
    7, 2, 9, 9, 7, 1, 2, 0, 8, 4, 4, 3, 3, 5, 7, 3, 2, 6, 5, 4, 8, 9,
    3, 8, 2, 3, 9, 1, 1, 9, 3, 2, 5, 9, 7, 4, 6, 3, 6, 6, 7, 3, 0, 5,
    8, 3, 6, 0, 4, 1, 4, 2, 8, 1, 3, 8, 8, 3, 0, 3, 2, 0, 3, 8, 2, 4,
    9, 0, 3, 7, 5, 8, 9, 8, 5, 2, 4, 3, 7, 4, 4, 1, 7, 0, 2, 9, 1, 3,
    2, 7, 6, 5, 6, 1, 8, 0, 9, 3, 7, 7, 3, 4, 4, 4, 0, 3, 0, 7, 0, 7,
    4, 6, 9, 2, 1, 1, 2, 0, 1, 9, 1, 3, 0, 2, 0, 3, 3, 0, 3, 8, 0, 1,
    9, 7, 6, 2, 1, 1, 0, 1, 1, 0, 0, 4, 4, 9, 2, 9, 3, 2, 1, 5, 1, 6,
    0, 8, 4, 2, 4, 4, 4, 8, 5, 9, 6, 3, 7, 6, 6, 9, 8, 3, 8, 9, 5, 2,
    2, 8, 6, 8, 4, 7, 8, 3, 1, 2, 3, 5, 5, 2, 6, 5, 8, 2, 1, 3, 1, 4,
    4, 9, 5, 7, 6, 8, 5, 7, 2, 6, 2, 4, 3, 3, 4, 4, 1, 8, 9, 3, 0, 3,
    9, 6, 8, 6, 4, 2, 6, 2, 4, 3, 4, 1, 0, 7, 7, 3, 2, 2, 6, 9, 7, 8,
    0, 2, 8, 0, 7, 3, 1, 8, 9, 1, 5, 4, 4, 1, 1, 0, 1, 0, 4, 4, 6, 8,
    2, 3, 2, 5, 2, 7, 1, 6, 2, 0, 1, 0, 5, 2, 6, 5, 2, 2, 7, 2, 1, 1,
    1, 6, 6, 0, 3, 9, 6, 6, 6, 5, 5, 7, 3, 0, 9, 2, 5, 4, 7, 1, 1, 0,
    5, 5, 7, 8, 5, 3, 7, 6, 3, 4, 6, 6, 8, 2, 0, 6, 5, 3, 1, 0, 9, 8,
    9, 6, 5, 2, 6, 9, 1, 8, 6, 2, 0, 5, 6, 4, 7, 6, 9, 3, 1, 2, 5, 7,
    0, 5, 8, 6, 3, 5, 6, 6, 2, 0, 1, 8, 5, 5, 8, 1, 0, 0, 7, 2, 9, 3,
    6, 0, 6, 5, 9, 8, 7, 6, 4, 8, 6, 1, 1, 7, 9, 1, 0, 4, 5, 3, 3, 4,
    8, 8, 5, 0, 3, 4, 6, 1, 1, 3, 6, 5, 7, 6, 8, 6, 7, 5, 3, 2, 4, 9,
    4, 4, 1, 6, 6, 8, 0, 3, 9, 6, 2, 6, 5, 7, 9, 7, 8, 7, 7, 1, 8, 5,
    5, 6, 0, 8, 4, 5, 5, 2, 9, 6, 5, 4, 1, 2, 6, 6, 5, 4, 0, 8, 5, 3,
    0, 6, 1, 4, 3, 4, 4, 4, 3, 1, 8, 5, 8, 6, 7, 6, 9, 7, 5, 1, 4, 5,
    6, 6, 1, 4, 0, 6, 8, 0, 0, 7, 0, 0, 2, 3, 7, 8, 7, 7, 6, 5, 9, 1,
    3, 4, 4, 0, 1, 7, 1, 2, 7, 4, 9, 4, 7, 0, 4, 2, 0, 5, 6, 2, 2, 3,
    0, 5, 3, 8, 9, 9, 4, 5, 6, 1, 3, 1, 4, 0, 7, 1, 1, 2, 7, 0, 0, 0,
    4, 0, 7, 8, 5, 4, 7, 3, 3, 2, 6, 9, 9, 3, 9, 0, 8, 1, 4, 5, 4, 6,
    6, 4, 6, 4, 5, 8, 8, 0, 7, 9, 7, 2, 7, 0, 8, 2, 6, 6, 8, 3, 0, 6,
    3, 4, 3, 2, 8, 5, 8, 7, 8, 5, 6, 9, 8, 3, 0, 5, 2, 3, 5, 8, 0, 8,
    9, 3, 3, 0, 6, 5, 7, 5, 7, 4, 0, 6, 7, 9, 5, 4, 5, 7, 1, 6, 3, 7,
    7, 5, 2, 5, 4, 2, 0, 2, 1, 1, 4, 9, 5, 5, 7, 6, 1, 5, 8, 1, 4, 0,
    0, 2, 5, 0, 1, 2, 6, 2, 2, 8, 5, 9, 4, 1, 3, 0, 2, 1, 6, 4, 7, 1,
    5, 5, 0, 9, 7, 9, 2, 5, 9, 2, 3, 0, 9, 9, 0, 7, 9, 6, 5, 4, 7, 3,
    7, 6, 1, 2, 5, 5, 1, 7, 6, 5, 6, 7, 5, 1, 3, 5, 7, 5, 1, 7, 8, 2,
    9, 6, 6, 6, 4, 5, 4, 7, 7, 9, 1, 7, 4, 5, 0, 1, 1, 2, 9, 9, 6, 1,
    4, 8, 9, 0, 3, 0, 4, 6, 3, 9, 9, 4, 7, 1, 3, 2, 9, 6, 2, 1, 0, 7,
    3, 4, 0, 4, 3, 7, 5, 1, 8, 9, 5, 7, 3, 5, 9, 6, 1, 4, 5, 8, 9, 0,
    1, 9, 3, 8, 9, 7, 1, 3, 1, 1, 1, 7, 9, 0, 4, 2, 9, 7, 8, 2, 8, 5,
    6, 4, 7, 5, 0, 3, 2, 0, 3, 1, 9, 8, 6, 9, 1, 5, 1, 4, 0, 2, 8, 7,
    0, 8, 0, 8, 5, 9, 9, 0, 4, 8, 0, 1, 0, 9, 4, 1, 2, 1, 4, 7, 2, 2,
    1, 3, 1, 7, 9, 4, 7, 6, 4, 7, 7, 7, 2, 6, 2, 2, 4, 1, 4, 2, 5, 4,
    8, 5, 4, 5, 4, 0, 3, 3, 2, 1, 5, 7, 1, 8, 5, 3, 0, 6, 1, 4, 2, 2,
    8, 8, 1, 3, 7, 5, 8, 5, 0, 4, 3, 0, 6, 3, 3, 2, 1, 7, 5, 1, 8, 2,
    9, 7, 9, 8, 6, 6, 2, 2, 3, 7, 1, 7, 2, 1, 5, 9, 1, 6, 0, 7, 7, 1,
    6, 6, 9, 2, 5, 4, 7, 4, 8, 7, 3, 8, 9, 8, 6, 6, 5, 4, 9, 4, 9, 4,
    5, 0, 1, 1, 4, 6, 5, 4, 0, 6, 2, 8, 4, 3, 3, 6, 6, 3, 9, 3, 7, 9,
    0, 0, 3, 9, 7, 6, 9, 2, 6, 5, 6, 7, 2, 1, 4, 6, 3, 8, 5, 3, 0, 6,
    7, 3, 6, 0, 9, 6, 5, 7, 1, 2, 0, 9, 1, 8, 0, 7, 6, 3, 8, 3, 2, 7,
    1, 6, 6, 4, 1, 6, 2, 7, 4, 8, 8, 8, 8, 0, 0, 7, 8, 6, 9, 2, 5, 6,
    0, 2, 9, 0, 2, 2, 8, 4, 7, 2, 1, 0, 4, 0, 3, 1, 7, 2, 1, 1, 8, 6,
    0, 8, 2, 0, 4, 1, 9, 0, 0, 0, 4, 2, 2, 9, 6, 6, 1, 7, 1, 1, 9, 6,
    3, 7, 7, 9, 2, 1, 3, 3, 7, 5, 7, 5, 1, 1, 4, 9, 5, 9, 5, 0, 1, 5,
    6, 6, 0, 4, 9, 6, 3, 1, 8, 6, 2, 9, 4, 7, 2, 6, 5, 4, 7, 3, 6, 4,
    2, 5, 2, 3, 0, 8, 1, 7, 7, 0, 3, 6, 7, 5, 1, 5, 9, 0, 6, 7, 3, 5,
    0, 2, 3, 5, 0, 7, 2, 8, 3, 5, 4, 0, 5, 6, 7, 0, 4, 0, 3, 8, 6, 7,
    4, 3, 5, 1, 3, 6, 2, 2, 2, 2, 4, 7, 7, 1, 5, 8, 9, 1, 5, 0, 4, 9,
    5, 3, 0, 9, 8, 4, 4, 4, 8, 9, 3, 3, 3, 0, 9, 6, 3, 4, 0, 8, 7, 8,
    0, 7, 6, 9, 3, 2, 5, 9, 9, 3, 9, 7, 8, 0, 5, 4, 1, 9, 3, 4, 1, 4,
    4, 7, 3, 7, 7, 4, 4, 1, 8, 4, 2, 6, 3, 1, 2, 9, 8, 6, 0, 8, 0, 9,
    9, 8, 8, 8, 6, 8, 7, 4, 1, 3, 2, 6, 0, 4, 7, 2 } ;

  const size_t nacc1 = 3 ;
  const double numacc1[3] = { 10000001, 10000003, 10000002 } ;

  const size_t nacc2 = 1001 ;
  double numacc2[1001] ;

  const size_t nacc3 = 1001 ;
  double numacc3[1001] ;

  const size_t nacc4 = 1001 ;
  double numacc4[1001] ;

  numacc2[0] = 1.2 ;
  numacc3[0] = 1000000.2 ;
  numacc4[0] = 10000000.2 ;

  for (i = 1 ; i < 1000  ; i += 2)
    {
      numacc2[i] = 1.1 ;
      numacc2[i+1] = 1.3 ;
      numacc3[i] = 1000000.1 ;
      numacc3[i+1] = 1000000.3 ;
      numacc4[i] = 10000000.1 ;
      numacc4[i+1] = 10000000.3 ;
    }

  gsl_ieee_env_setup ();

  {
    double mean = gsl_stats_mean (lew, 1, nlew);
    double sd = gsl_stats_sd (lew, 1, nlew);
    double lag1 = gsl_stats_lag1_autocorrelation (lew, 1, nlew);

    double expected_mean = -177.435000000000;
    double expected_sd = 277.332168044316;
    double expected_lag1 = -0.307304800605679;

    gsl_test_rel (mean, expected_mean, 1e-15, "lew gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-15, "lew gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-14, "lew autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (lottery, 1, nlottery);
    double sd = gsl_stats_sd (lottery, 1, nlottery);
    double lag1 = gsl_stats_lag1_autocorrelation (lottery, 1, nlottery);

    double expected_mean = 518.958715596330;
    double expected_sd = 291.699727470969;
    double expected_lag1 = -0.120948622967393;

    gsl_test_rel (mean, expected_mean, 1e-15, "lottery gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-15, "lottery gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-14, "lottery autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (mavro, 1, nmavro);
    double sd = gsl_stats_sd (mavro, 1, nmavro);
    double lag1 = gsl_stats_lag1_autocorrelation (mavro, 1, nmavro);

    double expected_mean = 2.00185600000000;
    double expected_sd = 0.000429123454003053;
    double expected_lag1 = 0.937989183438248;

    gsl_test_rel (mean, expected_mean, 1e-15, "mavro gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-13, "mavro gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-13, "mavro autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (michelson, 1, nmichelson);
    double sd = gsl_stats_sd (michelson, 1, nmichelson);
    double lag1 = gsl_stats_lag1_autocorrelation (michelson, 1, nmichelson);

    double expected_mean = 299.852400000000;
    double expected_sd = 0.0790105478190518;
    double expected_lag1 = 0.535199668621283;

    gsl_test_rel (mean, expected_mean, 1e-15, "michelson gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-13, "michelson gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-13, "michelson autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (pidigits, 1, npidigits);
    double sd = gsl_stats_sd (pidigits, 1, npidigits);
    double lag1 = gsl_stats_lag1_autocorrelation (pidigits, 1, npidigits);

    double expected_mean = 4.53480000000000;
    double expected_sd = 2.86733906028871;
    double expected_lag1 = -0.00355099287237972;

    gsl_test_rel (mean, expected_mean, 1e-14, "pidigits gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-15, "pidigits gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-14, "pidigits autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (numacc1, 1, nacc1);
    double sd = gsl_stats_sd (numacc1, 1, nacc1);
    double lag1 = gsl_stats_lag1_autocorrelation (numacc1, 1, nacc1);

    double expected_mean = 10000002;
    double expected_sd = 1;
    double expected_lag1 = -0.5;

    gsl_test_rel (mean, expected_mean, 1e-15, "acc1 gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-15, "acc1 gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-15, "acc1 autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (numacc2, 1, nacc2);
    double sd = gsl_stats_sd (numacc2, 1, nacc2);
    double lag1 = gsl_stats_lag1_autocorrelation (numacc2, 1, nacc2);

    double expected_mean = 1.2;
    double expected_sd = 0.1;
    double expected_lag1 = -0.999;

    gsl_test_rel (mean, expected_mean, 1e-15, "acc2 gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-15, "acc2 gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-10, "acc2 autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (numacc3, 1, nacc3);
    double sd = gsl_stats_sd (numacc3, 1, nacc3);
    double lag1 = gsl_stats_lag1_autocorrelation (numacc3, 1, nacc3);

    double expected_mean = 1000000.2;
    double expected_sd = 0.1;
    double expected_lag1 = -0.999;

    gsl_test_rel (mean, expected_mean, 1e-15, "acc3 gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-9, "acc3 gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-10, "acc3 autocorrelation") ;
  }

  {
    double mean = gsl_stats_mean (numacc4, 1, nacc4);
    double sd = gsl_stats_sd (numacc4, 1, nacc4);
    double lag1 = gsl_stats_lag1_autocorrelation (numacc4, 1, nacc4);

    double expected_mean = 10000000.2;
    double expected_sd = 0.1;
    double expected_lag1 = -0.999;

    gsl_test_rel (mean, expected_mean, 1e-15, "acc4 gsl_stats_mean") ;
    gsl_test_rel (sd, expected_sd, 1e-7, "acc4 gsl_stats_sd") ;
    gsl_test_rel (lag1, expected_lag1, 1e-10, "acc4 autocorrelation") ;
  }

  return 0;
}

ext: c

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
/*
Copyright (c) 2015, Patrick Weltevrede
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "psrsalsa.h"
#include <gsl/gsl_linalg.h>
int linalg_solve_matrix_eq_gauss_jordan(double *matrixa, double *matrixb, int n, int m, verbose_definition verbose)
{
  int *orig_pivot_row, *orig_pivot_column, *pivot_column_done;
  int pivotnr, pivot_col, pivot_row, rowi, coli, i;
  double *ptr1, *ptr2, tmpvalue;
  double largest_element_value;
  int currow;
  orig_pivot_row = malloc(n*sizeof(int));
  orig_pivot_column = malloc(n*sizeof(int));
  pivot_column_done = malloc(n*sizeof(int));
  if(orig_pivot_row == NULL || orig_pivot_column == NULL || pivot_column_done == NULL) {
    printerror(verbose.debug, "ERROR linalg_solve_matrix_eq_gauss_jordan: Memory allocation error.");
    return 1;
  }
  for(coli = 0; coli < n; coli++)
    pivot_column_done[coli] = 0;
  for(pivotnr = 0; pivotnr < n; pivotnr++) {
    largest_element_value = -1.0;
    for(rowi = 0; rowi < n; rowi++) {
      if(pivot_column_done[rowi] != 1) {
 for(coli = 0; coli < n; coli++) {
   if(pivot_column_done[coli] == 0) {
     tmpvalue = matrixa[rowi*n+coli];
     if(fabs(tmpvalue) >= largest_element_value) {
       largest_element_value = fabs(tmpvalue);
       pivot_row=rowi;
       pivot_col=coli;
     }
   }else if(pivot_column_done[coli] > 1) {
     printerror(verbose.debug, "ERROR linalg_solve_matrix_eq_gauss_jordan: The matrix equation is singular, no solution can be determined.");
     return 2;
   }
 }
      }
    }
    if(largest_element_value <= 0.0) {
      printerror(verbose.debug, "ERROR linalg_solve_matrix_eq_gauss_jordan: The matrix equation is singular, no solution can be determined.");
      return 2;
    }
    pivot_column_done[pivot_col] += 1;
    if(pivot_row != pivot_col) {
      ptr1 = matrixa + pivot_row*n;
      ptr2 = matrixa + pivot_col*n;
      for(i = 0; i < n; i++) {
 tmpvalue = *ptr1;
 *ptr1 = *ptr2;
 *ptr2 = tmpvalue;
      }
      ptr1 = matrixb + pivot_row*n;
      ptr2 = matrixb + pivot_col*n;
      for(i = 0; i < n; i++) {
 tmpvalue = *ptr1;
 *ptr1 = *ptr2;
 *ptr2 = tmpvalue;
      }
    }
    orig_pivot_row[pivotnr] = pivot_row;
    orig_pivot_column[pivotnr] = pivot_col;
    matrixa[pivot_col*n+pivot_col]=1.0;
    ptr1 = matrixa + pivot_col*n;
    for(i = 0; i < n; i++) {
      *ptr1 /= largest_element_value;
      ptr1 += 1;
    }
    ptr1 = matrixb + pivot_col*m;
    for(i = 0; i < m; i++) {
      *ptr1 /= largest_element_value;
      ptr1 += 1;
    }
    for(currow = 0; currow < n; currow++) {
      if(currow != pivot_col) {
 i = currow*n+pivot_col;
 tmpvalue = matrixa[i];
 matrixa[i] = 0.0;
 ptr1 = matrixa + n*currow;
 ptr2 = matrixa + n*pivot_col;
 for(i=1; i <= n; i++) {
   *ptr1 -= (*ptr2)*tmpvalue;
   ptr1++;
   ptr2++;
 }
 ptr1 = matrixb + currow*m;
 ptr2 = matrixb + pivot_col*m;
 for(i = 0; i < m; i++) {
   *ptr1 -= (*ptr2)*tmpvalue;
   ptr1++;
   ptr2++;
 }
      }
    }
  }
  for(coli = n-1; coli >= 0; coli--) {
    if(orig_pivot_row[coli] != orig_pivot_column[coli]) {
      for(rowi = 0; rowi <= n; rowi++) {
 tmpvalue = matrixa[rowi*n+orig_pivot_row[coli]];
 matrixa[rowi*n+orig_pivot_row[coli]] = matrixa[rowi*n+orig_pivot_column[coli]];
 matrixa[rowi*n+orig_pivot_column[coli]] = tmpvalue;
      }
    }
  }
  free(pivot_column_done);
  free(orig_pivot_row);
  free(orig_pivot_column);
  return 0;
}

max_stars_repo_name: David-McKenna/psrsalsa

max_issues_repo_name: David-McKenna/psrsalsa

max_forks_repo_name: David-McKenna/psrsalsa


text:
/* multilarge_nlinear/nielsen.c
 *
 * Copyright (C) 2016 Patrick Alken
 *
 * 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.
 */

/*
 * This module contains routines for updating the Levenberg-Marquardt
 * damping parameter on each iteration using Nielsen's method:
 *
 * [1] H. B. Nielsen, K. Madsen, Introduction to Optimization and
 *     Data Fitting, Informatics and Mathematical Modeling,
 *     Technical University of Denmark (DTU), 2010.
 *
 * 3 routines are needed to implement the update procedure:
 *
 * 1. init   - initialize parameter prior to iteration
 * 2. accept - update parameter after a step has been accepted
 * 3. reject - update parameter after a step has been rejected
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

#define LM_ONE_THIRD         (0.333333333333333)

static int nielsen_init(const gsl_matrix * JTJ, const gsl_vector * diag,
                        double * mu, long * nu);
static int nielsen_accept(const double rho, double * mu, long * nu);
static int nielsen_reject(double * mu, long * nu);

static int
nielsen_init(const gsl_matrix * JTJ, const gsl_vector * diag,
             double * mu, long * nu)
{
  const double mu0 = 1.0e-3;
  size_t i;
  double max = -1.0;

  *nu = 2;

  if (JTJ == NULL)
    {
      /* not using cholesky solver */
      *mu = mu0;
      return GSL_SUCCESS;
    }

  /* set mu = mu0 * max(diag(J~^T J~)), with J~ = J D^{-1} */

  for (i = 0; i < JTJ->size1; ++i)
    {
      double JTJi = gsl_matrix_get(JTJ, i, i);
      double di = gsl_vector_get(diag, i);
      double val = JTJi / (di * di);
      max = GSL_MAX(max, val);
    }

  *mu = mu0 * max;

  return GSL_SUCCESS;
}

static int
nielsen_accept(const double rho, double * mu, long * nu)
{
  double b;

  /* reset nu */
  *nu = 2;

  b = 2.0 * rho - 1.0;
  b = 1.0 - b*b*b;
  *mu *= GSL_MAX(LM_ONE_THIRD, b);

  return GSL_SUCCESS;
}

static int
nielsen_reject(double * mu, long * nu)
{
  *mu *= (double) *nu;

  /* nu := 2*nu */
  *nu <<= 1;

  return GSL_SUCCESS;
}


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef SIM_PHYSICS_DELTA_COEFF_AERO_DATABASE_H_
#define SIM_PHYSICS_DELTA_COEFF_AERO_DATABASE_H_

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#include <stdint.h>

#include <string>

#include "common/c_math/linalg.h"
#include "common/c_math/vec3.h"
#include "common/macros.h"
#include "sim/sim_params.h"
#include "sim/sim_types.h"
#include "system/labels.h"

class DeltaCoeffAeroDatabase {
  friend class AeroTest;

 public:
  explicit DeltaCoeffAeroDatabase(const std::string &filename);
  ~DeltaCoeffAeroDatabase();

  // Returns the nominal force-moment coefficient after adjusting for
  // delta coefficients.
  void AdjustForceMomentCoeff(double alpha, const Vec &flaps,
                              ForceMoment *force_moment_coeff) const;

 private:
  // Returns true if the database passes basic sign checks.
  bool IsValid() const;

  bool applicable_flaps_[kNumFlaps];

  double reynolds_number_;
  gsl_vector *alphads_, *deltads_;

  gsl_matrix *CX_[kNumFlaps];
  gsl_matrix *CY_[kNumFlaps];
  gsl_matrix *CZ_[kNumFlaps];
  gsl_matrix *CL_[kNumFlaps];
  gsl_matrix *CM_[kNumFlaps];
  gsl_matrix *CN_[kNumFlaps];

  DISALLOW_COPY_AND_ASSIGN(DeltaCoeffAeroDatabase);
};

#endif  // SIM_PHYSICS_DELTA_COEFF_AERO_DATABASE_H_


text:
//
//  Author  : github.com/luncliff (luncliff@gmail.com)
//  License : CC BY 4.0
//
//  Note
//      Test framework adapter
//
#pragma once

#include <gsl/gsl>

void _require_(bool expr);
void _require_(bool expr, gsl::czstring<> file, size_t line);
void _println_(gsl::czstring<> message);

void _fail_now_(gsl::czstring<> message, //
                gsl::czstring<> file = __FILE__, size_t line = __LINE__);

max_stars_repo_name: Farwaykorse/coroutine

max_issues_repo_name: Farwaykorse/coroutine

max_forks_repo_name: Farwaykorse/coroutine


text:
/* poly/zsolve_init.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_errno.h>

gsl_poly_complex_workspace *
gsl_poly_complex_workspace_alloc (size_t n)
{
  size_t nc ;

  gsl_poly_complex_workspace * w ;

  if (n == 0)
    {
      GSL_ERROR_VAL ("matrix size n must be positive integer", GSL_EDOM, 0);
    }

  w = (gsl_poly_complex_workspace *)
    malloc (sizeof(gsl_poly_complex_workspace));

  if (w == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for struct", GSL_ENOMEM, 0);
    }

  nc = n - 1;

  w->nc = nc;

  w->matrix = (double *) malloc (nc * nc * sizeof(double));

  if (w->matrix == 0)
    {
      free (w) ;       /* error in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for workspace matrix",
                        GSL_ENOMEM, 0);
    }

  return w ;
}

void
gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w)
{
  free(w->matrix) ;
  free(w);
}


text:
#pragma once
#include <Manager.h>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <limits>
#include <vector>
#include "../types/def.h"
#include <gsl\gsl>
#include "SFML/Window.hpp"
#include <utility>
#include <../_entity/Entity.h>
#include <../_entity/EntityManager.h>

#include "../_component/components/TransformComponent.h"
#include "../_component/components/CollisionComponent.h"

class ChunkManager : public Manager<ChunkManager> {
public:
    friend class CRSP<ChunkManager>;
    using ChunkMap = std::unordered_map<Chunk, EntityIdSet, PairIntIntHash>;
    using EntityMap = std::unordered_map<EntityId, ChunkSet>;
    using PositionCache = std::unordered_map<EntityId, std::pair<sf::Vector2f, sf::Vector2f>>;
    using ChunkLengthType = float;
    using GroupedEntities = std::unordered_map<Chunk, std::set<Entity*>, PairIntIntHash>;
private:
    ChunkManager();
    ~ChunkManager();

public:

    void start_up() override;
    void shut_down() override;

    std::vector<Chunk> get_relevant_chunks() const;

    ChunkSet get_chunks_of(EntityId id) const;

    ChunkSet calculate_chunks(Entity* entity) const;

    Chunk get_chunk_from_position(float x, float y) const;

    bool share_chunks(EntityId a, EntityId b) const;

    void update_entity_chunks();

    void update_single_entity(Entity* entity);

    void clear();

    /*
    There are not const because it might require to generate vectors for new chunks without entities
    */
    void draw_debug_chunks();

    void draw_debug_chunk_configuration();

    const EntityIdSet& get_entities_of_chunk(Chunk chunk);

    const std::set<Entity*>& get_colliding_entities_of_chunk(Chunk chunk);

    /**
    Returns a vector of vectors with the entities grouped by chunks, that is, in each vector
    we have all the entities that are in a particular chunk, entities can be repeated in more
    than one vector.
    */
    template<typename... T>
    GroupedEntities get_grouped_entities_with_components() {

        auto general = GroupedEntities{};

        for (auto & c : m_chunk_map) {

            if (c.second.size() == 0) continue;

            auto& entity_set = c.second;

            auto entities_in_chunk = std::vector<Entity*>{};

            for (auto& e : entity_set) {
                auto entity = EntityManager::get().get_entity(e);

                if (!entity->is_in_relevant_chunk()) continue;

                auto checks = { entity->has_component<T>()... };

                if (std::all_of(checks.begin(), checks.end(), [](bool i) {return i; })) {
                    entities_in_chunk.push_back(entity);
                }
            }

            general[c.first] = (std::move(entities_in_chunk));

        }
        return general;
    }

    /*
    Specific version of this function for the elements with a collision component
    that returns an internal member that we calculate on our "update_entity_chunks" instead
    of calculating it on request like we are doing in the function above.
    */
    template<>
    GroupedEntities get_grouped_entities_with_components<TransformComponent, CollisionComponent>() {
        return m_collision_components_cache;
    }

private:

    /**
    The length of the side of each square chunk
    */
    ChunkLengthType m_chunk_size{ 16 };

    /**
    Number of chunks outside the visible screen that we still consider relevant
    */
    unsigned int m_chunk_threshold{ 1 };

    /**
    This is used in case an entity enters a chunk exactly in the frame that we have to check
    something in that new chunk, we register the entity as being in the chunks that it
    really is plus "m_safe_threshold" more in each direction.
    */
    unsigned int m_safe_threshold{ 1 };

    ChunkMap m_chunk_map{};
    EntityMap m_entity_map{};
    PositionCache m_location_cache{};

    Chunk m_min_relevant_chunk{};
    Chunk m_max_relevant_chunk{};

    GroupedEntities m_collision_components_cache{};

    sf::Vector2f camera_previous_center{ std::numeric_limits<float>::infinity(),std::numeric_limits<float>::infinity() };
    sf::Vector2f camera_previous_size{ std::numeric_limits<float>::infinity(),std::numeric_limits<float>::infinity() };
    bool m_camera_changed{ true };

};

max_stars_repo_name: osor-io/NGene

max_issues_repo_name: osor-io/NGene

max_forks_repo_name: osor-io/NGene


text:
//IN method.
//Bilinear transformation (3D tensor of weights, 1D vector of biases)
//of Ni1 inputs from X1 and Ni2 inputs from X2 to No outputs in Y.
//This version uses CBLAS.

//Input X1 has Ni1 neurons, X2 has Ni2 neurons, and output Y has No neurons.

//If col-major: Y[o,l] = X1[:,l]' * W[:,:,o] * X2[:,l] + B[o]
//where:
//X1 has size Ni1 x L
//X2 has size Ni2 x L
//Y  has size No x L
//W  has size Ni1 x Ni2 x No
//B  has size No x 1

//If row-major: Y[l,o] = X1[l,:] * W[o,:,:] * X2[l,:]' + B[o]
//X1 has size L x Ni1
//X2 has size L x Ni2
//Y  has size L x No
//W  has size No x Ni2 x Ni1
//B  has size 1 x No

//This is equal to the outer product of X1 and X2 at each l (time point),
//and then take the weighted sum by each of the No matrices within W,
//producing No scalars, one for each output neuron in Y.

//I retain the for loop through L for compatibility with real-time streaming.

//This CBLAS version not working yet!!

#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int bilinear_cblas_s (float *Y, const float *X1, const float *X2, const float *W, const float *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L);
int bilinear_cblas_d (double *Y, const double *X1, const double *X2, const double *W, const double *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L);
int bilinear_cblas_c (float *Y, const float *X1, const float *X2, const float *W, const float *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L);
int bilinear_cblas_z (double *Y, const double *X1, const double *X2, const double *W, const double *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L);

int bilinear_cblas_s (float *Y, const float *X1, const float *X2, const float *W, const float *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L)
{
    const size_t N12 = Ni1*Ni2, Nw = N12*No;

    float *outer_prod;
    outer_prod = (float *)malloc(N12*sizeof(float));
    if (!outer_prod) { fprintf(stderr,"error in bilinear_cblas_s: problem with malloc. "); perror("malloc"); return 1; }

    for (size_t l=L; l>0u; --l, X1+=Ni1, X2+=Ni2, B-=No, W-=Nw)
    {
        //Re-zero outer_prod
        for (size_t n=N12; n>0u; --n) { outer_prod[n] = 0.0f; }

        //Outer product of X1, X2
        cblas_sger(CblasRowMajor,(int)Ni1,(int)Ni2,1.0f,X1,1,X2,1,outer_prod,(int)Ni2);

        //Weight by W for each output
        for (size_t o=No; o>0u; --o, ++B, ++Y, W+=N12)
        {
            *Y = *B + cblas_sdot((int)N12,W,1,outer_prod,1);
            //*Y = cblas_sdsdot((int)N12,*B,W,1,outer_prod,1);
        }
    }

    free(outer_prod);

    return 0;
}

int bilinear_cblas_d (double *Y, const double *X1, const double *X2, const double *W, const double *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L)
{
    const size_t N12 = Ni1*Ni2, Nw = N12*No;

    double *outer_prod;
    outer_prod = (double *)malloc(N12*sizeof(double));
    if (!outer_prod) { fprintf(stderr,"error in bilinear_cblas_d: problem with malloc. "); perror("malloc"); return 1; }

    for (size_t l=L; l>0u; --l, X1+=Ni1, X2+=Ni2, B-=No, W-=Nw)
    {
        //Re-zero outer_prod
        for (size_t n=N12; n>0u; --n) { outer_prod[n] = 0.0; }

        //Outer product of X1, X2
        cblas_dger(CblasRowMajor,(int)Ni1,(int)Ni2,1.0,X1,1,X2,1,outer_prod,(int)Ni2);

        //Weight by W for each output
        for (size_t o=No; o>0u; --o, ++B, ++Y, W+=N12)
        {
            *Y = *B + cblas_ddot((int)N12,W,1,outer_prod,1);
        }
    }

    free(outer_prod);

    return 0;
}

int bilinear_cblas_c (float *Y, const float *X1, const float *X2, const float *W, const float *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L)
{
    const size_t N12 = 2u*Ni1*Ni2, Nw = N12*No;
    const float oz[2] = {1.0f,0.0f};

    float *outer_prod;
    outer_prod = (float *)malloc(N12*sizeof(float));
    if (!outer_prod) { fprintf(stderr,"error in bilinear_cblas_c: problem with malloc. "); perror("malloc"); return 1; }

    for (size_t l=L; l>0u; --l, X1+=2u*Ni1, X2+=2u*Ni2, B-=2u*No, W-=Nw)
    {
        //Re-zero outer_prod
        for (size_t n=N12; n>0u; --n) { outer_prod[n] = 0.0f; }

        //Outer product of X1, X2
        cblas_cgeru(CblasRowMajor,(int)Ni1,(int)Ni2,oz,X1,1,X2,1,outer_prod,(int)Ni2);

        //Weight by W for each output
        for (size_t o=No; o>0u; --o, W+=N12)
        {
            cblas_cdotu_sub((int)N12,W,1,outer_prod,1,Y);
            *Y++ += *B++; *Y++ += *B++;
        }
    }

    free(outer_prod);

    return 0;
}

int bilinear_cblas_z (double *Y, const double *X1, const double *X2, const double *W, const double *B, const size_t Ni1, const size_t Ni2, const size_t No, const size_t L)
{
    const size_t N12 = 2u*Ni1*Ni2, Nw = N12*No;
    const double oz[2] = {1.0,0.0};

    double *outer_prod;
    outer_prod = (double *)malloc(N12*sizeof(double));
    if (!outer_prod) { fprintf(stderr,"error in bilinear_cblas_z: problem with malloc. "); perror("malloc"); return 1; }

    for (size_t l=L; l>0u; --l, X1+=2u*Ni1, X2+=2u*Ni2, B-=2u*No, W-=Nw)
    {
        //Re-zero outer_prod
        for (size_t n=N12; n>0u; --n) { outer_prod[n] = 0.0; }

        //Outer product of X1, X2
        cblas_zgeru(CblasRowMajor,(int)Ni1,(int)Ni2,oz,X1,1,X2,1,outer_prod,(int)Ni2);

        //Weight by W for each output
        for (size_t o=No; o>0u; --o, W+=N12)
        {
            cblas_zdotu_sub((int)N12,W,1,outer_prod,1,Y);
            *Y++ += *B++; *Y++ += *B++;
        }
    }

    free(outer_prod);

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#ifndef _BLAS_UTIL_H_
#define _BLAS_UTIL_H_

#ifdef __BLAS_LEGACY__
#include <math.h>
#include "utils/cblas.h"
//#include <lapacke.h> //! hasn't used
#define _D2_MALLOC_SCALAR(x)       (SCALAR *) malloc( (x) *sizeof(SCALAR))
#define _D2_MALLOC_INT(x)       (int *) malloc( (x) *sizeof(int))
#define _D2_MALLOC_SIZE_T(x)       (size_t *) malloc( (x) *sizeof(size_t))
#define _D2_CALLOC_SCALAR(x)       (SCALAR *) calloc( (x) , sizeof(SCALAR))
#define _D2_CALLOC_INT(x)       (int *) calloc( (x) , sizeof(int))
#define _D2_CALLOC_SIZE_T(x)       (size_t *) calloc( (x) , sizeof(size_t))
#define _D2_FREE(x)         free(x)

#elif defined __APPLE__
#include <Accelerate/Accelerate.h>
#define _D2_MALLOC_SCALAR(x)       (SCALAR *) malloc( (x) *sizeof(SCALAR))
#define _D2_MALLOC_INT(x)       (int *) malloc( (x) *sizeof(int))
#define _D2_MALLOC_SIZE_T(x)       (size_t *) malloc( (x) *sizeof(size_t))
#define _D2_CALLOC_SCALAR(x)       (SCALAR *) calloc( (x) , sizeof(SCALAR))
#define _D2_CALLOC_INT(x)       (int *) calloc( (x) , sizeof(int))
#define _D2_CALLOC_SIZE_T(x)       (size_t *) calloc( (x) , sizeof(size_t))
#define _D2_FREE(x)         free(x)

#elif defined __USE_MKL__
#include <mkl.h>
#define _D2_MALLOC_SCALAR(x)       (SCALAR *) mkl_malloc( (x) *sizeof(SCALAR), 16)
#define _D2_MALLOC_INT(x)       (int *) mkl_malloc( (x) *sizeof(int), 16)
#define _D2_MALLOC_SIZE_T(x)       (size_t *) mkl_malloc( (x) *sizeof(size_t), 16)
#define _D2_CALLOC_SCALAR(x)       (SCALAR *) mkl_calloc( (x) , sizeof(SCALAR), 16)
#define _D2_CALLOC_INT(x)       (int *) mkl_calloc( (x) , sizeof(int), 16)
#define _D2_CALLOC_SIZE_T(x)       (size_t *) mkl_calloc( (x) , sizeof(size_t), 16)
#define _D2_FREE(x)         mkl_free(x)

#endif

#endif /* _BLAS_UTIL_H_ */


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <assert.h>

#include "core_allvars.h"
#include "core_proto.h"

void init_galaxy(int p, int halonr)
{
  int j, step;

	assert(halonr == Halo[halonr].FirstHaloInFOFgroup);

  Gal[p].Type = 0;

  Gal[p].GalaxyNr = GalaxyCounter;
  GalaxyCounter++;

  Gal[p].HaloNr = halonr;
  Gal[p].MostBoundID = Halo[halonr].MostBoundID;
  Gal[p].SnapNum = Halo[halonr].SnapNum - 1;

  Gal[p].mergeType = 0;
  Gal[p].mergeIntoID = -1;
  Gal[p].mergeIntoSnapNum = -1;
  Gal[p].dT = -1.0;

  for(j = 0; j < 3; j++)
  {
    Gal[p].Pos[j] = Halo[halonr].Pos[j];
    Gal[p].Vel[j] = Halo[halonr].Vel[j];
  }

  Gal[p].Len = Halo[halonr].Len;
  Gal[p].Vmax = Halo[halonr].Vmax;
  Gal[p].Vvir = get_virial_velocity(halonr);
  Gal[p].Mvir = get_virial_mass(halonr);
  Gal[p].Rvir = get_virial_radius(halonr);

  Gal[p].deltaMvir = 0.0;

  Gal[p].ColdGas = 0.0;
  Gal[p].StellarMass = 0.0;
  Gal[p].BulgeMass = 0.0;
  Gal[p].HotGas = 0.0;
  Gal[p].EjectedMass = 0.0;
  Gal[p].BlackHoleMass = 0.0;
  Gal[p].ICS = 0.0;

  Gal[p].MetalsColdGas = 0.0;
  Gal[p].MetalsStellarMass = 0.0;
  Gal[p].MetalsBulgeMass = 0.0;
  Gal[p].MetalsHotGas = 0.0;
  Gal[p].MetalsEjectedMass = 0.0;
  Gal[p].MetalsICS = 0.0;

  for(step = 0; step < STEPS; step++)
  {
    Gal[p].SfrDisk[step] = 0.0;
    Gal[p].SfrBulge[step] = 0.0;
    Gal[p].SfrDiskColdGas[step] = 0.0;
    Gal[p].SfrDiskColdGasMetals[step] = 0.0;
    Gal[p].SfrBulgeColdGas[step] = 0.0;
    Gal[p].SfrBulgeColdGasMetals[step] = 0.0;
  }

  Gal[p].DiskScaleRadius = get_disk_radius(halonr, p);
  Gal[p].MergTime = 999.9;
  Gal[p].Cooling = 0.0;
  Gal[p].Heating = 0.0;
  Gal[p].r_heat = 0.0;
  Gal[p].QuasarModeBHaccretionMass = 0.0;
  Gal[p].TimeSinceMajorMerger = -1.0;
  Gal[p].TimeSinceMinorMerger = -1.0;
  Gal[p].OutflowRate = 0.0;
	Gal[p].TotalSatelliteBaryons = 0.0;

	// infall properties
  Gal[p].infallMvir = -1.0;
  Gal[p].infallVvir = -1.0;
  Gal[p].infallVmax = -1.0;

  //Raouf New things add to code
  Gal[p].Qjet = 0.0;
  Gal[p].Rcocoon = 0.0;
  Gal[p].Rshocked = 0.0;
  Gal[p].t_AGN_returne = 0.0;
  Gal[p].t_AGN_on = 0.0;
  Gal[p].Tshocked = 0.0;
  Gal[p].Mshocked = 0.0;
  for(j = 0; j < 7; j++)
  {
      Gal[p].RadioLuminosity[j] = 0.0;
  }
  for(j = 0; j < 600; j++)
  {
      Gal[p].RadioLuminosity_lifetime[j] = -999.0;
  }
  for(j = 0; j < 600; j++)
  {
      Gal[p].Rshocked_lifetime[j] = -999.0;
  }

  Gal[p].RadioAGNaccretionRate = 0.0;
  Gal[p].rho_zero_Makino = 0.0;
  Gal[p].rho_zero_Capelo = 0.0;
  Gal[p].rho_zero_iso = 0.0;
  Gal[p].b_gas = 0.0;
  Gal[p].Rs = 0.0;
  Gal[p].concentration = 0.0;
  Gal[p].Temp_Gas = 0.0;
  Gal[p].Lx_bol = 0.0;
  Gal[p].R_index = 0.0;
  Gal[p].Q_index = 0.0;
  Gal[p].R_cool = 0.0;
  Gal[p].fcool = 0.0;
  Gal[p].t_static = 0.0;
  Gal[p].t_AGN_off = 0.0;
  Gal[p].time_to_next_on = 0.0;
  Gal[p].delta = 0.0;
  Gal[p].t_cool_Makino = 0.0;
}

double get_disk_radius(int halonr, int p)
{
  double SpinMagnitude, SpinParameter;

	if(Gal[p].Vvir > 0.0 && Gal[p].Rvir > 0.0)
	{
		// See Mo, Shude & White (1998) eq12, and using a Bullock style lambda.
		SpinMagnitude = sqrt(Halo[halonr].Spin[0] * Halo[halonr].Spin[0] +
			Halo[halonr].Spin[1] * Halo[halonr].Spin[1] + Halo[halonr].Spin[2] * Halo[halonr].Spin[2]);

		// trim the extreme tail of the spin distribution for more a realistic r_s
		if(SpinMagnitude > 1.5)
			SpinMagnitude = 1.5;

		SpinParameter = SpinMagnitude / (1.414 * Gal[p].Vvir * Gal[p].Rvir);
		return (SpinParameter / 1.414) * Gal[p].Rvir;
	}
	else
		return 0.1 * Gal[p].Rvir;

}

double get_metallicity(double gas, double metals)
{
  double metallicity;

  if(gas > 0.0 && metals > 0.0)
  {
    metallicity = metals / gas;
    if(metallicity < 1.0)
      return metallicity;
    else
      return 1.0;
  }
  else
    return 0.0;

}

double dmax(double x, double y)
{
  if(x > y)
    return x;
  else
    return y;
}

double get_virial_mass(int halonr)
{
  if(halonr == Halo[halonr].FirstHaloInFOFgroup && Halo[halonr].Mvir >= 0.0)
    return Halo[halonr].Mvir;   /* take spherical overdensity mass estimate */
  else
    return Halo[halonr].Len * PartMass;
}

double get_virial_velocity(int halonr)
{
	double Rvir;

	Rvir = get_virial_radius(halonr);

  if(Rvir > 0.0)
		return sqrt(G * get_virial_mass(halonr) / Rvir);
	else
		return 0.0;
}

double get_virial_radius(int halonr)
{
  // return Halo[halonr].Rvir;  // Used for Bolshoi

  double zplus1, hubble_of_z_sq, rhocrit, fac;

  zplus1 = 1 + ZZ[Halo[halonr].SnapNum];
  hubble_of_z_sq =
    Hubble * Hubble *(Omega * zplus1 * zplus1 * zplus1 + (1 - Omega - OmegaLambda) * zplus1 * zplus1 +
    OmegaLambda);

  rhocrit = 3 * hubble_of_z_sq / (8 * M_PI * G);
  fac = 1 / (200 * 4 * M_PI / 3.0 * rhocrit);

  return cbrt(get_virial_mass(halonr) * fac);
}

max_stars_repo_name: mojtabaraouf/sage

max_issues_repo_name: mojtabaraouf/sage

max_forks_repo_name: mojtabaraouf/sage