text:
/* multifit/multilinear.c
 *
 * Copyright (C) 2000, 2007, 2010 Brian Gough
 * Copyright (C) 2013, 2015 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.
 */

//#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 "linear_common.c"

static int multifit_linear_svd (const gsl_matrix * X,
                                const int balance,
                                gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear (const gsl_matrix * X,
                     const gsl_vector * y,
                     gsl_vector * c,
                     gsl_matrix * cov,
                     double *chisq, gsl_multifit_linear_workspace * work)
{
  size_t rank;
  int status = gsl_multifit_linear_tsvd(X, y, GSL_DBL_EPSILON, c, cov, chisq, &rank, work);

  return status;
}

/*
gsl_multifit_linear_tsvd()
  Solve linear least squares system with truncated SVD
Inputs: X     - least squares matrix, n-by-p
        y     - right hand side vector, n-by-1
        tol   - tolerance for singular value truncation; if
                s_j <= tol * s_0
                then it is discarded from series expansion
        c     - (output) solution vector, p-by-1
        cov   - (output) covariance matrix, p-by-p
        chisq - (output) cost function chi^2
        rank  - (output) effective rank (number of singular values
                used in solution)
        work  - workspace
*/

int
gsl_multifit_linear_tsvd (const gsl_matrix * X,
                          const gsl_vector * y,
                          const double tol,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double * chisq,
                          size_t * rank,
                          gsl_multifit_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (y->size != n)
    {
      GSL_ERROR("number of observations in y does not match matrix",
                GSL_EBADLEN);
    }
  else if (p != c->size)
    {
      GSL_ERROR ("number of parameters c does not match matrix",
                 GSL_EBADLEN);
    }
  else if (tol <= 0)
    {
      GSL_ERROR ("tolerance must be positive", GSL_EINVAL);
    }
  else
    {
      int status;
      double rnorm = 0.0, snorm;

      /* compute balanced SVD */
      status = gsl_multifit_linear_bsvd (X, work);
      if (status)
        return status;

      status = multifit_linear_solve (X, y, tol, -1.0, rank,
                                      c, &rnorm, &snorm, work);

      *chisq = rnorm * rnorm;

      /* variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */
      {
        double r2 = rnorm * rnorm;
        double s2 = r2 / (double)(n - *rank);
        size_t i, j;
        gsl_matrix_view QSI = gsl_matrix_submatrix(work->QSI, 0, 0, p, p);
        gsl_vector_view D = gsl_vector_subvector(work->D, 0, p);

        for (i = 0; i < p; i++)
          {
            gsl_vector_view row_i = gsl_matrix_row (&QSI.matrix, i);
            double d_i = gsl_vector_get (&D.vector, i);

            for (j = i; j < p; j++)
              {
                gsl_vector_view row_j = gsl_matrix_row (&QSI.matrix, j);
                double d_j = gsl_vector_get (&D.vector, j);
                double s;

                gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);

                gsl_matrix_set (cov, i, j, s * s2 / (d_i * d_j));
                gsl_matrix_set (cov, j, i, s * s2 / (d_i * d_j));
              }
          }
      }

      return status;
    }
}

/*
gsl_multifit_linear_svd()
  Perform SVD decomposition of the matrix X and store in work without
balancing
*/

int
gsl_multifit_linear_svd (const gsl_matrix * X,
                         gsl_multifit_linear_workspace * work)
{
  /* do not balance by default */
  int status = multifit_linear_svd(X, 0, work);

  return status;
}

/*
gsl_multifit_linear_bsvd()
  Perform SVD decomposition of the matrix X and store in work with
balancing
*/

int
gsl_multifit_linear_bsvd (const gsl_matrix * X,
                          gsl_multifit_linear_workspace * work)
{
  int status = multifit_linear_svd(X, 1, work);

  return status;
}

size_t
gsl_multifit_linear_rank(const double tol, const gsl_multifit_linear_workspace * work)
{
  double s0 = gsl_vector_get (work->S, 0);
  size_t rank = 0;
  size_t j;

  for (j = 0; j < work->p; j++)
    {
      double sj = gsl_vector_get (work->S, j);

      if (sj > tol * s0)
        ++rank;
    }

  return rank;
}

/* Estimation of values for given x */

int
gsl_multifit_linear_est (const gsl_vector * x,
                         const gsl_vector * c,
                         const gsl_matrix * cov, double *y, double *y_err)
{

  if (x->size != c->size)
    {
      GSL_ERROR ("number of parameters c does not match number of observations x",
         GSL_EBADLEN);
    }
  else if (cov->size1 != cov->size2)
    {
      GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
    }
  else if (c->size != cov->size1)
    {
      GSL_ERROR ("number of parameters c does not match size of covariance matrix cov",
         GSL_EBADLEN);
    }
  else
    {
      size_t i, j;
      double var = 0;

      gsl_blas_ddot(x, c, y);       /* y = x.c */

      /* var = x' cov x */

      for (i = 0; i < x->size; i++)
        {
          const double xi = gsl_vector_get (x, i);
          var += xi * xi * gsl_matrix_get (cov, i, i);

          for (j = 0; j < i; j++)
            {
              const double xj = gsl_vector_get (x, j);
              var += 2 * xi * xj * gsl_matrix_get (cov, i, j);
            }
        }

      *y_err = sqrt (var);

      return GSL_SUCCESS;
    }
}

/*
gsl_multifit_linear_rcond()
  Return reciprocal condition number of LS matrix;
gsl_multifit_linear_svd() must first be called to
compute the SVD of X and its reciprocal condition number
*/

double
gsl_multifit_linear_rcond (const gsl_multifit_linear_workspace * w)
{
  return w->rcond;
}

/*
gsl_multifit_linear_residuals()
  Compute vector of residuals from fit
Inputs: X - design matrix
        y - rhs vector
        c - fit coefficients
        r - (output) where to store residuals
*/

int
gsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
                               const gsl_vector *c, gsl_vector *r)
{
  if (X->size1 != y->size)
    {
      GSL_ERROR
        ("number of observations in y does not match rows of matrix X",
         GSL_EBADLEN);
    }
  else if (X->size2 != c->size)
    {
      GSL_ERROR ("number of parameters c does not match columns of matrix X",
                 GSL_EBADLEN);
    }
  else if (y->size != r->size)
    {
      GSL_ERROR ("number of observations in y does not match number of residuals",
                 GSL_EBADLEN);
    }
  else
    {
      /* r = y - X c */
      gsl_vector_memcpy(r, y);
      gsl_blas_dgemv(CblasNoTrans, -1.0, X, c, 1.0, r);

      return GSL_SUCCESS;
    }
} /* gsl_multifit_linear_residuals() */

/* Perform a SVD decomposition on the least squares matrix X = U S Q^T
 *
 * Inputs: X       - least squares matrix
 *         balance - 1 to perform column balancing
 *         work    - workspace
 *
 * Notes:
 * 1) On output,
 *    work->A contains the matrix U
 *    work->Q contains the matrix Q
 *    work->S contains the vector of singular values
 * 2) The matrix X may have smaller dimensions than the workspace
 *    in the case of stdform2() - but the dimensions cannot be larger
 * 3) On output, work->n and work->p are set to the dimensions of X
 * 4) On output, work->rcond is set to the reciprocal condition number of X
 */

static int
multifit_linear_svd (const gsl_matrix * X,
                     const int balance,
                     gsl_multifit_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (n > work->nmax || p > work->pmax)
    {
      GSL_ERROR("observation matrix larger than workspace", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
      gsl_matrix_view Q = gsl_matrix_submatrix(work->Q, 0, 0, p, p);
      gsl_matrix_view QSI = gsl_matrix_submatrix(work->QSI, 0, 0, p, p);
      gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
      gsl_vector_view xt = gsl_vector_subvector(work->xt, 0, p);
      gsl_vector_view D = gsl_vector_subvector(work->D, 0, p);

      /* Copy X to workspace,  A <= X */

      gsl_matrix_memcpy (&A.matrix, X);

      /* Balance the columns of the matrix A if requested */

      if (balance)
        {
          gsl_linalg_balance_columns (&A.matrix, &D.vector);
        }
      else
        {
          gsl_vector_set_all (&D.vector, 1.0);
        }

      /* decompose A into U S Q^T */
      gsl_linalg_SV_decomp_mod (&A.matrix, &QSI.matrix, &Q.matrix,
                                &S.vector, &xt.vector);

      /* compute reciprocal condition number rcond = smin / smax */
      {
        double smin, smax;
        gsl_vector_minmax(&S.vector, &smin, &smax);
        work->rcond = smin / smax;
      }

      work->n = n;
      work->p = p;

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: agroimpacts/imager

max_issues_repo_name: agroimpacts/imager

max_forks_repo_name: agroimpacts/imager


text:
/* multimin/fminimizer.c
 *
 * 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 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multimin.h>

gsl_multimin_fminimizer *
gsl_multimin_fminimizer_alloc (const gsl_multimin_fminimizer_type * T,
                               size_t n)
{
  int status;

  gsl_multimin_fminimizer *s =
    (gsl_multimin_fminimizer *) malloc (sizeof (gsl_multimin_fminimizer));

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

  s->type = T;

  s->x = gsl_vector_calloc (n);

  if (s->x == 0)
    {
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
    }

  s->state = malloc (T->size);

  if (s->state == 0)
    {
      gsl_vector_free (s->x);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for minimizer state",
                     GSL_ENOMEM, 0);
    }

  status = (T->alloc) (s->state, n);

  if (status != GSL_SUCCESS)
    {
      free (s->state);
      gsl_vector_free (s->x);
      free (s);

      GSL_ERROR_VAL ("failed to initialize minimizer state", GSL_ENOMEM, 0);
    }

  return s;
}

int
gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * s,
                             gsl_multimin_function * f,
                             const gsl_vector * x,
                             const gsl_vector * step_size)
{
  if (s->x->size != f->n)
    {
      GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN);
    }

  if (x->size != f->n || step_size->size != f->n)
    {
      GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN);
    }

  s->f = f;

  gsl_vector_memcpy (s->x,x);

  return (s->type->set) (s->state, s->f, s->x, &(s->size), step_size);
}

void
gsl_multimin_fminimizer_free (gsl_multimin_fminimizer * s)
{
  (s->type->free) (s->state);
  free (s->state);
  gsl_vector_free (s->x);
  free (s);
}

int
gsl_multimin_fminimizer_iterate (gsl_multimin_fminimizer * s)
{
  return (s->type->iterate) (s->state, s->f, s->x, &(s->size), &(s->fval));
}

const char *
gsl_multimin_fminimizer_name (const gsl_multimin_fminimizer * s)
{
  return s->type->name;
}

gsl_vector *
gsl_multimin_fminimizer_x (const gsl_multimin_fminimizer * s)
{
  return s->x;
}

double
gsl_multimin_fminimizer_minimum (const gsl_multimin_fminimizer * s)
{
  return s->fval;
}

double
gsl_multimin_fminimizer_size (const gsl_multimin_fminimizer * s)
{
  return s->size;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
** TFCE heuristic
**
** M. Kuhlmann, MPI-KYB, Sept 2015
*/

#include <viaio/VImage.h>
#include <viaio/Vlib.h>
#include <viaio/mu.h>
#include <via/via.h>

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

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

#define ABS(x) ((x) > 0 ? (x) : -(x))

VImage tfce(VImage t_image,VImage map,int nstrata)
{
  int b,r,c,i,j,k;
  size_t nslices = VImageNBands(t_image);
  size_t nrows  = VImageNRows(t_image);
  size_t ncols  = VImageNColumns(t_image);
  double tiny = 1e-16;

  /* if pixel map is not given, compute it */
  if(!map) {
    VImage mask = VCreateImage(nslices,nrows,ncols,VBitRepn);
    VFillImage(mask,VAllBands,0);

    size_t nvox = 0;
    float u;
    for (b=0; b<nslices; b++) {
      for (r=0; r<nrows; r++) {
	for (c=0; c<ncols; c++) {
	  u = VPixel(t_image,b,r,c,VFloat);
	  if (ABS(u) >= tiny) {
	    VPixel(mask,b,r,c,VBit) = 1;
	    nvox++;
	  }
	}
      }
    }

    if(nvox == 0)
      VError("TFCE: Number of voxels!=0 is zero.");
    /* voxel addresses */
    map = VCreateImage(1,3,nvox,VShortRepn);
    if (map == NULL) VError(" error allocating addr map");
    VFillImage(map,VAllBands,0);
    i = 0;
    for (b=0; b<nslices; b++) {
      for (r=0; r<nrows; r++) {
	for (c=0; c<ncols; c++) {
	  if (VGetPixel(mask,b,r,c) < 0.5) continue;
	  VPixel(map,0,0,i,VShort) = b;
	  VPixel(map,0,1,i,VShort) = r;
	  VPixel(map,0,2,i,VShort) = c;
	  i++;
	}
      }
    }
    VDestroyImage(mask);
  }

  /* stratification */
  /* set z strata */
  size_t nvox = VImageNColumns(map);
  double z;

  /* get an idea of z-value distribution */
  double zmax=0;
  for (i=0; i<nvox; i++) {
    b = VPixel(map,0,0,i,VShort);
    r = VPixel(map,0,1,i,VShort);
    c = VPixel(map,0,2,i,VShort);
    z = VPixel(t_image,b,r,c,VFloat);
    if (z > zmax) zmax = z;
  }

  /* ini strata */
  double step = zmax / (double)nstrata;
  double *stratum = (double *) VCalloc(nstrata+1,sizeof(double));
  for (j=0; j<=nstrata; j++) stratum[j] = step*j;

  float *clustersize = (float *) VCalloc(nvox,sizeof(float));
  if (!clustersize) VError(" err allocating clustersize");
  float *cluster_cont = (float *) VCalloc(nvox,sizeof(float));
  if (!cluster_cont) VError(" err allocating cluster_cont");
  for (i=0; i<nvox; i++) {
    clustersize[i] = 0;
    cluster_cont[i] = 0;
  }

  /* output, create TFCE image */
  VImage tfce  = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  if (!tfce) VError(" err allocating tfce image");
  VFillImage(tfce,VAllBands,0);

  /* ini images */
  VImage bin_image = VCreateImage(nslices,nrows,ncols,VBitRepn);
  if (!bin_image) VError(" err allocating bin_image image");
  VImage label_image = VCreateImage(nslices,nrows,ncols,VShortRepn);
  if (!label_image) VError(" err allocating label_image image");

  /* for each stratum */
  double *strat_sq = (double *)VCalloc(nstrata,sizeof(double));
  for (j=0; j<nstrata; j++)
    strat_sq[j] = stratum[j] * stratum[j];
  for (j=0; j<nstrata; j++) {
    VFillImage(bin_image,VAllBands,0);
    for (i=0; i<nvox; i++) {
      b = VPixel(map,0,0,i,VShort);
      r = VPixel(map,0,1,i,VShort);
      c = VPixel(map,0,2,i,VShort);
      z = VPixel(t_image,b,r,c,VFloat);
      if (z > stratum[j]) {
	VPixel(bin_image,b,r,c,VBit) = 1;
      }
    }

    /* get connected components (clusters) */
    long nl=0;
    VLabelImage3d(bin_image,label_image,26,VShortRepn,&nl);
    if (nl < 1) continue;   /* no voxels in this stratum found */

    /* compute cluster sizes */
    int max_k = 0;
    for (i=0; i<nvox; i++) clustersize[i] = 0;
    for (i=0; i<nvox; i++) {
      b = VPixel(map,0,0,i,VShort);
      r = VPixel(map,0,1,i,VShort);
      c = VPixel(map,0,2,i,VShort);
      k = VPixel(label_image,b,r,c,VShort);
      if (k > 0) clustersize[k]++;
      if (k > max_k) max_k = k;
    }

    /* compute contribution of the cluster to the voxels in it */
    for (i=0; i<=max_k; i++)
      cluster_cont[i] = sqrt(clustersize[i]) * strat_sq[j];

    /* add cluster contribution to each voxel */
    for (i=0; i<nvox; i++) {
      b = VPixel(map,0,0,i,VShort);
      r = VPixel(map,0,1,i,VShort);
      c = VPixel(map,0,2,i,VShort);
      k = VPixel(label_image,b,r,c,VShort);
      VPixel(tfce,b,r,c,VFloat) += cluster_cont[k];
    }
  }

  /* thickness correction */
  for (i=0; i<nvox; i++) {
    b = VPixel(map,0,0,i,VShort);
    r = VPixel(map,0,1,i,VShort);
    c = VPixel(map,0,2,i,VShort);
    k = VPixel(label_image,b,r,c,VShort);
    VPixel(tfce,b,r,c,VFloat) *= step;
  }

  VDestroyImage(bin_image);
  VDestroyImage(label_image);
  VFree(clustersize);
  VFree(cluster_cont);
  VFree(stratum);
  VFree(strat_sq);

  return tfce;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira

 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, 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.

 $Id$
*/

#include <config.h>

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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_combination.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_chebyshev.h>
#include <gsl/gsl_deriv.h>

#include "string_f.h"

#include <fortran_types.h>

/* Numerical threshold for oct_bessel_k0 and oct_bessel_k1 */
#define  BESSEL_K_THRES  1.0e2

/* ---------------------- Interface to GSL functions ------------------------ */

/* Mathematical Functions */
double FC_FUNC_(oct_asinh, OCT_ASINH)
     (const double *x)
{
  return gsl_asinh(*x);
}

/* Special Functions */
double FC_FUNC_(oct_gamma, OCT_GAMMA)
     (const double *x)
{
  return gsl_sf_gamma(*x);
}

double FC_FUNC_(oct_hypergeometric, OCT_HYPERGEOMETRIC)
     (const double *a, const double*b, const double *x)
{
  return gsl_sf_hyperg_U(*a, *b, *x);
}

double FC_FUNC_(oct_incomplete_gamma, OCT_INCOMPLETE_GAMMA)
     (const double *a, const double *x)
{
  return gsl_sf_gamma_inc_Q(*a, *x);
}

double FC_FUNC_(oct_sph_bessel, OCT_SPH_BESSEL)
     (const fint *l, const double*x)
{
  return gsl_sf_bessel_jl(*l, *x);
}

double FC_FUNC_(oct_bessel, OCT_BESSEL)
     (const fint *n, const double *x)
{
  return gsl_sf_bessel_Jn(*n, *x);
}

double FC_FUNC_(oct_bessel_in, OCT_BESSEL_IN)
     (const fint *n, const double *x)
{
  return gsl_sf_bessel_In(*n, *x);
}

double FC_FUNC_(oct_bessel_j0, OCT_BESSEL_J0)
     (const double *x)
{
  return gsl_sf_bessel_J0(*x);
}

double FC_FUNC_(oct_bessel_j1, OCT_BESSEL_J1)
     (const double *x)
{
  return gsl_sf_bessel_J1(*x);
}

double FC_FUNC_(oct_bessel_k0, OCT_BESSEL_K0)
     (const double *x)
{
  if( *x > BESSEL_K_THRES )
  {
    return 0.0e0;
  } else {
    return gsl_sf_bessel_K0(*x);
  }
}

double FC_FUNC_(oct_bessel_k1, OCT_BESSEL_K1)
     (const double *x)
{
  if( *x > BESSEL_K_THRES )
  {
    return 0.0e0;
  } else {
    return gsl_sf_bessel_K1(*x);
  }
}

/* the GSL headers specify double x, not const double x */
double FC_FUNC_(oct_erfc, OCT_ERFC)
     (const double *x)
{
  /* avoid floating invalids in the asymptotic limit */
  if(*x >  20.0) return  0.0;
  if(*x < -20.0) return  2.0;
  /* otherwise call gsl */
  return gsl_sf_erfc(*x);
}

/* the GSL headers specify double x, not const double x */
double FC_FUNC_(oct_erf, OCT_ERF)
     (const double *x)
{
  /* avoid floating invalids in the asymptotic limit */
  if(*x >  20.0) return  1.0;
  if(*x < -20.0) return -1.0;
  /* otherwise call gsl */
  return gsl_sf_erf(*x);
}

double FC_FUNC_(oct_legendre_sphplm, OCT_LEGENDRE_SPHPLM)
     (const fint *l, const int *m, const double *x)
{
  return gsl_sf_legendre_sphPlm(*l, *m, *x);
}

double FC_FUNC_(oct_sine_integral, OCT_SINE_INTEGRAL)
     (const double *x)
{
  return gsl_sf_Si(*x);
}

/* generalized Laguerre polynomials */
double FC_FUNC_(oct_sf_laguerre_n, OCT_SF_LAGUERRE_N)
     (const int *n, const double *a, const double *x)
{
  return gsl_sf_laguerre_n(*n, *a, *x);
}

/* Vectors and Matrices */

/* Permutations */

/* Combinations */

void FC_FUNC_(oct_combination_init, OCT_COMBINATION_INIT)
     (gsl_combination **c, const fint *n, const fint *k)
{
  *c = gsl_combination_calloc (*n, *k);
}

void FC_FUNC_(oct_combination_next, OCT_COMBINATION_NEXT)
     (gsl_combination **c, fint *next)
{
  *next = gsl_combination_next (((gsl_combination *)(*c)));
}

void FC_FUNC_(oct_get_combination, OCT_GET_COMBINATION)
     (gsl_combination **c, fint *comb)
{
  int i;
  for (i=0;i< ((gsl_combination *)(*c))->k; i++) {
    comb[i]=(fint)((gsl_combination *)(*c))->data[i];
  }
}

void FC_FUNC_(oct_combination_end, OCT_COMBINATION_END)
     (gsl_combination **c)
{
  gsl_combination_free (((gsl_combination *)(*c)));
}

/* Linear Algebra */

/* Random Number Generation */
void FC_FUNC_(oct_ran_init, OCT_RAN_INIT)
     (gsl_rng **r)
{
  gsl_rng_env_setup();
  *r = gsl_rng_alloc(gsl_rng_default);
}

void FC_FUNC_(oct_ran_end, OCT_RAN_END)
     (gsl_rng **r)
{
  gsl_rng_free(*r);
}

/* Random Number Distributions */
double FC_FUNC_(oct_ran_gaussian, OCT_RAN_GAUSSIAN)
     (const gsl_rng **r, const double *sigma)
{
  return gsl_ran_gaussian(*r, *sigma);
}

double FC_FUNC_(oct_ran_flat, OCT_RAN_FLAT)
     (const gsl_rng **r, const double *a, const double *b)
{
  return gsl_ran_flat(*r, *a, *b);
}

/* Interpolation */
void FC_FUNC_(oct_spline_end, OCT_SPLINE_END)
     (void **spl, void **acc)
{
  gsl_spline_free((gsl_spline *)(*spl));
  gsl_interp_accel_free((gsl_interp_accel *)(*acc));
}

void FC_FUNC_(oct_spline_fit, OCT_SPLINE_FIT)
		 (const fint *nrc, const double *x, const double *y, void **spl, void **acc)
{
  /* the GSL headers actually specify size_t instead of const int for nrc */
  *acc = (void *)gsl_interp_accel_alloc();
  *spl = (void *)gsl_spline_alloc(gsl_interp_cspline, *nrc);
  gsl_spline_init((gsl_spline *)(*spl), x, y, *nrc);
  fflush(stdout);
}

double FC_FUNC_(oct_spline_eval, OCT_SPLINE_EVAL)
     (const double *x, const void **spl, void **acc)
{
  /* the GSL headers specify double x instead of const double x */
  return gsl_spline_eval((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
}

void FC_FUNC_(oct_spline_eval_array, OCT_SPLINE_EVAL_ARRAY)
     (const fint * nn, double *xf, const void **spl, void **acc)
{
  int ii;
  for(ii = 0; ii < *nn; ii++){
    xf[ii] = gsl_spline_eval((gsl_spline *)(*spl), xf[ii], (gsl_interp_accel *)(*acc));
  }
}

void FC_FUNC_(oct_spline_eval_array4, OCT_SPLINE_EVAL_ARRAY4)
     (const fint * nn, float *xf, const void **spl, void **acc)
{
  int ii;
  for(ii = 0; ii < *nn; ii++){
    xf[ii] = (float) gsl_spline_eval((gsl_spline *)(*spl), (double) xf[ii], (gsl_interp_accel *)(*acc));
  }
}

/* use a stride of 2 to store into just the real part of a Fortran complex array */
void FC_FUNC_(oct_spline_eval_arrayz, OCT_SPLINE_EVAL_ARRAYZ)
     (const fint * nn, double *xf, const void **spl, void **acc)
{
  int ii;
  for(ii = 0; ii < *nn; ii++){
    xf[ii*2] = gsl_spline_eval((gsl_spline *)(*spl), xf[ii*2], (gsl_interp_accel *)(*acc));
  }
}

/* use a stride of 2 to store into just the real part of a Fortran complex array */
void FC_FUNC_(oct_spline_eval_arrayc, OCT_SPLINE_EVAL_ARRAYC)
     (const fint * nn, float *xf, const void **spl, void **acc)
{
  int ii;
  for(ii = 0; ii < *nn; ii++){
    xf[ii*2] = (float) gsl_spline_eval((gsl_spline *)(*spl), (double) xf[ii*2], (gsl_interp_accel *)(*acc));
  }
}

/* This function returns the number of points with which a spline
	 was constructed (the size component of the gsl_spline struct). */
fint FC_FUNC_(oct_spline_npoints, OCT_SPLINE_NPOINTS)
     (const void **spl)
{
  return (fint)((gsl_spline *)(*spl))->size;
}

/* This function places in the x array the x values of a given spline spl*/
void FC_FUNC_(oct_spline_x, OCT_SPLINE_X)
     (const void **spl, double *x)
{
  int size, i;

  size = (int)((gsl_spline *)(*spl))->size;
  for(i=0; i<size; i++)
    x[i] = ((gsl_spline *)(*spl))->x[i];
}

/* This function places in the y array the y values of a given spline spl*/
void FC_FUNC_(oct_spline_y, OCT_SPLINE_Y)
     (const void **spl, double *y)
{
  int size, i;

  size = (int)((gsl_spline *)(*spl))->size;
  for(i=0; i<size; i++)
    y[i] = ((gsl_spline *)(*spl))->y[i];
}

/* Returns the integral of the spline stored in spl, between a and b */
double FC_FUNC_(oct_spline_eval_integ, OCT_SPLINE_EVAL_INTEG)
     (const void **spl, const double *a, const double *b, void **acc)
{
  /* the GSL headers specify double a, double b */
  return gsl_spline_eval_integ((gsl_spline *)(*spl), *a, *b, (gsl_interp_accel *)(* acc));
}

/* Performs the derivative of a spline */
double FC_FUNC_(oct_spline_eval_der, OCT_SPLINE_EVAL_DER)
     (const double *x, const void **spl, void **acc)
{
  /* the GSL headers specify double x */
  return gsl_spline_eval_deriv((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
}

/* Performs the second derivative of a spline */
double FC_FUNC_(oct_spline_eval_der2, OCT_SPLINE_EVAL_DER2)
     (const double *x, const void **spl, void **acc)
{
  /* the GSL headers specify double x */
  return gsl_spline_eval_deriv2((gsl_spline *)(*spl), *x, (gsl_interp_accel *)(*acc));
}

void FC_FUNC_(oct_strerror, OCT_STRERROR)
     (const fint *err, STR_F_TYPE res STR_ARG1)
{
  const char *c;

  c = gsl_strerror(*err);
  TO_F_STR1(c, res);
}

/* Chebyshev Approximations */
/*
void FC_FUNC_(oct_chebyshev_coeffs, OCT_CHEBYSHEV_COEFFS)
     (gsl_complex *coeffs, int *order)
{
  int i;
  double f (double x, void *p){return cos(x);}
  double g (double x, void *p){return -sin(x);}
  gsl_cheb_series *cs = gsl_cheb_alloc (*order);
  gsl_function F;
  F.function = f;
  F.params = 0;
  gsl_cheb_init (cs, &F, -1.0, 1.0);
  for(i=0; i<=12; i++){GSL_SET_REAL(&coeffs[i], (*cs).c[i]);}
  F.function = g;
  F.params = 0;
  gsl_cheb_init (cs, &F, -1.0, 1.0);
  for(i=0; i<=12; i++){GSL_SET_IMAG(&coeffs[i], (*cs).c[i]);}
}
*/

 /* Numerical Derivatives.
    The following is an interface to the GSL support for numerical derivatives,
    in particular the gsl_deriv_central function. It computes a four points
    approximation to the derivative of a function, supplying an error estimation. */

 /* This is a type used to communicate with Fortran; func_nd is the type of the
    interface to a Fortran subroutine that calculates the value of the function.
    The first argument is the function argument, whereas the second is the function
    value. For convenience reasons, it is wrapped around the "param_nd_t" struct. */
typedef void (*func_nd)(double*, double*);
typedef struct{
  func_nd func;
} param_nd_t;

 /* This is the function that is called by the GSL function gsl_deriv_central. It
    receives as first argument the function argument, and as second argument
    a pointer to a params data type (a GSL data type), which in this case should
    be a pointer to a param_nd_t data type, where the address of the Fortran
    subroutine is. */
double function_oct_numerical_derivative (double x, void * params)
{
  double res;
  param_nd_t * p;

  p = (param_nd_t *) params;
  p->func(&x, &res);
  return res;
}
 /* This is the function that should be called by Fortran. The interface is defined
    in loct_math.F90 file. */
void FC_FUNC_(oct_numerical_derivative, OCT_NUMERICAL_DERIVATIVE)
     (const double *x, const double *h, double *result, double *abserr, const func_nd f)
{
  gsl_function F;
  param_nd_t p;

  p.func = f;
  F.function = &function_oct_numerical_derivative;
  F.params = (void *) &p;
  /* the GSL headers specify double x, double h */
  gsl_deriv_central (&F, *x, *h, result, abserr);
  return;
}

max_stars_repo_name: neelravi/octopus

max_issues_repo_name: neelravi/octopus

max_forks_repo_name: neelravi/octopus


text:
$NetBSD: patch-Src___dotblas.c,v 1.1 2015/01/27 05:05:30 dbj Exp $

--- Src/_dotblas.c.orig	2006-07-24 20:11:35.000000000 +0000
+++ Src/_dotblas.c
@@ -12,7 +12,11 @@ static char module_doc[] =
 #include "Python.h"
 #include "libnumarray.h"
 #include "arrayobject.h"
+#ifdef __APPLE__
+#include <Accelerate/Accelerate.h>
+#else
 #include <cblas.h>
+#endif

 #include <stdio.h>


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"

void
cblas_srot (const int N, float *X, const int incX, float *Y, const int incY,
            const float c, const float s)
{
#define BASE float
#include "source_rot.h"
#undef BASE
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
#ifndef LIC_METADATA_DEFINITION
#define LIC_METADATA_DEFINITION

#include <gsl/gsl>

#include "format/ClrMetadataFormat.h"

namespace lic
{

class Assembly;

class Metadata
{
public:
    Metadata(Assembly& assembly, MetadataTable table, size_t rid, gsl::byte* row);
    ~Metadata();

protected:
    Assembly& assembly;
    MetadataTable table;
    size_t rid;
    gsl::byte* row;
};

}

#endif // !LIC_METADATA_DEFINITION

ext: h

max_stars_repo_name: roberthusak/lic

max_issues_repo_name: roberthusak/lic

max_forks_repo_name: roberthusak/lic


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 dir_710838b0_50ee_47d4_bb24_95e5ded5be0b_h
#define dir_710838b0_50ee_47d4_bb24_95e5ded5be0b_h

#include <direct.h>
#include <io.h>
#include <gslib/type.h>

#ifdef _UNICODE

#define _gs_getcwd _wgetcwd
#define _gs_chdir _wchdir
#define _gs_mkdir _wmkdir
#define _gs_rmdir _wrmdir
#define _gs_finddata_t _wfinddata_t
#define _gs_findfirst _wfindfirst
#define _gs_findnext _wfindnext
#define _gs_findclose _findclose

#else

#define _gs_getcwd _getcwd
#define _gs_chdir _chdir
#define _gs_mkdir _mkdir
#define _gs_rmdir _rmdir
#define _gs_finddata_t _finddata_t
#define _gs_findfirst _findfirst
#define _gs_findnext _findnext
#define _gs_findclose _findclose

#endif

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include "globals.h"
#include "peano.h"
#include <gsl/gsl_heapsort.h>

peanoKey Peano_Key ( const double x, const double y, const double z );
static void reorder_particles();

void Print_Int_Bits128 ( const peanoKey val )
{
    for ( int i = 127; i >= 0; i-- ) {

        printf ( "%llu", ( long long ) ( ( val & ( ( peanoKey ) 1 << i ) ) >> i ) );

        if ( i % 3 == 0 && i != 0 ) {
            printf ( "." );
        }
    }
    printf ( "\n" );
    fflush ( stdout );

    return ;
}

void Print_Int_Bits128r ( const peanoKey val )
{
    for ( int i = 127; i >= 0; i-- ) {

        printf ( "%llu", ( long long ) ( ( val & ( ( peanoKey ) 1 << i ) ) >> i ) );

        if ( i % 3 - 2 == 0 && i != 0 ) {
            printf ( "." );
        }
    }
    printf ( "\n" );
    fflush ( stdout );

    return ;
}

int compare_peanoKeys ( const void *a, const void *b )
{
    const peanoKey *x = ( const peanoKey * ) a;
    const peanoKey *y = ( const peanoKey * ) b;

    return ( int ) ( *x > *y ) - ( *x < *y );
}

static peanoKey *Keys = NULL;
static size_t *Idx = NULL;

void Sort_Particles_By_Peano_Key()
{
    if ( Keys == NULL ) {
        Keys = malloc ( Param.Npart * sizeof ( *Keys ) );
    } else {
        memset ( Keys, 0, Param.Npart * sizeof ( *Keys ) );
    }

    if ( Idx == NULL ) {
        Idx = malloc ( Param.Npart * sizeof ( *Idx ) );
    } else {
        memset ( Idx, 0, Param.Npart * sizeof ( *Idx ) );
    }

    #pragma omp parallel for
    for ( int ipart = 0; ipart < Param.Npart; ipart++ ) {

        double px = P[ipart].Pos[0] / Problem.Boxsize[0];
        double py = P[ipart].Pos[1] / Problem.Boxsize[1];
        double pz = P[ipart].Pos[2] / Problem.Boxsize[2];

        P[ipart].Key = Keys[ipart] = Peano_Key ( px, py, pz );
    }

    gsl_heapsort_index ( Idx, Keys, Param.Npart, sizeof ( *Keys ), &compare_peanoKeys );

    reorder_particles();

    return ;
}

static void reorder_particles()
{
    for ( int i = 0; i < Param.Npart; i++ ) {

        if ( Idx[i] == i ) {
            continue;
        }

        int dest = i;

        struct ParticleData Ptmp = P[i];
        struct GasParticleData Sphtmp = SphP[i];

        int src = Idx[i];

        for ( ;; ) {

            P[dest] = P[src];
            SphP[dest] = SphP[src];

            Idx[dest] = dest;

            dest = src;

            src = Idx[dest];

            if ( src == i ) {
                break;
            }
        }

        P[dest] = Ptmp;
        SphP[dest] = Sphtmp;

        Idx[dest] = dest;

    } // for i

    return ;
}

peanoKey Peano_Key ( const double x, const double y, const double z )
{
    Assert ( x >= 0 && x <= 1, "X coordinate of out range [0,1] have %g", x );
    Assert ( y >= 0 && y <= 1, "Y coordinate of out range [0,1] have %g", y );
    Assert ( z >= 0 && z <= 1, "Z coordinate of out range [0,1] have %g", z );

    const uint64_t m = 1UL << 63; // = 2^63;

    uint64_t X[3] = { y * m, z * m, x * m };

    /* Inverse undo */

    for ( uint64_t q = m; q > 1; q >>= 1 ) {

        uint64_t P = q - 1;

        if ( X[0] & q ) {
            X[0] ^= P;    // invert
        }

        for ( int i = 1; i < 3; i++ ) {

            if ( X[i] & q ) {

                X[0] ^= P; // invert

            } else {

                uint64_t t = ( X[0] ^ X[i] ) & P;

                X[0] ^= t;
                X[i] ^= t;

            } // exchange
        }
    }

    /* Gray encode (inverse of decode) */

    for ( int i = 1; i < 3; i++ ) {
        X[i] ^= X[i - 1];
    }

    uint64_t t = X[2];

    for ( int i = 1; i < 64; i <<= 1 ) {
        X[2] ^= X[2] >> i;
    }

    t ^= X[2];

    for ( int i = 1; i >= 0; i-- ) {
        X[i] ^= t;
    }

    /* branch free bit interleave of transpose array X into key */

    peanoKey key = 0;

    X[1] >>= 1;
    X[2] >>= 2;    // lowest bits not important

    for ( int i = 0; i < N_PEANO_TRIPLETS + 1; i++ ) {

        uint64_t col = ( ( X[0] & 0x8000000000000000 )
                         | ( X[1] & 0x4000000000000000 )
                         | ( X[2] & 0x2000000000000000 ) ) >> 61;

        key <<= 3;

        X[0] <<= 1;
        X[1] <<= 1;
        X[2] <<= 1;

        key |= col;
    }

    key <<= 2;

    return key;
}

/* This constructs the peano key with reversed triplet order. The order in the
 * triplets however is the same ! Also level zero is carried explicitely
 * to ease tree construction. */

peanoKey Reversed_Peano_Key ( const double x, const double y, const double z )
{
    Assert ( x >= 0 && x <= 1, "X coordinate of out range [0,1] have %g", x );
    Assert ( y >= 0 && y <= 1, "Y coordinate of out range [0,1] have %g", y );
    Assert ( z >= 0 && z <= 1, "Z coordinate of out range [0,1] have %g", z );

    const uint64_t m = 1UL << 63; // = 2^63;

    uint64_t X[3] = { y * m, z * m, x * m };

    /* Inverse undo */

    for ( uint64_t q = m; q > 1; q >>= 1 ) {

        uint64_t P = q - 1;

        if ( X[0] & q ) {
            X[0] ^= P;    // invert
        }

        for ( int i = 1; i < 3; i++ ) {

            if ( X[i] & q ) {

                X[0] ^= P; // invert

            } else {

                uint64_t t = ( X[0] ^ X[i] ) & P;

                X[0] ^= t;
                X[i] ^= t;

            } // exchange
        }
    }

    /* Gray encode (inverse of decode) */

    for ( int i = 1; i < 3; i++ ) {
        X[i] ^= X[i - 1];
    }

    uint64_t t = X[2];

    for ( int i = 1; i < 64; i <<= 1 ) {
        X[2] ^= X[2] >> i;
    }

    t ^= X[2];

    for ( int i = 1; i >= 0; i-- ) {
        X[i] ^= t;
    }

    /* branch free reversed (!) bit interleave of transpose array X into key */

    peanoKey key = 0;

    X[0] >>= 18;
    X[1] >>= 19;
    X[2] >>= 20;    // lowest bits not important

    for ( int i = 0; i < N_PEANO_TRIPLETS + 1; i++ ) {

        uint64_t col = ( ( X[0] & 0x4 ) | ( X[1] & 0x2 ) | ( X[2] & 0x1 ) );

        key <<= 3;

        key |= col;

        X[0] >>= 1;
        X[1] >>= 1;
        X[2] >>= 1;
    }

    key <<= 3; // include level 0

    return key;
}

void test_peanokey()
{
    const double box[3]  = { 1.0, 1, 1};
    double a[3] = { 0 };
    int order = 1;
    float delta = 1 / pow ( 2.0, order );
    int n = roundf ( 1 / delta );

    for ( int i = 0; i < n; i++ ) {
        for ( int j = 0; j < n; j++ ) {
            for ( int k = 0; k < n; k++ ) {

                a[0] = ( i + 0.5 ) * delta / box[0];
                a[1] = ( j + 0.5 ) * delta / box[1];
                a[2] = ( k + 0.5 ) * delta / box[2];

                peanoKey stdkey =  Peano_Key ( a[0], a[1], a[2] );

                printf ( "%g %g %g %llu  \n", a[0], a[1], a[2],
                         ( long long ) stdkey );

                Print_Int_Bits128 ( stdkey );

                printf ( "\n" );
            }
        }
    }

    return ;
}

max_stars_repo_name: elehcim/WVTICs

max_issues_repo_name: elehcim/WVTICs

max_forks_repo_name: elehcim/WVTICs


text:
/* diff/diff.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 David Morrison
 *
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_diff.h>

int
gsl_diff_backward (const gsl_function * f,
                   double x, double *result, double *abserr)
{
  /* Construct a divided difference table with a fairly large step
     size to get a very rough estimate of f''.  Use this to estimate
     the step size which will minimize the error in calculating f'. */

  int i, k;
  double h = GSL_SQRT_DBL_EPSILON;
  double a[3], d[3], a2;

  /* Algorithm based on description on pg. 204 of Conte and de Boor
     (CdB) - coefficients of Newton form of polynomial of degree 2. */

  for (i = 0; i < 3; i++)
    {
      a[i] = x + (i - 2.0) * h;
      d[i] = GSL_FN_EVAL (f, a[i]);
    }

  for (k = 1; k < 4; k++)
    {
      for (i = 0; i < 3 - k; i++)
        {
          d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
        }
    }

  /* Adapt procedure described on pg. 282 of CdB to find best value of
     step size. */

  a2 = fabs (d[0] + d[1] + d[2]);

  if (a2 < 100.0 * GSL_SQRT_DBL_EPSILON)
    {
      a2 = 100.0 * GSL_SQRT_DBL_EPSILON;
    }

  h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));

  if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
    {
      h = 100.0 * GSL_SQRT_DBL_EPSILON;
    }

  *result = (GSL_FN_EVAL (f, x) - GSL_FN_EVAL (f, x - h)) / h;
  *abserr = fabs (10.0 * a2 * h);

  return GSL_SUCCESS;
}

int
gsl_diff_forward (const gsl_function * f,
                  double x, double *result, double *abserr)
{
  /* Construct a divided difference table with a fairly large step
     size to get a very rough estimate of f''.  Use this to estimate
     the step size which will minimize the error in calculating f'. */

  int i, k;
  double h = GSL_SQRT_DBL_EPSILON;
  double a[3], d[3], a2;

  /* Algorithm based on description on pg. 204 of Conte and de Boor
     (CdB) - coefficients of Newton form of polynomial of degree 2. */

  for (i = 0; i < 3; i++)
    {
      a[i] = x + i * h;
      d[i] = GSL_FN_EVAL (f, a[i]);
    }

  for (k = 1; k < 4; k++)
    {
      for (i = 0; i < 3 - k; i++)
        {
          d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
        }
    }

  /* Adapt procedure described on pg. 282 of CdB to find best value of
     step size. */

  a2 = fabs (d[0] + d[1] + d[2]);

  if (a2 < 100.0 * GSL_SQRT_DBL_EPSILON)
    {
      a2 = 100.0 * GSL_SQRT_DBL_EPSILON;
    }

  h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));

  if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
    {
      h = 100.0 * GSL_SQRT_DBL_EPSILON;
    }

  *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h;
  *abserr = fabs (10.0 * a2 * h);

  return GSL_SUCCESS;
}

int
gsl_diff_central (const gsl_function * f,
                  double x, double *result, double *abserr)
{
  /* Construct a divided difference table with a fairly large step
     size to get a very rough estimate of f'''.  Use this to estimate
     the step size which will minimize the error in calculating f'. */

  int i, k;
  double h = GSL_SQRT_DBL_EPSILON;
  double a[4], d[4], a3;

  /* Algorithm based on description on pg. 204 of Conte and de Boor
     (CdB) - coefficients of Newton form of polynomial of degree 3. */

  for (i = 0; i < 4; i++)
    {
      a[i] = x + (i - 2.0) * h;
      d[i] = GSL_FN_EVAL (f, a[i]);
    }

  for (k = 1; k < 5; k++)
    {
      for (i = 0; i < 4 - k; i++)
        {
          d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
        }
    }

  /* Adapt procedure described on pg. 282 of CdB to find best
     value of step size. */

  a3 = fabs (d[0] + d[1] + d[2] + d[3]);

  if (a3 < 100.0 * GSL_SQRT_DBL_EPSILON)
    {
      a3 = 100.0 * GSL_SQRT_DBL_EPSILON;
    }

  h = pow (GSL_SQRT_DBL_EPSILON / (2.0 * a3), 1.0 / 3.0);

  if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
    {
      h = 100.0 * GSL_SQRT_DBL_EPSILON;
    }

  *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x - h)) / (2.0 * h);
  *abserr = fabs (100.0 * a3 * h * h);

  return GSL_SUCCESS;
}


text:
#if !defined(lowMachFlow_h)
#define lowMachFlow_h
#include <petsc.h>

// Define the test field number
typedef enum { VTEST, QTEST, WTEST } LowMachFlowTestFields;

typedef enum { STROUHAL, REYNOLDS, FROUDE, PECLET, HEATRELEASE, GAMMA, PTH, MU, K, CP, BETA, GRAVITY_DIRECTION, TOTAL_LOW_MACH_FLOW_PARAMETERS } LowMachFlowParametersTypes;
PETSC_EXTERN const char* lowMachFlowParametersTypeNames[TOTAL_LOW_MACH_FLOW_PARAMETERS + 1];

typedef enum {VEL, PRES, TEMP, TOTAL_LOW_MACH_FLOW_FIELDS} LowMachFlowFields;

typedef enum {MOM, MASS, ENERGY, TOTAL_LOW_MACH_SOURCE_FIELDS} LowMachSourceFields;

// Low Mach Kernels
PETSC_EXTERN void LowMachFlow_vIntegrandTestFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[],
    PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]);

PETSC_EXTERN void LowMachFlow_vIntegrandTestGradientFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                                             const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                                             PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]);

PETSC_EXTERN void LowMachFlow_wIntegrandTestFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[],
                                                     PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]);

PETSC_EXTERN void LowMachFlow_wIntegrandTestGradientFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                                             const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                                             PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]);

PETSC_EXTERN void LowMachFlow_qIntegrandTestFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[],
                                                     PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]);

PETSC_EXTERN void LowMachFlow_g0_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                  PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g1_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g0_qT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g0_qT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g1_qT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g0_vT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g0_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g1_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g3_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g2_vp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g0_wu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g0_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g1_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

PETSC_EXTERN void LowMachFlow_g3_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                                    PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);

#endif

max_stars_repo_name: mtmcgurn-buffalo/ablate

max_issues_repo_name: mtmcgurn-buffalo/ablate

max_forks_repo_name: mtmcgurn-buffalo/ablate


text:
/*
 * Test the effective difference of a single bit changing on the input.
 * Evaluate each bit position on the first hash against all the bit
 * positions on the other hash for any less than random distribution.
 *
 * This tests for the "Bit Independence Criterion" of this hash.
 */

#include "hash_bbt.h"

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <stdint.h>
#include <stdio.h>
#include <string.h>

#define HASH_PARAMS bbt_table_1
extern struct bbt_hash_patterns HASH_PARAMS;

#define HASH_INPUT_SZ 32
//#define SAMPLE_SZ 1048576
//#define SAMPLE_SZ 100
#define SAMPLE_SZ 65536

#define HASH_BITS_COUNT (8*sizeof(bbt_hash_t))
#define COUNT_MATRIX_SZ (HASH_BITS_COUNT * HASH_BITS_COUNT)
#define BIT_ROTATE(HASH,SHIFT) ((bbt_hash_t)((HASH)<<(SHIFT)) | (bbt_hash_t)((HASH)>>(BBT_HASH_WIDTH - (SHIFT))))

int main(int argc, char **argv) {

	if (argc < 2) {
		printf("Usage: %d\n", argc);
		printf("Usage: %s <bit position> <numeric random seed>\n", argv[0]);
		return 1;
	}

	unsigned argPos = 1;
	unsigned long int bitPos = strtoul(argv[argPos++], NULL, 10);
	if (bitPos > ((8*HASH_INPUT_SZ)-1)) {
		fprintf(stderr, "ERROR: bitPos too large. Beyond end of buffer length.\n");
		return 1;
	}
	unsigned long int seed = strtoul(argv[argPos++], NULL, 10);

	unsigned char hashInput[HASH_INPUT_SZ];
	gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937); /* Mersenne Twister */
	gsl_rng_set(rng, seed);

	bbt_hash_ctxt hc;
	bbt_hash_init(&hc, &(HASH_PARAMS));

	unsigned bitIndex = bitPos/8;
	unsigned bitShift = 7 - (bitPos%8);
	bbt_hash_t bitMask = ((bbt_hash_t)1) << bitShift;

	bbt_hash_t bitMatrix[HASH_BITS_COUNT];
	unsigned countMatrix[COUNT_MATRIX_SZ];
	memset(countMatrix, 0, COUNT_MATRIX_SZ * sizeof(unsigned));

	for (unsigned sampleCount = 0; sampleCount < SAMPLE_SZ; sampleCount++) {
		unsigned char *p = hashInput;
		for (unsigned i = 0; i < HASH_INPUT_SZ; i++) {
			*p = gsl_rng_get(rng) & 0xFF;
			p++;
		}

		bbt_hash_calc(&hc, hashInput, HASH_INPUT_SZ);
		bbt_hash_t hashA = bbt_hash_getHash(&hc);
		bbt_hash_reset(&hc);

		hashInput[bitIndex] = hashInput[bitIndex] ^ bitMask;

		bbt_hash_calc(&hc, hashInput, HASH_INPUT_SZ);
		bbt_hash_t hashB = bbt_hash_getHash(&hc);
		bbt_hash_reset(&hc);

		for (unsigned i = 0; i < HASH_BITS_COUNT; i++) {
			bitMatrix[i] = hashB ^ BIT_ROTATE(hashA, i);
		}

		for (unsigned i = 0; i < HASH_BITS_COUNT; i++) {
			bbt_hash_t hashC = bitMatrix[i];
			for (unsigned j = 0; j < HASH_BITS_COUNT; j++) {
				if ((hashC & 1) != 0) {
					countMatrix[j * HASH_BITS_COUNT + i] ++;
				}
				hashC = hashC >> 1;
			}
		}
	}

	printf("cat(\"Bit change at index %u and shift %u.\\n\")\n", bitIndex, bitShift);
	printf("sampleSize <- %u\n", SAMPLE_SZ);
	printf("cat(\"Sample size is %u.\\n\")\n", SAMPLE_SZ);
	printf("cat(\"Half of sample size is %u.\\n\")\n", SAMPLE_SZ/2);
	printf("data <- c(\n");
	for (unsigned i = 0; i < HASH_BITS_COUNT; i++) {
		for (unsigned j = 0; j < HASH_BITS_COUNT; j++) {
			printf("%u, ", countMatrix[j * HASH_BITS_COUNT + i]);
		}
		printf("\n");
	}
	printf("0)\n");
	printf("dataMatrix <- matrix(data[1:%lu], ncol=%lu, byrow=T)\n", COUNT_MATRIX_SZ, HASH_BITS_COUNT);
	printf("cat(\"Summary of the distribution of counts:\\n\")\n");
	printf("print(summary(as.vector(dataMatrix)))\n");
	printf("p <- dataMatrix / sampleSize\n");
	printf("q <- 1-p\n");
	printf("h <- -p*log(p, base=2) -q*log(q, base=2)\n");
	printf("cat(\"Summary of the distribution of H (entropy):\\n\")\n");
	printf("print(summary(as.vector(h)))\n");

	gsl_rng_free(rng);
	return 0;
}

max_stars_repo_name: jonkanderson/BitBalancedTableHash

max_issues_repo_name: jonkanderson/BitBalancedTableHash

max_forks_repo_name: jonkanderson/BitBalancedTableHash


text:
/*
  Copyright [2017-2020] [IBM Corporation]
  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 __BUFFER_MGR_H__
#define __BUFFER_MGR_H__

#ifdef __cplusplus

#include <common/chksum.h>
#include <common/delete_copy.h>
#include <common/logging.h>
#include <common/moveable_ptr.h>
#include <common/utils.h> /* MiB, UNLIKELY */
#include <api/components.h>
#include <api/fabric_itf.h>
#include <api/kvstore_itf.h>
#include <api/registrar_memory_direct.h> /* Opaque_memory_region */
#include <sys/mman.h>

#include "mcas_config.h"
#include "memory_registered.h"
#include <gsl/pointers> /* not_null */
#include <cassert>
#include <cstring>      /* memset */
#include <memory>       /* unique_ptr */

namespace
{
  inline auto alloc_base(std::size_t len) -> gsl::not_null<void *>
  {
    auto b = ::aligned_alloc(MiB(2), len);
    if (b == nullptr) {
      throw std::bad_alloc();
    }
    ::memset(b, 0, len);
    return gsl::not_null<void *>(b);
  }
}

namespace mcas
{
template <class Memory>
class Buffer_manager : private common::log_source {

  static constexpr const char *_cname       = "Buffer_manager";

public:
  static constexpr size_t DEFAULT_BUFFER_COUNT = NUM_SHARD_BUFFERS;
  static constexpr size_t BUFFER_LEN           = MiB(2); /* corresponds to huge page see below */
  using memory_registered_t                    = memory_registered<Memory>;

  using memory_region_t = typename Memory::memory_region_t;

  struct iov_mem_lock
  {
  private:
    common::moveable_ptr<::iovec> _iov;
  public:
    iov_mem_lock(::iovec *iov_)
      : _iov(iov_)
    {
      ::madvise(_iov->iov_base, _iov->iov_len, MADV_HUGEPAGE);
      ::mlock(_iov->iov_base, _iov->iov_len);
    }
    iov_mem_lock(iov_mem_lock &&) noexcept = default;
    ~iov_mem_lock()
    {
      if ( _iov )
      {
        ::munlock(_iov->iov_base, _iov->iov_len);
      }
    }
  };

  /* Although client tried to hide it with a reinterpret_cast, buffer_base is
   used as a component::memory_region_t */

  struct buffer_base : public component::Registrar_memory_direct::Opaque_memory_region, private common::log_source
  {
  private:
    static constexpr const char *_cname = "buffer_base";
    void *_base;
    std::size_t _original_length;
    memory_registered_t _region;
    const unsigned magic;

  public:
    buffer_base(unsigned    debug_level_,
             Memory * transport_,
             void *      base_,
             size_t      length_
    )
      : common::log_source(debug_level_)
      , _base(base_)
      , _original_length(length_)
      , _region(memory_registered_t(debug_level_, transport_, base_, length_, 0, 0))
      , magic(0xC0FFEE)
    {
      if ( length_ <= 1 )
      {
        throw std::domain_error("buffer length too small");
      }
      CPLOG(3, "%s::%s %p transport %p region %p", _cname, __func__, common::p_fmt(this),
            common::p_fmt(transport()), common::p_fmt(region()));
    }

    void *get_desc() const { return _region.desc(); }

    DELETE_COPY(buffer_base);

  protected:
    ~buffer_base()
    {
      CPLOG(3, "%s::%s %p transport %p region %p", _cname, __func__, common::p_fmt(this),
            common::p_fmt(transport()), common::p_fmt(region()));
    }

  public:
    size_t original_length() const { return _original_length; }
    inline gsl::not_null<void *> base() const { return _base; }
    auto region() const { return _region.mr(); }
    Memory * transport() const { return _region.transport(); }

#if 0
    /* unused */
    inline bool check_magic() const { return magic == 0xC0FFEE; }
    unsigned int crc32() const { return common::chksum32(_base, _length); }
#endif
  };

  struct free_deleter
  {
    void operator()(void *v_) { ::free(v_); }
  };

  struct buffer_internal
    : private std::unique_ptr<void, free_deleter>
    , public buffer_base
  {
    ::iovec iov[2];
    void *desc[2];
    iov_mem_lock _ml;
    using completion_t = void (*)(void *, buffer_internal *);
    completion_t   completion_cb;
    void *         value_adjunct;
    inline void set_completion(completion_t completion_) { completion_cb = completion_; }
    void set_completion(completion_t completion_, void *value_adjunct_)
    {
      completion_cb = completion_;
      value_adjunct = value_adjunct_;
    }
    buffer_internal(unsigned debug_level_,
             Memory *         transport_,
             size_t              length_
    )
      : std::unique_ptr<void, free_deleter>(alloc_base(length_))
      , buffer_base(debug_level_, transport_, std::unique_ptr<void, free_deleter>::get(), length_)
      , iov{{this->base(), this->original_length()}, {nullptr, 0}}
      , desc{this->get_desc(), nullptr}
      , _ml(&this->iov[0])
      , completion_cb(nullptr)
      , value_adjunct(nullptr)
    {
    }

    DELETE_COPY(buffer_internal);

    void set_length(size_t s) { iov->iov_len = s; }

    inline void reset_length()
    {
      assert(this->original_length() > 1);
      this->iov[0].iov_len  = this->original_length();
      value_adjunct = nullptr;
      completion_cb = nullptr;
    }

    size_t length() const { return iov[0].iov_len; }

    unsigned int crc32() const { return common::chksum32(iov->iov_base, iov->iov_len); }
  };

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"  // missing initializers
  Buffer_manager(unsigned debug_level_,
                 Memory *transport,
                 size_t buffer_count = DEFAULT_BUFFER_COUNT)
    : common::log_source(debug_level_),
      _buffer_count(buffer_count),
      _transport(transport)
  {
    for (unsigned i = 0; i < _buffer_count; i++) {
      const auto len  = BUFFER_LEN;
      _buffers.emplace_back(std::make_unique<buffer_internal>(debug_level_, _transport, len));
      _free.push_back(_buffers.back().get());
    }
    PLOG("%s %p allocated %lu buffers", __func__, common::p_fmt(this), buffer_count);
  }
#pragma GCC diagnostic pop

  Buffer_manager(Buffer_manager &&) noexcept = default;

  ~Buffer_manager()
  {
  }

  using completion_t = void (*)(void *, buffer_internal *);

  gsl::not_null<buffer_internal *> allocate(completion_t completion_)
  {
    if (UNLIKELY(_free.empty())) throw Program_exception("Buffer_manager: no shard buffers remaining");
    gsl::not_null<buffer_internal *> iob = _free.back();
    _free.pop_back();
    CPLOG(3, "%s::%s %p (%lu free)", _cname, __func__, common::p_fmt(iob), _free.size());
    iob->reset_length();
    iob->set_completion(completion_);
    return iob;
  }

  void free(gsl::not_null<buffer_internal *> iob)
  {
    CPLOG(3, "%s::%s %p (%lu free)", _cname, __func__, common::p_fmt(iob), _free.size());
    iob->reset_length();
    iob->set_completion(nullptr);
    _free.push_back(iob);
  }

private:
  // static constexpr size_t buffer_len() { return BUFFER_LEN; }

  using pool_t = component::IKVStore::pool_t;
  using key_t  = std::uint64_t;

  const size_t                           _buffer_count;
  gsl::not_null<Memory *>             _transport;
  std::vector<std::unique_ptr<buffer_internal>> _buffers;
  std::vector<buffer_internal *>        _free;
};
}  // namespace mcas

#endif
#endif

max_stars_repo_name: fQuinzan/mcas

max_issues_repo_name: fQuinzan/mcas

max_forks_repo_name: fQuinzan/mcas


text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		       High Level Game Framework            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <random>
#include <halley/utils/utils.h>
#include <gsl/span>

namespace Halley {
	typedef std::mt19937 RandomGenType;

	class Random {
	public:
		static Random& getGlobal();

		Random();
		Random(long seed);
		Random(char* bytes, size_t nBytes);

		template <typename T> T get(T max) { return get(T(0), max); }
		template <typename T> T get(T min, T max);

		template <typename T>
		T getInt(T min, T max)
		{
			if (min > max) {
				std::swap(min, max);
			}
			std::uniform_int_distribution<T> dist(min, max);
			return dist(generator);
		}

		template <typename T>
		T getFloat(T min, T max)
		{
			if (min > max) std::swap(min, max);
			std::uniform_real_distribution<T> dist(min, max);
			return dist(generator);
		}

		void getBytes(gsl::span<gsl::byte> dst);
		void setSeed(long seed);
		void setSeed(char* bytes, size_t nBytes);

	private:
		RandomGenType generator;
	};

	template <typename T>
	struct _getHack {
		static T get(Random* rng, T min, T max) { return rng->getInt(min, max); }
	};

	template <>
	struct _getHack<float> {
		static float get(Random *rng, float min, float max) { return rng->getFloat(min, max); }
	};

	template <>
	struct _getHack<double> {
		static double get(Random *rng, double min, double max) { return rng->getFloat(min, max); }
	};

	template <typename T>
	T Random::get(T min, T max)
	{
		return _getHack<T>::get(this, min, max);
	}

}

max_stars_repo_name: Healthire/halley

max_issues_repo_name: Healthire/halley

max_forks_repo_name: Healthire/halley


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2017 Couchbase, Inc
 *
 *   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.
 */

#pragma once

#include <gsl_p/span.h>

namespace phosphor {

/**
 * Pure virtual base class for receiving stats from phosphor internals
 *
 * Methods on this class will be called with the key and value of each stat.
 * Callback implementations MUST NOT re-enter the TraceLog or TraceBuffer
 * upon which it used as locks may be held when the callback is invoked.
 *
 * Example usage:
 *
 *     class MyStatsCallback : public phosphor::StatsCallback {
 *          // Implement callback methods
 *     } callback;
 *
 *     phosphor::TraceLog::getInstance().getStats(callback);
 *
 *     // Stash data as required for application
 *     auto data = callback.getData();
 *
 * Implementations should note that phosphor makes no guarantees about
 * atomicity of the stats with respect to each other.
 */
class StatsCallback {
public:
    /**
     * Utility template method for passing a string literal as the first
     * argument instead of a string span
     * @param s String to convert to a span
     * @param value Value to forward
     */
    template <size_t N, typename T>
    void operator()(const char (&s)[N], T&& value) {
        this->operator()(gsl_p::make_span(s), std::forward<T>(value));
    }

    virtual void operator()(gsl_p::cstring_span key,
                            gsl_p::cstring_span value) = 0;
    virtual void operator()(gsl_p::cstring_span key, bool value) = 0;
    virtual void operator()(gsl_p::cstring_span key, size_t value) = 0;
    virtual void operator()(gsl_p::cstring_span key, ssize_t value) = 0;
    virtual void operator()(gsl_p::cstring_span key, double value) = 0;
};

}

max_stars_repo_name: Chippiewill/Phosphor

max_issues_repo_name: Chippiewill/Phosphor

max_forks_repo_name: Chippiewill/Phosphor