text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mersenne.h"
#include <complex.h>
#include <lapacke.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327
#endif

/* Lattice size and dimensions */
#define NT 48
#define NX 6

#define ND 2
#define NDIRS (2*ND)

#define TUP 0
#define XUP 1
#define TDN 2
#define XDN 3

#define VOLUME (NT*NX)

//One staggered fermion
#define Nf 2

//#define ANTISYMMETRIC //Antisymmetric boundaries
//#define SYMMETRIC     //implemented in Thirring_hop
#define OPENX       //open in space, (anti)symmetric in time

#define CG_ACCURACY 1e-30
#define CG_MAX_ITER 10000

/* Simulation parameters */
double m;
double g;
double mu;

/* Neighbour index arrays, to be filled at the beginning
 */
int *tup,*xup,*tdn,*xdn;

double ***A;
int    ***eta;   //Staggered eta matrix

/* Functions for fetching neighboring coordinates */
static inline int tdir(int t, int dir){
  if( dir == TUP ) return tup[t];
  if( dir == TDN ) return tdn[t];
  return(t);
}

static inline int xdir(int x, int dir){
  if( dir == XUP ) return xup[x];
  if( dir == XDN ) return xdn[x];
  return(x);
}

/* Opposite of a direction */
static inline int opp_dir(int dir){
  return ( dir + ND ) % NDIRS;
}

double gauge_action(){
    double S = 0;
    for( int t=0;t<NT;t++) for( int x=0;x<NX;x++) for( int dir=0; dir<ND; dir++ ) {
        S+= 1.-cos(A[t][x][dir]);
    }
    return Nf/g*S;
}

void update_gauge(){
    for( int t=0;t<NT;t++) for( int x=0;x<NX;x++) for( int dir=0; dir<ND; dir++ ) {
        double s1 = -Nf/g*cos(A[t][x][dir]);
        double new = 2*M_PI*mersenne()-M_PI;
        double s2 = -Nf/g*cos(new);
        double edS = exp(-s2+s1);
        //printf("New A=%f at (%d,%d), direction %d, exp(-deltaS)=%g\n", A[t][x][dir],t,x,dir, edS);

        if( mersenne() < edS ) {
            A[t][x][dir] = new;
            //printf("Accepted\n");
        }
    }
}

void fermion_matrix( _Complex double * M ){
  static double expmu, expmmu;
  static int init=1;
  if(init){
    expmu = exp(mu);
    expmmu = exp(-mu);
    init = 0;
  }

  for( int i=0;i<VOLUME*VOLUME;i++)  M[i] = 0;

  for( int t=0;t<NT;t++) for( int x=0;x<NX;x++){
    M[(NX*t + x) + (NX*t+x)*VOLUME] = m;
    int t2 = tup[t];
    int index = (NX*t + x) + (NX*t2+x)*VOLUME;
    double sA = sin(A[t][x][0]);
    double cA = cos(A[t][x][0]);
    if( t2 > t ) M[index] = 0.5 * (cA + sA*I) * eta[t][x][0] * expmu ;
    else M[index] = -0.5 * (cA + sA*I) * eta[t][x][0] * expmu ;

    t2 = tdn[t];
    sA = sin(A[t2][x][0]);
    cA = cos(A[t2][x][0]);
    index = (NX*t + x) + (NX*t2+x)*VOLUME;
    if( t2 > t ) M[index] = 0.5 * (cA - sA*I) * eta[t2][x][0] * expmmu ;
    else M[index] = -0.5 * (cA - sA*I) * eta[t2][x][0] * expmmu ;

    int x2 = xup[x];
    index = (NX*t + x) + (NX*t+x2)*VOLUME;
    sA = sin(A[t][x][1]);
    cA = cos(A[t][x][1]);
    if( x2 > x ) M[index] = 0.5 * (cA + sA*I) * eta[t][x][1] ;
    else M[index] = -0.5 * (cA + sA*I) * eta[t][x][1] ;

    x2 = xdn[x];
    index = (NX*t + x) + (NX*t+x2)*VOLUME;
    sA = sin(A[t][x2][1]);
    cA = cos(A[t][x2][1]);
    if( x2 > x ) M[index] = 0.5 * (cA - sA*I) * eta[t][x2][1] ;
    else M[index] = -0.5 * (cA - sA*I) * eta[t][x2][1] ;
  }
}

double complex determinant(){
    double complex * M; //The matrix
    M = (_Complex double *)malloc( VOLUME*VOLUME*sizeof(_Complex double) );
    fermion_matrix(M);
    int n=VOLUME;
    int info;
    int *ipiv;
    ipiv = malloc( n*sizeof(int) );
    double complex det = 1;

    /*for( int i=0;i<VOLUME;i++) { for( int j=0;j<VOLUME;j++){
        printf(" %4.2f ", creal( M[i+j*VOLUME] ) );
      }
      printf("\n");
    }
    printf("\n");
    */

    LAPACK_zgetrf( &n, &n, M, &n, ipiv, &info );

    if( info > 0 ) printf("error %d\n", info);

    for(int i=0; i<n; i++) {
      det *= M[i*n+i];
      if( ipiv[i] != i+1 ) det*=-1;
    }

    /*
    int lwork=n*n;
    _Complex double *work;
    work = malloc( lwork*sizeof(_Complex double) );

    LAPACK_zgetri(&n, M, &n, ipiv, work, &lwork, &info);
        free(work);

    for( int i=0;i<VOLUME;i++) { for( int j=0;j<VOLUME;j++){
        printf(" %4.2f ", creal( M[i+j*VOLUME] ) );
      }
      printf("\n");
    }
    printf("\n");
    */

    free(M);
    free(ipiv);
    return( det );
}

int n_average;
void measure(){
  /* Measurements */
  static int n=0;
  static double M = 0;
  static double complex det = 0;
  static double sign = 1;

  double complex this_det = determinant();
  det += this_det;
  double this_sign = (creal(this_det) > 0) ? 1 : -1;
  sign += this_sign;

  for( int t=0;t<NT;t++) for( int x=0;x<NX;x++) for( int dir=0; dir<ND; dir++) {
    M += this_sign*A[t][x][dir];
  }

  n++;

  if(n%n_average == 0){
     printf("Magnetisation %g\n",M/(n_average*VOLUME));
     printf("Determinant Re %g Im %g\n",creal(det)/(n_average),cimag(det)/(n_average));
     printf("Sign %g\n",sign/(n_average));
     sign = 0; det = 0; M = 0;
  }

}

int main(void){

  #ifdef DEBUG
  feenableexcept(FE_INVALID | FE_OVERFLOW);
  #endif

  long seed;
  int n_loops, n_measure;

  /* Read in the input */
  printf(" Number of updates : ");
  scanf("%d",&n_loops);

  printf(" Updates / measurement : ");
  scanf("%d",&n_measure);

  printf(" Average over configurations : ");
  scanf("%d",&n_average);

  printf(" m : ");
  scanf("%lf",&m);

  printf(" g : ");
  scanf("%lf",&g);

  printf(" mu : ");
  scanf("%lf",&mu);

  printf(" Random number : ");
  scanf("%ld",&seed);
  seed_mersenne( seed );

  /* "Warm up" the rng generator */
  for ( int i=0; i<543210; i++) mersenne();

  printf(" \n++++++++++++++++++++++++++++++++++++++++++\n");
  printf(" 2D quenched Thirring model, ( %d , %d ) lattice\n", NT, NX );
  printf(" %d updates per measurements\n", n_measure );
  printf(" m %f \n", m);
  printf(" g %f \n", g);
  printf(" mu %f \n", mu);
  printf(" Random seed %ld\n", seed );

  eta = malloc( NT*sizeof(int *) );
  A = malloc( NT*sizeof(double *) );
  for (int t=0; t<NT; t++){
    A[t] = malloc( (NX+1)*sizeof(double *) );
    eta[t] = malloc( (NX+1)*sizeof(int *) );
    for (int x=0; x<NX+1; x++) {
      eta[t][x] = malloc( ND*sizeof(int) );
      A[t][x] = malloc( ND*sizeof(double) );
    }
  }
  xup = malloc( (NX+1)*sizeof(int) );
  xdn = malloc( (NX+1)*sizeof(int) );
  tup = malloc( NT*sizeof(int) );
  tdn = malloc( NT*sizeof(int) );

  /* fill up the index array */
  for ( int i=0; i<NT; i++ ) {
    tup[i] = (i+1) % NT;
    tdn[i] = (i-1+NT) % NT;
  }
  for (int i=0; i<NX+1; i++) {
    xup[i] = (i+1) % NX;
    xdn[i] = (i-1+NX) % NX;
  }

  for (int t=0; t<NT; t++) for (int x=0; x<NX; x++) {
    for( int dir=0; dir<ND; dir++ ) A[t][x][dir] = 0;
    eta[t][x][1] = 1;
    if( x%2 == 0 ){
      eta[t][x][0] = 1;
    } else {
      eta[t][x][0] = -1;
    }
#ifdef OPENX
    eta[t][NX][0] = eta[t][NX][1] = 0;
#endif
  }

  for (int i=1; i<n_loops+1; i++) {
    update_gauge();

    //printf("Updated gauge\n");
    //printf("Action %g\n", gauge_action());

    if((i%n_measure)==0){
      /* Statistics */
      measure();
    }
  }

  return 0;
}


text:
/*
    This code derives from the R wavelets package and it is modified by Davide
    Albanese <albanese@fbk.it>.

    The Python interface is written by  Davide Albanese <albanese@fbk.it>.
    (C) 2009 mlpy Developers.

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

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

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

#include <Python.h>
#include <numpy/arrayobject.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_math.h>

#define SQRT_2 1.4142135623730951

void
uwt_forward (double *V, int n, int j,
	      double *h, double *g, int l,
	      double *Wj, double *Vj)
{
  int t, k, z;
  double k_div;

  for(t = 0; t < n; t++)
    {
      k = t;
      Wj[t] = h[0] * V[k];
      Vj[t] = g[0] * V[k];

      for(z = 1; z < l; z++)
	{
	  k -= (int) pow(2, (j - 1));

	  k_div = -k / (double) n;

	  if(k < 0)
	    k += (int) ceil(k_div) * n;

	  Wj[t] += h[z] * V[k];
	  Vj[t] += g[z] * V[k];
	}
    }
}

void
uwt_backward (double *W, double *V, int j,
	       int n, double *h, double *g,
	       int l, double *Vj)
{
  int t, k, z;
  double k_div;

  for(t = 0; t < n; t++)
    {
      k = t;
      Vj[t] = h[0] * W[k] + g[0] * V[k];

      for(z = 1; z < l; z++)
	{
	  k += (int) pow(2, (j - 1));
	  k_div = (double) k / (double) n;

	  if(k >= n)
	    k -= (int) floor(k_div) * n;

	  Vj[t] += h[z] * W[k] + g[z] * V[k];
	}
    }
}

static PyObject *uwt_uwt(PyObject *self, PyObject *args, PyObject *keywds)
{
  PyObject *x = NULL; PyObject *xa = NULL;
  PyObject *Xa = NULL;
  int levels = 0;
  char wf;

  int i, k, j, n, J;

  double *_x;
  double *_X;

  double *v;
  double *wj, *vj;
  double *h, *g;

  npy_intp Xa_dims[2];

  gsl_wavelet *wave;

  /* Parse Tuple*/
  static char *kwlist[] = {"x", "wf", "k", "levels", NULL};
  if (!PyArg_ParseTupleAndKeywords(args, keywds, "Oci|i", kwlist, &x, &wf, &k, &levels))
    return NULL;

  xa = PyArray_FROM_OTF(x, NPY_DOUBLE, NPY_IN_ARRAY);
  if (xa == NULL) return NULL;

  n = (int) PyArray_DIM(xa, 0);
  _x = (double *) PyArray_DATA(xa);

  switch (wf)
    {
    case 'd':
      wave = gsl_wavelet_alloc (gsl_wavelet_daubechies_centered, k);
      break;

    case 'h':
      wave = gsl_wavelet_alloc (gsl_wavelet_haar_centered, k);
      break;

    case 'b':
      wave = gsl_wavelet_alloc (gsl_wavelet_bspline_centered, k);
      break;

    default:
      PyErr_SetString(PyExc_ValueError, "wavelet family is not valid");
      return NULL;
    }

  h = (double *) malloc(wave->nc * sizeof(double));
  g = (double *) malloc(wave->nc * sizeof(double));

  for(i=0; i<wave->nc; i++)
    {
      h[i] = wave->h1[i] / SQRT_2;
      g[i] = wave->g1[i] / SQRT_2;
    }

  if (levels == 0)
    J = (int) floor(log(((n-1) / (wave->nc-1)) + 1) / log(2));
  else
    J = levels;

  Xa_dims[0] = (npy_intp) (2 * J);
  Xa_dims[1] = PyArray_DIM(xa, 0);
  Xa = PyArray_SimpleNew(2, Xa_dims, NPY_DOUBLE);
  _X = (double *) PyArray_DATA(Xa);

  v = _x;
  for(j=0; j<J; j++)
    {
      wj = _X +(j * n);
      vj = _X +((j + J) * n);
      uwt_forward(v, n, j+1, g, h, wave->nc, wj,  vj);
      v = vj;
    }

  gsl_wavelet_free(wave);
  free(h);
  free(g);
  Py_DECREF(xa);

  return Py_BuildValue("N", Xa);
}

static PyObject *uwt_iuwt(PyObject *self, PyObject *args, PyObject *keywds)
{
  PyObject *X = NULL; PyObject *Xa = NULL;

  PyObject *xa = NULL;

  char wf;
  int i, k, n, J;
  double *w1, *v1;

  double *_X;
  double *_x;
  double *h, *g;

  npy_intp xa_dims[1];

  gsl_wavelet *wave;

  /* Parse Tuple*/
  static char *kwlist[] = {"X", "wf", "k", NULL};
  if (!PyArg_ParseTupleAndKeywords(args, keywds, "Oci", kwlist, &X, &wf, &k))
    return NULL;

  Xa = PyArray_FROM_OTF(X, NPY_DOUBLE, NPY_IN_ARRAY);
  if (Xa == NULL) return NULL;

  n = (int) PyArray_DIM(Xa, 1);
  J = ((int) PyArray_DIM(Xa, 0)) / 2;

  _X = (double *) PyArray_DATA(Xa);

  switch (wf)
    {
    case 'd':
      wave = gsl_wavelet_alloc (gsl_wavelet_daubechies, k);
      break;

    case 'h':
      wave = gsl_wavelet_alloc (gsl_wavelet_haar, k);
      break;

    case 'b':
      wave = gsl_wavelet_alloc (gsl_wavelet_bspline, k);
      break;

    default:
      PyErr_SetString(PyExc_ValueError, "wavelet family is not valid");
      return NULL;
    }

  h = (double *) malloc(wave->nc * sizeof(double));
  g = (double *) malloc(wave->nc * sizeof(double));

  for(i=0; i<wave->nc; i++)
    {
      h[i] = wave->h2[i] / SQRT_2;
      g[i] = wave->g2[i] / SQRT_2;
    }

  w1 = _X;
  v1 = _X + (J * n);

  xa_dims[0] = (npy_intp) n;
  xa = PyArray_SimpleNew(1, xa_dims, NPY_DOUBLE);
  _x = (double *) PyArray_DATA(xa);

  uwt_backward (w1, v1, 1, n, g, h, wave->nc, _x);

  gsl_wavelet_free(wave);
  free(h);
  free(g);
  Py_DECREF(Xa);

  return Py_BuildValue("N", xa);
}

/* Doc strings: */
static char module_doc[]  = "Undecimated Wavelet Transform Module";

static char uwt_uwt_doc[] =
  "Undecimated Wavelet Tranform\n\n"
  ":Parameters:\n"
  "   x : 1d array_like object (the length is restricted to powers of two)\n"
  "      data\n"
  "   wf : string ('d': daubechies, 'h': haar, 'b': bspline)\n"
  "      wavelet family\n"
  "   k : int\n"
  "      member of the wavelet family\n\n"
  "         * daubechies: k = 4, 6, ..., 20 with k even\n"
  "         * haar: the only valid choice of k is k = 2\n"
  "         * bspline: k = 103, 105, 202, 204, 206, 208, 301, 303, 305 307, 309\n\n"
  "   levels : int\n"
  "      level of the decomposition (J).\n"
  "      If levels = 0 this is the value J such that the length of X\n"
  "      is at least as great as the length of the level J wavelet filter,\n"
  "      but less than the length of the level J+1 wavelet filter.\n"
  "      Thus, j <= log_2((n-1)/(l-1)+1), where n is the length of x\n\n"
  ":Returns:\n"
  "   X : 2d numpy array (2J * len(x))\n"
  "      misaligned scaling and wavelet coefficients::\n\n"
  "         [[wavelet coefficients W_1]\n"
  "          [wavelet coefficients W_2]\n"
  "                        :\n"
  "          [wavelet coefficients W_J]\n"
  "          [scaling coefficients V_1]\n"
  "          [scaling coefficients V_2]\n"
  "                       :\n"
  "          [scaling coefficients V_J]]"
  ;

static char uwt_iuwt_doc[] =
  "Inverse Undecimated Wavelet Tranform\n\n"
  ":Parameters:\n"
  "   X : 2d array_like object (the length is restricted to powers of two)\n"
  "      misaligned scaling and wavelet coefficients\n"
  "   wf : string ('d': daubechies, 'h': haar, 'b': bspline)\n"
  "      wavelet family\n"
  "   k : int\n"
  "      member of the wavelet family\n\n"
  "         * daubechies: k = 4, 6, ..., 20 with k even\n"
  "         * haar: the only valid choice of k is k = 2\n"
  "         * bspline: k = 103, 105, 202, 204, 206, 208, 301, 303, 305 307, 309\n\n"
  ":Returns:\n"
  "   x : 1d numpy array\n"
  "      data"
  ;

/* Method table */
static PyMethodDef uwt_methods[] = {
  {"uwt",
   (PyCFunction)uwt_uwt,
   METH_VARARGS | METH_KEYWORDS,
   uwt_uwt_doc},
  {"iuwt",
   (PyCFunction)uwt_iuwt,
   METH_VARARGS | METH_KEYWORDS,
   uwt_iuwt_doc},
  {NULL, NULL, 0, NULL}
};

#if PY_MAJOR_VERSION >= 3

static struct PyModuleDef moduledef = {
  PyModuleDef_HEAD_INIT,
  "_uwt",
  module_doc,
  -1,
  uwt_methods,
  NULL, NULL, NULL, NULL
};

PyObject *PyInit__uwt(void)
{
  PyObject *m;
  m = PyModule_Create(&moduledef);
  if (!m) {
    return NULL;
  }

  import_array();

  return m;
}

#else

PyMODINIT_FUNC init_uwt(void)
{
  PyObject *m;

  m = Py_InitModule3("_uwt", uwt_methods, module_doc);
  if (m == NULL) {
    return;
  }

  import_array();
}

#endif

ext: c


text:
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#ifdef _MSC_VER				/* identifies this as a Microsoft compiler */
#define _USE_MATH_DEFINES	/* added RX2011 */
#endif
#include <math.h>
#include <time.h>
#include <limits.h>
#include <gsl/gsl_multimin.h>
//#include "/Users/tischler/dev/GNU/include/gsl/gsl_multimin.h"
#include "Euler.h"

#ifdef DEBUG_ON
#define DEBUG 1				/* a value of 5 gets everything */
#endif

#define MIN(A,B)	( (A)<(B)?(A):(B) )
#define MAX(A,B)	( (A)>(B)?(A):(B) )

double my_f (const gsl_vector *v, void *params);

/******************************************************************************************
 * return the error between a lattice rotated by Euler angles and the measured spots
 * thus the returned value is a minimum for the optimimal Euler angles
 */
double my_f(
const gsl_vector *v,	/* vector with current value of Euler angles */
void	*pattern)		/* pointer to the parameters, in this case a pointer to the pattern structure */
{
	double  G[3];						/* calculated G vector for each hkl */
	double  recip[3][3];				/* reciprocal space */
	double  M_Euler[3][3];				/* an Euler matrix */
	size_t  N;							/* local copy of number of measured spots to check */
	double a,b,g;						/* local copy of the Euler angles */
	struct patternOfOneGrain *p;		/* local value of pattern */
	double  rms;						/* rms value of error, the returned value */
	size_t  i;							/* loop index */
	double  hkl[3];						/* float value of hkl */
	double  dot;

	/* p = pattern; */
	p = (struct patternOfOneGrain *)pattern;
	N = p->Ni;											/* local value for convenience */
	a = gsl_vector_get(v,0);							/* the Euler angles */
	b = gsl_vector_get(v,1);
	g = gsl_vector_get(v,2);
	MatrixCopy33(recip,p->xtal.recip);					/* copy unrotated reciprocal lattice into recip */

	EulerMatrix(a,b,g,M_Euler);							/* make the rotation matrix M_Euler from Euler angles */
	MatrixMultiply33(M_Euler,recip,recip);				/* rotate recip by Euler angles */
	for (i=0,rms=0.0;i<N;i++) {
		VECTOR_COPY3(hkl,p->hkls[i]);					/* convert integer hkl to float for following MatrixMultiply31 */
if (hkl[0]!=hkl[0]) {
fprintf(stderr,"hkl[0] is NaN\n");
}
		MatrixMultiply31(recip,hkl,G);
if (G[0]!=G[0]) {
fprintf(stderr,"G[0] is NaN\n");
}
		normalize3(G);									/* normalized gvec of hkl[i] rotated by Euler angles */
		dot = dot3(G,p->Ghat[i]);
		dot = MIN(dot,1.0);
		rms += 1. - dot;								/* (1-dot) ~ 0.5*(delta angle)^2*/
	}
	rms = sqrt(rms/(double)N);							/* rms of angle errors */
	return rms;
}

/* When calling optimizeEulerAngles(), be sure to set pattern.xtal:

	e.g.
	pattern.xtal.a = pattern.xtal.b = pattern.xtal.c = 4.05;		// just default to Al for now
	pattern.xtal.alpha = pattern.xtal.beta = pattern.xtal.gamma = M_PI/2.0;

	then call:
		FillRecipInLattice(&(pattern.lattice));
	to set the reciprocal lattice in the lattice structure
*/
/* this routine optimizes approximate Euler angles using the gsl simplex method. */
int optimizeEulerAngles(
double  startStep,					/* starting step size (radians), make this a bit larger than the error */
double  epsAbs,						/* tolerance on the answer (=0.001) */
long	maxIter,					/* maximum number of iterations, stop after this many regardless, perhaps 100 */
struct	patternOfOneGrain *pattern) /* provides G^'s and hkl's for one fitted pattern, and the lattice parameters */
{
	size_t np = 3;					/* number of dimensions of the function, 3 Euler angles */

	/* T tells what kind of minimizer we are using, select the simplex method */
	const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
	gsl_multimin_fminimizer *s = NULL;		/* this is the minimizer */
	gsl_vector *ss, *x;						/* pointers to vectors,  ss is stepsize, x is point */
	gsl_multimin_function minex_func;		/* struct defining parts needed for minimizing */
	size_t iter=0, i;
	int status;								/* status of returned gsl function, 0 is OK */
	double size;							/* size of minimizer, sort of maximum error */

	double  G[3];							/* calculated G vector for each hkl */
	double  recip[3][3];					/* reciprocal space */
	double  M_Euler[3][3];					/* an Euler matrix */
	double  hkl[3];							/* float value of hkl */
	double  dot;

	if (!pattern) { fprintf(stderr,"pattern in NULL on entry to optimizeEulerAngles()\n"); return 1; }
	if (pattern->Ni<3) { fprintf(stderr,"less than 3 pairs of Ghat & hkls on entry to optimizeEulerAngles()\n"); return 1; }
	if (!(pattern->Ghat) || !(pattern->hkls)) { fprintf(stderr,"Ghat or hkls are NULL on entry to optimizeEulerAngles()\n"); return 1; }
	if ((pattern->xtal.a)<=0 || (pattern->xtal.b)<=0 || (pattern->xtal.c)<=0) { fprintf(stderr,"invalid lattice in optimizeEulerAngles()\n"); return 1; }

	/* Initial vertex step size vector */
	ss = gsl_vector_alloc(np);				/* step sizes in x and y */

	/* Set all step sizes to startStep, the starting step size */
	gsl_vector_set_all(ss, startStep);		/* set ss[0] = ss[1]  = ss[2] = startStep (radian) */

	/* Starting point */
	x = gsl_vector_alloc(np);				/* allocate space for a vector of length 3 */
	gsl_vector_set(x, 0, pattern->alpha);   /* set to starting guess, input (alpha,beta,gamma) */
	gsl_vector_set(x, 1, pattern->beta);
	gsl_vector_set(x, 2, pattern->gamma);

	/* Initialize the minimizer */
	minex_func.f = &my_f;					/* the function that returns the value */
	minex_func.n = np;						/* number of parameters */
	minex_func.params = (void *)pattern;	/* a generic pointer that points to something with info for my_f() */
	/* Note: my_func.df, and my_func.fdf are not needed since simplex method does not use derivatives */

	s = gsl_multimin_fminimizer_alloc(T, np);	/* allocate for a 2d simplex minimizer named 's' */

	/* Set starting point for the minimizer:
		for the minimizer 's',
		working on a function 'minex_func'
		starting from the initial guess x[2]
		using step sizes of ss[2]  */
	gsl_multimin_fminimizer_set(s, &minex_func, x, ss);

#if (DEBUG)		/******************************************************************************/
	fprintf(fout, "start minimization with Euler angles (%15.10f, %15.10f, %15.10f) (deg.)\n", \
		pattern->alpha*180/M_PI,pattern->beta*180/M_PI,pattern->gamma*180/M_PI);
#endif			/******************************************************************************/

	/* Iterate the minimizer.  This looop iterates the minimization method until it is good enough */
	iter = 0;
	do {
		iter++;
		status = gsl_multimin_fminimizer_iterate(s);	/* do one iteration of the minimizer */
		if (status) break;								/* status=0 is OK, check here for an error */

		size = gsl_multimin_fminimizer_size (s); 		/* get 'size' of minimizer */
		status = gsl_multimin_test_size (size, epsAbs);	/* test size against epsAbs */

#if (DEBUG>1)		/******************************************************************************/
		fprintf(fout,"%5ld  f(", iter);
		for (i=0; i<np; i++) {
			fprintf(fout,"%.8f", gsl_vector_get(s->x, i));
			if (i<np-1) fprintf(fout,", ");
		}
		fprintf(fout,") = %15.8g size = %.5g\n", s->fval, size);
#endif			/******************************************************************************/
	} while (status == GSL_CONTINUE && iter < (size_t)maxIter); /* stop after maxIter regardless of size */

#if (DEBUG)		/******************************************************************************/
	if (status == GSL_SUCCESS) fprintf(fout,"converged to minimum successfully\n");
	else fprintf(fout,"did not converge successfully, status = %d\n",status);
	fprintf(fout, "  optimization changed Euler angles by (%15.10f, %15.10f, %15.10f) (deg.)\n", \
		(gsl_vector_get(s->x,0)-pattern->alpha)*180/M_PI, (gsl_vector_get(s->x,1)-pattern->beta)*180/M_PI, \
		(gsl_vector_get(s->x,2)-pattern->gamma)*180/M_PI);
#endif			/******************************************************************************/

	pattern->alpha = gsl_vector_get(s->x,0);			/* set to final values, resultant (alpha,beta,gamma) returned */
	pattern->beta = gsl_vector_get(s->x,1);
	pattern->gamma = gsl_vector_get(s->x,2);

#if (DEBUG)		/******************************************************************************/
	fprintf(fout, "after minimization, Euler angles are (%15.10f, %15.10f, %15.10f) (deg.)\n", \
		pattern->alpha*180/M_PI,pattern->beta*180/M_PI,pattern->gamma*180/M_PI);
#endif			/******************************************************************************/

	/* recalculate the err[] list in pattern */
	MatrixCopy33(recip,pattern->xtal.recip);		/* copy unrotated reciprocal lattice into recip */
	EulerMatrix(pattern->alpha,pattern->beta,pattern->gamma,M_Euler);/* rotation matrix from Euler angles */
	MatrixMultiply33(M_Euler,recip,recip);			/* rotate recip by Euler angles */
	for (i=0;i<(size_t)(pattern->Ni);i++) {
		VECTOR_COPY3(hkl,pattern->hkls[i]);			/* convert integer hkl to float for following MatrixMultiply31 */
		MatrixMultiply31(recip,hkl,G);
		normalize3(G);								/* normalized gvec of hkl[i] rotated by Euler angles */
		dot = dot3(G,pattern->Ghat[i]);				/* note, Ghat here is the measured, not calculated G */
		(pattern->err)[i] = acos(MIN(dot,1.0));
	}

	gsl_vector_free(x);									/* free vector of x,y positions */
	gsl_vector_free(ss);								/* free vector of step sizes */
	gsl_multimin_fminimizer_free (s);					/* free the minimizer */
	return status;
}

max_stars_repo_name: carterbox/cold

max_issues_repo_name: carterbox/cold

max_forks_repo_name: carterbox/cold


text:
// MIT License
//
// Copyright (c) 2020 SunnyCase
//
// 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
#include "chinodef.h"
#include <atomic>
#include <chino/error.h>
#include <chino/result.h>
#include <gsl/gsl-lite.hpp>
#include <string_view>

namespace chino::io
{
enum class std_handles
{
    in,
    out,
    err
};

enum class create_disposition
{
    create_always = 2,
    create_new = 1,
    open_always = 4,
    open_existing = 3,
    truncate_existing = 5
};

result<void, error_code> alloc_console() noexcept;
handle_t get_std_handle(std_handles type) noexcept;

result<handle_t, error_code> open(const insert_lookup_object_options &options, create_disposition create_disp = create_disposition::open_existing) noexcept;
result<size_t, error_code> read(handle_t file, gsl::span<gsl::byte> buffer) noexcept;
result<void, error_code> write(handle_t file, gsl::span<const gsl::byte> buffer) noexcept;
result<void, error_code> close(handle_t file) noexcept;
}

ext: h

max_stars_repo_name: chino-os/chino-os

max_issues_repo_name: dotnetGame/chino-os

max_forks_repo_name: dotnetGame/chino-os


text:
#ifndef TEST_ALL
#define TEST_ALL
//#include <cblas.h>

#include "test_matrix_util.h"
#include "test_spn.h"

int test_all(int result);

#endif

max_stars_repo_name: ThayaFluss/cnl

max_issues_repo_name: ThayaFluss/cnl

max_forks_repo_name: ThayaFluss/cnl


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

void
cblas_ssyr2k (const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
              const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
              const float alpha, const float *A, const int lda,
              const float *B, const int ldb, const float beta, float *C,
              const int ldc)
{
#define BASE float
#include "source_syr2k_r.h"
#undef BASE
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/**
  * Copyright 2019 José Manuel Abuín Mosquera <josemanuel.abuin@usc.es>
  *
  * This file is part of Matrix Market Suite.
  *
  * Matrix Market Suite 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.
  *
  * Matrix Market Suite 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 Matrix Market Suite. If not, see <http://www.gnu.org/licenses/>.
  */

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

#include "ConjugateGradient.h"
#include "ConjugateGradientSolver.h"
#include "ConjugateGradientSolverBasic.h"

/* Include redundant execution header. */
#include "../../../include/ourRMTlib.h"

void usageConjugateGradient(){

	fprintf(stderr, "\n");
	fprintf(stderr, "Usage: MM-Suite ConjugateGradient [options] <input-matrix> <input-vector>\n");
	fprintf(stderr, "Algorithm options:\n\n");
	fprintf(stderr, "       -i INT        Iteration number. Default: number of matrix rows * 2\n");
	fprintf(stderr, "       -b            Uses basic operations instead of optimized BLAS libraries. Default: False\n");
	fprintf(stderr, "\nInput/output options:\n\n");
	fprintf(stderr, "       -o STR        Output file name. Default: stdout\n");
	fprintf(stderr, "       -r            Input format is row per line. Default: False\n");
	fprintf(stderr, "\nPerformance options:\n\n");
	fprintf(stderr, "       -t INT        Number of threads to use in OpenBLAS. Default: 1\n");
	fprintf(stderr, "\n");

}

//int ConjugateGradient(int argc, char *argv[]) {
int main(int argc, char *argv[]) {

	int 			ret_code;
	int 			option;

	unsigned long 		*II;
	unsigned long 		*J;
	double 			*A;

	unsigned long 		M;
	unsigned long 		N;
	unsigned long long 	nz;

	double 			*b;
	unsigned long 		M_Vector;
	unsigned long 		N_Vector;
	unsigned long long 	nz_vector;

	char			*outputFileName = NULL;
	int			iterationNumber = 0;

	char			*inputMatrixFile = NULL;
	char			*inputVectorFile = NULL;

	int			numThreads = 1;

	int			inputFormatRow = 0;
	int			basicOps = 0;
	int			F=0;

	while ((option = getopt(argc, argv,"bro:i:t:f:")) >= 0) {
		switch (option) {
			case 'o' :
				//free(outputFileName);

				outputFileName = (char *) malloc(sizeof(char)*strlen(optarg)+1);
				strcpy(outputFileName,optarg);

				break;

			case 'i' :
				iterationNumber = atoi(optarg);
				break;

			case 'r':
				inputFormatRow = 1;
				break;
			case 'b':
				basicOps = 1;
				break;

			case 't':
				numThreads = atoi(optarg);
				break;

			case 'f':
				F = atoi(optarg);
				break;

			default: break;
		}

	}

	if ((optind + 2 > argc) || (optind + 3 <= argc)) {
		usageConjugateGradient();
		return -1;
	}

	//openblas_set_num_threads(numThreads);

	if(outputFileName == NULL) {
		outputFileName = (char *) malloc(sizeof(char)*7);
		sprintf(outputFileName,"stdout");
	}

	start_timer();

	inputMatrixFile = (char *)malloc(sizeof(char)*strlen(argv[optind])+1);
	inputVectorFile = (char *)malloc(sizeof(char)*strlen(argv[optind+1])+1);

	strcpy(inputMatrixFile,argv[optind]);
	strcpy(inputVectorFile,argv[optind+1]);

	//Read matrix

	//Read matrix
	if(inputFormatRow){

		if(!readDenseCoordinateMatrixRowLine(inputMatrixFile,&II,&J,&A,&M,&N,&nz)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return -1;
		}

		//writeDenseVector("stdout",values,M,N,nz);

	}
	else {
		if(!readDenseCoordinateMatrix(inputMatrixFile,&II,&J,&A,&M,&N,&nz)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return -1;
		}
	}

	//Read vector
	if(!readDenseVector(inputVectorFile, &b,&M_Vector,&N_Vector,&nz_vector)){
		fprintf(stderr, "[%s] Can not read Vector\n",__func__);
		return -1;
	}

    //double *y=(double *) malloc(nz_vector * sizeof(double));
    fprintf(stderr,"[%s] Solving system using conjugate gradient method\n",__func__);
    double t_real = realtime();

	// Vector to store result
	double *x=(double *) calloc(nz_vector,sizeof(double));

	if(basicOps){
		ret_code = ConjugateGradientSolverBasic(II,J,A,M,N,nz,b,M_Vector,N_Vector,nz_vector, x, iterationNumber);
	}
	else{
		# ifdef ENABLE_RMT
			  activateRMT("f-L-lp-lp-2d-l-l-l-1dC-l-l-l-1dC-i-i-ipC", &ConjugateGradientSolver, 14, II, J, A, M, N, M, N, nz , b, N, M_Vector, N_Vector,
																	  nz_vector, x, nz_vector, iterationNumber, F, &ret_code);

		# else
			ret_code = ConjugateGradientSolver(II,J,A,M,N,nz,b,M_Vector,N_Vector,nz_vector, x, iterationNumber, F, &ret_code);
		# endif

	}

	fprintf(stderr, "\n[%s] Time spent in Conjugate Gradient: %.6f sec\n", __func__, realtime() - t_real);

	if(ret_code){

		writeDenseVector(outputFileName, x,M_Vector,N_Vector,nz_vector);

	}
	else{
		fprintf(stderr,"[%s] Error executing ConjugateGradientSolver\n",__func__);
		return -1;

	}

	end_timer();

	return 0;
}

max_stars_repo_name: sanemarslan/RMTLib

max_issues_repo_name: sanemarslan/RMTLib

max_forks_repo_name: sanemarslan/RMTLib


text:
/*
 * @Description  : Phase Class store the information of a phase in a temperature period
 */
#ifndef PHASES_H
#define PHASES_H
#include <string>
#include "VTypes.h"
#include <gsl/gsl_spline.h>

typedef std::vector<gsl_spline*> VSP;

class Phase
{
private:
    static int NPhases;
    int key; // Phase ID
    int DimX; // Dimension of X
    int DimData; // Dimension of T
    VVD X; // Store the Field information at every temperature
    VD T; // Temperature
    VVD dXdT; // The derivative of X respect to Temperature
    SI low_trans;
    SI high_trans;

    gsl_interp_accel *acc;
    VSP inters;

public:
    Phase();
    Phase(int _key, VVD _X, VD _T, VVD _dXdT);
    Phase(const Phase &Ph);
    Phase(const Phase &ph1, const Phase &ph2); //Merge two phases // ! When using this constructor, please make sure one of the two phases are redundant. We don't do such check
    ~Phase();

    void SetPhase(int _key, VVD _X, VD _T, VVD _dXdT);
    void SetUpInterpolation();
    void FreeInterpolation();
    VD valAt(double _T);
    int GetKey() const {return key;}
    double GetTmin() const {return T.front();}
    double GetTmax() const {return T.back();}
    VD GetXatTmin() const {VD Xmin; for(int i = 0; i < DimX; i++ ) {Xmin.push_back(X[i].front());} return Xmin;}
    VD GetXatTmax() const {VD Xmax; for(int i = 0; i < DimX; i++ ) {Xmax.push_back(X[i].back());} return Xmax;}
    void addLowTrans(int _key) {low_trans.insert(_key);}
    void eraseLowTrans(int _key) {low_trans.erase(_key);}
    void addHighTrans(int _key) {high_trans.insert(_key);}
    void eraseHighTrans(int _key) {high_trans.erase(_key);}
    void addLinkFrom(Phase *other_phase);
    std::string repr();
    std::string repr_full();

    VVD GetX() const {return X;}
    VVD GetdXdT() const {return dXdT;}
    VD GetT() const {return T;}
    SI GetLowTrans() const {return low_trans;}
    SI GetHighTrans() const {return high_trans;}
    static int GetNPhases() {return NPhases;}
};

#endif


text:
/**
 * \file include/shg/numalg.h
 * Numerical algorithms.
 */

#ifndef SHG_NUMALG_H
#define SHG_NUMALG_H

#include <functional>
#include <limits>
#include <stdexcept>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_multimin.h>
#include <shg/fcmp.h>
#include <shg/matrix.h>
#include <shg/vector.h>

namespace SHG {

/**
 * Degree of polynomial in one variable.
 *
 * \returns \f$\max \{i: a_i \neq 0, i = 0, \ldots, m - 1 \}\f$ if it
 * exists, otherwise -1.
 *
 * \throws std::invalid_argument if a.size() == 0 or a.size() > 32768.
 */
template <template <class> class Indexed_container, class T>
int degree_of_polynomial(const Indexed_container<T>& a, T eps);

/**
 * Calculates real roots of polynomial a[n] x^n + ... + a[2] x^2 +
 * a[1] x^1 + a[0]. Calculated roots are put in x[0] < x[1] < ... <
 * x[m - 1], where m = x.size().
 *
 * \throws std::invalid_argument if degree of polynomial is less than
 * 1
 * \throws std::runtime_error if the GSL routine fails
 */
void real_roots(const Vecdouble& a, Vecdouble& x);

template <template <class> class Indexed_container, class T>
int degree_of_polynomial(const Indexed_container<T>& a, T eps) {
     static_assert(std::numeric_limits<int>::max() >= 32767,
                   "According to the standard, INT_MAX must be at "
                   "least 32767.");
     auto k = a.size();
     if (k == 0 || k > 32768u)
          throw std::invalid_argument(
               "invalid polynomial in degree_of_polynomial");
     while (k-- > 0)
          if (fane(a[k], T(0), eps))
               return k;
     return -1;
}

/**
 * Solves the equation \f$a_0 + a_1x + \ldots + a_{n - 1}x^{n - 1} =
 * 0\f$, \f$n \geq 1\f$.
 */
void solve_polynomial(const Vecdouble& a, Veccomplex& x);

/**
 * Solves system of linear equations.
 */
void solve_linear(Matdouble& a, Vecdouble& x);

/** Light wrapper for GSL minimizer. */
class Minimizer_base {
public:
     Minimizer_base(const gsl_min_fminimizer_type* T);
     Minimizer_base(const Minimizer_base&) = delete;
     Minimizer_base& operator=(const Minimizer_base&) = delete;
     virtual ~Minimizer_base();
     gsl_min_fminimizer* get() const;

private:
     gsl_min_fminimizer* s_;
};

/** Wrapper for GSL minimizer. */
class Minimizer : Minimizer_base {
public:
     Minimizer(const gsl_min_fminimizer_type* T);
     Minimizer(const Minimizer&) = delete;
     Minimizer& operator=(const Minimizer&) = delete;

     bool is_set() const { return is_set_; }

     void set(gsl_function* f, double x_minimum, double x_lower,
              double x_upper);

     int iterate(int max_iter, double epsabs, double epsrel);

     int iter() const { return iter_; }

     double x_minimum() const { return x_minimum_; }
     double x_lower() const { return x_lower_; }
     double x_upper() const { return x_upper_; }
     double f_minimum() const { return f_minimum_; }
     double f_lower() const { return f_lower_; }
     double f_upper() const { return f_upper_; }

private:
     gsl_min_fminimizer* s_{get()};
     int iter_{0};
     bool is_set_{false};
     double x_minimum_{};
     double x_lower_{};
     double x_upper_{};
     double f_minimum_{};
     double f_lower_{};
     double f_upper_{};
};

/**
 * Uniform search for minimum of the function of one variable.
 *
 * The method search() searches uniformly the interval \f$[a, b]\f$
 * for minima of the function \f$f\f$. The method is inspired by
 * \cite bronsztejn-siemiendiajew-musiol-muhlig-2004, section
 * 18.2.4.1, point 1, page 946.
 *
 * First, \f$n\f$ is calculated as \f[ n = \min \{ 2, \lceil
 * \frac{2}{\epsilon} (b - a) \rceil \}, \f] and then the interval
 * \f$[a, b]\f$ is divided into \f$n\f$ intervals \f$[x_k, x_{k +
 * 1}]\f$, \f$k = 0, \ldots n - 1\f$, where \f$x_k = a + kh\f$, \f$h =
 * (b - a) / n\f$. The values of \f$f\f$ are calculated in each
 * \f$x_k\f$ and all values \f[ (x_k, x_{k + 1}, x_{k + 2}, f(x_k),
 * f(x_{k + 1}), f(x_{k + 2})) \f] such that \f[ f(x_k) > f(x_{k + 1})
 * < f(x_{k + 2}) \f] are reported in \f$\var{result}\f$. The length
 * of the interval \f$[x_k, x_{k + 2}]\f$ is not greater than
 * \f$\epsilon\f$. The function is evaluated \f$n + 1\f$times.
 *
 * \throws std::invalid_argument if on entry \f$a \geq b\f$ or
 * \f$\epsilon \leq 0\f$.
 */
struct Uniform_search_for_minimum {
     struct Result {
          double x_lower;
          double x_minimum;
          double x_upper;
          double f_lower;
          double f_minimum;
          double f_upper;
     };
     std::vector<Result> result{};
     void search(std::function<double(double)> f, double a, double b,
                 double eps);
};

/** Light wrapper for GSL minimizer for functions of several
    variables. */
class Multimin_fminimizer_base {
public:
     Multimin_fminimizer_base(const gsl_multimin_fminimizer_type* T,
                              size_t n);
     Multimin_fminimizer_base(const Multimin_fminimizer_base&) =
          delete;
     Multimin_fminimizer_base& operator=(
          const Multimin_fminimizer_base&) = delete;
     virtual ~Multimin_fminimizer_base();
     gsl_multimin_fminimizer* get() const;

private:
     gsl_multimin_fminimizer* s_;
};

/** Wrapper for GSL minimizer for functions of several variables. */
class Multimin_fminimizer : public Multimin_fminimizer_base {
public:
     Multimin_fminimizer(const gsl_multimin_fminimizer_type* T,
                         size_t n);
     Multimin_fminimizer(const Multimin_fminimizer&) = delete;
     Multimin_fminimizer& operator=(const Multimin_fminimizer&) =
          delete;
     virtual ~Multimin_fminimizer();

     bool is_set() const { return is_set_; }

     int iter() const { return iter_; }

     void set(gsl_multimin_function* f, const std::vector<double>& x,
              const std::vector<double>& step);
     int iterate(int max_iter, double eps);

     const std::vector<double>& x_minimum() const {
          return x_minimum_;
     }
     double f_minimum() const { return f_minimum_; }

private:
     gsl_multimin_fminimizer* s_{get()};
     int iter_{0};
     bool is_set_{false};
     std::size_t n_;
     gsl_vector* x_;
     gsl_vector* ss_;
     std::vector<double> x_minimum_;
     double f_minimum_;
};

}  // namespace SHG

#endif

max_stars_repo_name: shgalus/shg

max_issues_repo_name: shgalus/shg

max_forks_repo_name: shgalus/shg


text:
/* randist/shuffle.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Inline swap and copy functions for moving objects around */

static inline
void swap (void * base, size_t size, size_t i, size_t j)
{
  register char * a = size * i + (char *) base ;
  register char * b = size * j + (char *) base ;
  register size_t s = size ;

  if (i == j)
    return ;

  do
    {
      char tmp = *a;
      *a++ = *b;
      *b++ = tmp;
    }
  while (--s > 0);
}

static inline void
copy (void * dest, size_t i, void * src, size_t j, size_t size)
{
  register char * a = size * i + (char *) dest ;
  register char * b = size * j + (char *) src ;
  register size_t s = size ;

  do
    {
      *a++ = *b++;
    }
  while (--s > 0);
}

/* Randomly permute (shuffle) N indices

   Supply an array x[N] with nmemb members, each of size size and on
   return it will be shuffled into a random order.  The algorithm is
   from Knuth, SemiNumerical Algorithms, v2, p139, who cites Moses and
   Oakford, and Durstenfeld */

void
gsl_ran_shuffle (const gsl_rng * r, void * base, size_t n, size_t size)
{
  size_t i ;

  for (i = n - 1; i > 0; i--)
    {
      size_t j = gsl_rng_uniform_int(r, i+1); /* originally (i + 1) * gsl_rng_uniform (r) */

      swap (base, size, i, j) ;
    }
}

int
gsl_ran_choose (const gsl_rng * r, void * dest, size_t k, void * src,
                 size_t n, size_t size)
{
  size_t i, j = 0;

  /* Choose k out of n items, return an array x[] of the k items.
     These items will prevserve the relative order of the original
     input -- you can use shuffle() to randomize the output if you
     wish */

  if (k > n)
    {
      GSL_ERROR ("k is greater than n, cannot sample more than n items",
                 GSL_EINVAL) ;
    }

  for (i = 0; i < n && j < k; i++)
    {
      if ((n - i) * gsl_rng_uniform (r) < k - j)
        {
          copy (dest, j, src, i, size) ;
          j++ ;
        }
    }

  return GSL_SUCCESS;
}

void
gsl_ran_sample (const gsl_rng * r, void * dest, size_t k, void * src,
                size_t n, size_t size)
{
  size_t i, j = 0;

  /* Choose k out of n items, with replacement */

  for (i = 0; i < k; i++)
    {
      j = gsl_rng_uniform_int (r, n);  /* originally n * gsl_rng_uniform (r) */

      copy (dest, i, src, j, size) ;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/** @file */
#ifndef __CCL_CORE_H_INCLUDED__
#define __CCL_CORE_H_INCLUDED__

#include <stdbool.h>
#include <stdio.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

CCL_BEGIN_DECLS

/**
 * Struct containing the parameters defining a cosmology
 */
typedef struct ccl_parameters {

  // Densities: CDM, baryons, total matter, neutrinos, curvature
  double Omega_c; /**< Density of CDM relative to the critical density*/
  double Omega_b; /**< Density of baryons relative to the critical density*/
  double Omega_m; /**< Density of all matter relative to the critical density*/
  double Omega_k; /**< Density of curvature relative to the critical density*/
  double sqrtk; /**< Square root of the magnitude of curvature, k */ //TODO check
  int k_sign; /**<Sign of the curvature k */

  // Dark Energy
  double w0;
  double wa;

  // Hubble parameters
  double H0;
  double h;

  // Neutrino properties

  double Neff; // Effective number of relativistic neutrino species in the early universe.
  int N_nu_mass; // Number of species of neutrinos which are nonrelativistic today
  double N_nu_rel;  // Number of species of neutrinos which are relativistic  today
  double *mnu;  // total mass of massive neutrinos (This is a pointer so that it can hold multiple masses.)
  double sum_nu_masses; // sum of the neutrino masses.
  double Omega_n_mass; // Omega_nu for MASSIVE neutrinos
  double Omega_n_rel; // Omega_nu for MASSLESS neutrinos

  //double Neff_partial[CCL_MAX_NU_SPECIES];
  //double mnu[CCL_MAX_NU_SPECIES];

  // Primordial power spectra
  double A_s;
  double n_s;

  // Radiation parameters
  double Omega_g;
  double T_CMB;

  // BCM baryonic model parameters
  double bcm_log10Mc;
  double bcm_etab;
  double bcm_ks;

  // Derived parameters
  double sigma8;
  double Omega_l;
  double z_star;

  //Modified growth rate
  bool has_mgrowth;
  int nz_mgrowth;
  double *z_mgrowth;
  double *df_mgrowth;
} ccl_parameters;

/**
 * Struct containing references to gsl splines for distance and acceleration calculations
 */
typedef struct ccl_data{
  // These are all functions of the scale factor a.

  // Distances are defined in Mpc
  double growth0;
  gsl_spline * chi;
  gsl_spline * growth;
  gsl_spline * fgrowth;
  gsl_spline * E;
  gsl_spline * achi;

  // All these splines use the same accelerator so that
  // if one calls them successively with the same a value
  // they will be much faster.
  gsl_interp_accel *accelerator;
  gsl_interp_accel *accelerator_achi;
  gsl_interp_accel *accelerator_m;
  gsl_interp_accel *accelerator_d;
  //TODO: it seems like we're not really using this accelerator, and we should
  gsl_interp_accel *accelerator_k;

  // Function of Halo mass M

  gsl_spline * logsigma;
  gsl_spline * dlnsigma_dlogm;

  // splines for halo mass function
  gsl_spline * alphahmf;
  gsl_spline * betahmf;
  gsl_spline * gammahmf;
  gsl_spline * phihmf;
  gsl_spline * etahmf;

  // These are all functions of the wavenumber k and the scale factor a.
  gsl_spline2d * p_lin;
  gsl_spline2d * p_nl;
  double k_min_lin; //k_min  [1/Mpc] <- minimum wavenumber that the power spectrum has been computed to
  double k_min_nl;
  double k_max_lin;
  double k_max_nl;
} ccl_data;

/**
 * Sturct containing references to instances of the above structs, and boolean flags of precomputed values.
 */
typedef struct ccl_cosmology
{
  ccl_parameters    params;
  ccl_configuration config;
  ccl_data          data;

  bool computed_distances;
  bool computed_growth;
  bool computed_power;
  bool computed_sigma;
  bool computed_hmfparams;

  int status;
  //this is optional - less tedious than tracking all numerical values for status in error handler function
  char status_message[500];

  // other flags?
} ccl_cosmology;

// Label for whether you are passing a pointer to a sum of neutrino masses or a pointer to a list of 3 masses.
typedef enum ccl_mnu_convention {
  ccl_mnu_list=0,   // you pass a list of three neutrino masses
  ccl_mnu_sum =1,  // sum, defaults to splitting with normal hierarchy
  ccl_mnu_sum_inverted = 2, //sum, split with inverted hierarchy
  ccl_mnu_sum_equal = 3, //sum, split into equal masses
  // More options could be added here
} ccl_mnu_convention;

// Initialization and life cycle of objects
void ccl_cosmology_read_config(void);
ccl_cosmology * ccl_cosmology_create(ccl_parameters params, ccl_configuration config);

/* Internal function to set the status message safely. */
void ccl_cosmology_set_status_message(ccl_cosmology * cosmo, const char * status_message, ...);

// User-facing creation routines
/**
 * Create a cosmology
 * @param Omega_c Omega_c
 * @param Omega_b Omega_b
 * @param Omega_k Omega_k
 * @param Neff Number of relativistic neutrino species in the early universe
 * @param mnu neutrino mass, either sum or list of length 3
 * @param mnu_type determines neutrino mass convention (ccl_mnu_list, ccl_mnu_sum, ccl_mnu_sum_inverted, ccl_mnu_sum_equal)
 * @param w0 Dark energy EoS parameter
 * @param wa Dark energy EoS parameter
 * @param h Hubble constant in units of 100 km/s/Mpc
 * @param norm_pk the normalization of the power spectrum, either A_s or sigma8
 * @param n_s the power-law index of the power spectrum
 * @param bcm_log10Mc log10 cluster mass, one of the parameters of the BCM model
 * @param bcm_etab ejection radius parameter, one of the parameters of the BCM model
 * @param bcm_ks wavenumber for the stellar profile, one of the parameters of the BCM model
 * @param nz_mgrowth the number of redshifts where the modified growth is provided
 * @param zarr_mgrowth the array of redshifts where the modified growth is provided
 * @param dfarr_mgrowth the modified growth function vector provided
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * For specific cases see documentation for ccl_error.c
 * @return void
 */
ccl_parameters ccl_parameters_create(double Omega_c, double Omega_b, double Omega_k,
                     double Neff, double* mnu, ccl_mnu_convention mnu_type,
                     double w0, double wa, double h, double norm_pk,
                     double n_s, double bcm_log10Mc, double bcm_etab, double bcm_ks,
                     int nz_mgrowth,double *zarr_mgrowth,
                     double *dfarr_mgrowth, int *status);

/* ------- ROUTINE: ccl_parameters_create_flat_lcdm --------
INPUT: some cosmological parameters needed to create a flat LCDM model
TASK: call ccl_parameters_create to produce an LCDM model
*/
ccl_parameters ccl_parameters_create_flat_lcdm(double Omega_c, double Omega_b, double h,
					       double norm_pk, double n_s, int *status);

/**
 * Free a parameters struct
 * @param params ccl_parameters struct
 * @return void
 */
void ccl_parameters_free(ccl_parameters * params);

/**
 * Write a cosmology parameters object to a file in yaml format, .
 * @param params Cosmological parameters
 * @param filename Name of file to create and write
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * @return void
 */
void ccl_parameters_write_yaml(ccl_parameters * params, const char * filename, int * status);

/**
 * Read a cosmology parameters object from a file in yaml format, .
 * @param filename Name of existing file to read from
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 * @return cosmo Cosmological parameters
 */
ccl_parameters ccl_parameters_read_yaml(const char * filename, int *status);

/**
 * Free a cosmology struct
 * @param cosmo Cosmological parameters
 * @return void
 */
void ccl_cosmology_free(ccl_cosmology * cosmo);

CCL_END_DECLS

#endif

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
/*                                                                                */
/* rt - retention time (predictor) v4.0.20160301                                                      */
/*                                                                                */
/* (c) Magnus Palmblad, Division of Ion Physics, Uppsala University, 2002-                            */
/*                                                                                */
/*                                                                                */
/* Version 4.0 of rt now also accepts pepXML input instead of just tab-delimited files, recognizes    */
/* the input format automatically and can generate output (retention coefficients) either as a tab-   */
/* delimited text file or an XML structure. Known bugs: This version does not extract modifications   */
/* from the pepXML files and only works with the pepXML files generated from PeptideProphet.          */
/* The program does not perform detailed checks on the input, i.e. peptide lengths and syntax errors. */
/*                                                                                */
/* Usage: rt -i <training set in tab delimited or pepXML format> [-f <tab>|<XML> output format)]>     */
/*                                                                                */
/* compile with gcc -o rt rt3.c -lgsl -lgslcblas -lpepXML -LpepXMLLib                                 */
/*                                                                                */
/*                                                                                */
/* The mzIdentML support is extremely "quick-and-dirty" but work with tested idconvert-ed pepXML      */
/*                                                                                */

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include "pepXMLLib/pepXMLReader.h"

#define AMINO_ACIDS "ARNDCEQGHILKMFPSTWYVmsty"
#define MAX_PEPTIDE_LENGTH 100

int main(int argc, char *argv[])
{
  char line[250];
  char infile[85], outfile[85]; /* for reading input and writing output */
  char slask[85];

  typedef struct {
    char sequence[MAX_PEPTIDE_LENGTH];
    int composition[24];
    float retention;
  } training_set_type;

  training_set_type *training_set;

  char output[256];
  char output_format[256];
  char input_format[20];

  char *p, q[100], **PEP;

  pmsms_pipeline_analysis pepXML_file;
  pdelegate_list dlgl;
  pdelegate_type dlg;

  FILE *inp, *outp;

  int a, b, i, j, k, l, m, n, seqlen, n_training,length;
  double xi, yi, ei, chisq;
  gsl_matrix *X, *cov;
  gsl_vector *y, *w, *c;
  gsl_multifit_linear_workspace *work;

  /* parsing command line parameters */

  if ((argc<2)||(argc>5))
    {
      printf("usage: rt -i <training set in tab delimited or PeptideProphet pepXML format> [-f <tab>|<XML> output format)] (type rt --help for more information\n");
      return -1;
    }

  if( (strcmp(argv[1],"--help")==0) || (strcmp(argv[1],"-help")==0) || (strcmp(argv[1],"-h")==0) ) /* want help? */
    {
      printf("rt - (c) Magnus Palmblad 2002-2012\n\nusage: rt -i <training set in tab-delimited or pepXML format (after PeptideProphet)> [-f <tab>|<XML> output format)]\nrt trains a linear model for predicting retention time, or in fact any additive property as function of amino acid composition. The input should either be a tab-delimited text file with peptide sequences in the first column and a numerical value (retention time) in the second column, or a pepXML file generated by PeptideProphet in the Trans-Proteomic Pipeline. For more information, see http://www.ms-utils.org/rt or e-mail magnus.palmblad@gmail.com\n");
      return 0;
    }

  for(i=1;i<argc;i++) {
    if( (argv[i][0]=='-') && (argv[i][1]=='i') ) strcpy(infile,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( (argv[i][0]=='-') && (argv[i][1]=='f') ) strcpy(output_format,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
  }

  if( (toupper(output_format[0])!='T') && (toupper(output_format[0])!='X') && (toupper(output_format[0])!='C') ) {
    printf("usage: rt -i <training set in tab delimited or pepXML format> [-f <tab>|<XML> output format)]\n");
    return -1;
  }

  strcpy(input_format,"tab"); /* default to tab-delimited input */

  /* reading training set file */

  //printf("checking input format and scanning training set in file %s...", infile); fflush(stdout);
  n=0;
  if ((inp = fopen(infile, "r"))==NULL) {printf("error opening training set file\n");return -1;}
  fgets(line, 250, inp);
  if(line[0]=='<') {strcpy(input_format,"pepXML"); n++;}
  fgets(line, 250, inp);
  while(line[i]) {line[i]=toupper(line[i]); i++;}
  if(strstr(line,"MZIDENTML")!=NULL) {strcpy(input_format,"mzIdentML"); n++;}

  while (fgets(line, 250, inp) != NULL) n++;

  PEP=(char**)malloc(sizeof(char*)*n);
  for(i=0;i<10000;i++) PEP[i]=(char*)malloc(sizeof(char)*200);

  if(strcmp(input_format,"tab")==0)
    {
      n_training=n;

      // printf("done\nallocating memory for %i peptides in training set...", n_training); fflush(stdout);
      training_set=(training_set_type*)malloc(sizeof(training_set_type)*n_training);

      // printf("done\nreading training set in tab delimited format (peptide tab retention time) from file %s...", infile); fflush(stdout);
      if ((inp = fopen(infile, "r"))==NULL) {printf("error opening training set file\n");return -1;}
      n=0;

      while (fgets(line, 250, inp) != NULL)
	{
	  sscanf(line,"%s\t%f",training_set[n].sequence,&(training_set[n].retention));
	  n++;
	}
      n_training=n;
      fclose(inp);
      // printf("done\n");
    }

  if(strcmp(input_format,"pepXML")==0)
    {
      // printf("done\nscanning training set in pepXML format from file %s...", infile); fflush(stdout);
      n=0;
      pepXML_file = read_pepxml_file(infile, 0, 0, NULL);
      for (i=0; i<pepXML_file->run_summary_count; i++) {
	for (j=0; j<pepXML_file->run_summary_array[i].spectrum_query_count; j++) n++;
      }
      n_training=n;

      // printf("done\nallocating memory for %i peptides in training set...", n_training); fflush(stdout);
      training_set=(training_set_type*)malloc(sizeof(training_set_type)*n_training);

      n=0;
      // printf("reading training set in pepXML format from file %s...", infile); fflush(stdout);
      for (i=0; i<pepXML_file->run_summary_count; i++) {
	for (j=0; j<pepXML_file->run_summary_array[i].spectrum_query_count; j++) {
	  strcpy(training_set[n].sequence, pepXML_file->run_summary_array[i].spectrum_query_array[j].search_result_array[0].search_hit_array[0].peptide);
	  training_set[n].retention = pepXML_file->run_summary_array[i].spectrum_query_array[j].retention_time_sec;
	  n++;
	}
      }
    }

  if(strcmp(input_format,"mzIdentML")==0)
    {
      //printf("done\nscanning training set in mzIdentML format from file %s...", infile); fflush(stdout);
      n=0;
      if ((inp = fopen(infile, "r"))==NULL) {printf("error opening training set file\n");return -1;}
      while (fgets(line, 250, inp) != NULL)
	{
	  if(strstr(line,"PeptideEvidence")!=NULL) break;
	  if(strstr(line,"<Peptide id")!=NULL)
	    {
	      //PEP[a]=(char*)malloc(sizeof(char)*strlen(q));
	      p=strtok(line,"\""); p=strtok('\0',"\"");
	      a=atoi(p+4);
	      //printf("\n%i = ",a); fflush(stdout);
	      p=strtok('\0',"\"");
	      if(strstr(p,"<PeptideSequence>")!=NULL)
		{
		  p=strchr(p+10,'>')+1;
		  i=strspn(p,AMINO_ACIDS);
		  strncpy(q,p,i); q[i]='\0';
		  // printf("PEP_%i=%s\n",a,q); fflush(stdout);
		  strcpy(PEP[a],q);
		}
	      n++;
	    }
	}

      n_training=n;

      // printf("done\nallocating memory for %i peptides in training set...", n_training); fflush(stdout);
      training_set=(training_set_type*)malloc(sizeof(training_set_type)*n_training*10);

      n=0;
      while (fgets(line, 250, inp) != NULL)
	{
	  if(strstr(line,"peptide_ref")!=NULL)
	    {
	      p=strstr(line,"peptide_ref");
	      a=atoi(strtok(p+17,"\""));
	    }
	  if(strstr(line,"name=\"scan start time\" value")!=NULL)
	    {
	      p=strstr(line,"name=\"scan start time\" value");
	      i=strcspn(p+31,"\"");
	      strncpy(q,p+30,i); q[i]='\0';
	      // if(strlen(PEP[a])>6) {printf("%i, %s, %s\n",a,PEP[a],q); fflush(stdout);}
	      if(strlen(PEP[a])>6) {strcpy(training_set[n].sequence,PEP[a]); training_set[n].retention=atof(q); n++;}
	    }
	}
    }

  n_training=n;
  //return 0;

  /* allocating memory for training */

  X=gsl_matrix_alloc(n_training,25); /* 20 normal amino acids, pS, pT, pY, oxidized M and V0 */
  y=gsl_vector_alloc(n_training);
  w=gsl_vector_alloc(n_training);
  c=gsl_vector_alloc(25);
  cov=gsl_matrix_alloc(25,25);

  for(i=0;i<25;i++) gsl_vector_set(c,i,1.0);

  /* calculate training set X matrix (amino acid composition for each peptide in training set) */

  for(n=0;n<n_training;n++)
    {
      gsl_vector_set(y,n,training_set[n].retention);
      gsl_vector_set(w,n,0.36);

      for(i=0;i<24;i++) training_set[n].composition[i]=0;
      {
	for(i=0;i<strlen(training_set[n].sequence);i++)
	  {
	    a=24-strlen(strchr(AMINO_ACIDS,training_set[n].sequence[i]));
	    training_set[n].composition[a]++;
	  }
      }

      for(i=0;i<24;i++) gsl_matrix_set(X,n,i,training_set[n].composition[i]);
      gsl_matrix_set(X,n,24,1); /* to get V0 - the "offset" between runs */
    }

  work=gsl_multifit_linear_alloc(n_training,25);
  gsl_multifit_wlinear(X,w,y,c,cov,&chisq,work);
  gsl_multifit_linear_free(work);

  if(toupper(output_format[0])=='T') /* output as tab-delimited text */
   {
     // printf("\nbest fit to experimental retention data:\n");
     for(i=0;i<24;i++) printf("%c\t%f\n",AMINO_ACIDS[i],gsl_vector_get(c,(i)));
     printf("O\t%f\n\n",gsl_vector_get(c,(24)));

     return 0;
   }

  if(toupper(output_format[0])=='C') /* output CSV */
   {
     // printf("\nbest fit to experimental retention data:\n");
     for(i=0;i<24;i++) printf("%c,%f\n",AMINO_ACIDS[i],gsl_vector_get(c,(i)));
     printf("O,%f\n\n",gsl_vector_get(c,(24)));

     return 0;
   }

  if(toupper(output_format[0])=='X') /* output XML */
   {
     // printf("\nretenion model in XML structure:\n\n"); /* comment out in deployed versions */

     printf("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
     printf("<rt_model training_file=\"%s\" rt_version=\"4.0\" input_format=\"%s\">\n",infile,input_format);
     for(i=0;i<19;i++) printf(" <retention_coefficient amino_acid=\"%c\" value=\"%f\"/>\n",AMINO_ACIDS[i],gsl_vector_get(c,(i)));
     printf(" <retention_coefficient amino_acid=\"m\" modification=\"oxidation\" value=\"%f\"/>\n",gsl_vector_get(c,(19)));
     for(i=21;i<24;i++) printf(" <retention_coefficient amino_acid=\"%c\" modification=\"phosphorylation\" value=\"%f\"/>\n",AMINO_ACIDS[i],gsl_vector_get(c,(i)));
     printf(" <retention_coefficient amino_acid=\"O\" modification=\"constant offset\" value=\"%f\"/>\n",gsl_vector_get(c,(24)));
     printf("</rt_model>\n");

     return 0;
   }
}

max_stars_repo_name: veitveit/Automatic-Workflow-Composition

max_issues_repo_name: veitveit/Automatic-Workflow-Composition

max_forks_repo_name: veitveit/Automatic-Workflow-Composition


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

#include <gslib/std.h>

__gslib_begin__

enum json_tag
{
    jst_array,
    jst_table,
    jst_pair,
    jst_value,
};

class json_node_value;

class __gs_novtable json_node abstract
{
public:
    virtual ~json_node() {}
    virtual json_tag get_tag() const = 0;
    virtual const string& get_name() const = 0;
    virtual json_node* duplicate() const = 0;   /* delete */
};

struct json_node_hash
#if defined(_MSC_VER) && (_MSC_VER < 1914)
    : public std::unary_function<json_node*, size_t>
#endif
{
public:
    size_t operator()(const json_node* p) const
    {
        assert(p);
        return string_hash(p->get_name());
    }
};

struct json_node_equalto
#if defined(_MSC_VER) && (_MSC_VER < 1914)
    : public std::binary_function<json_node*, json_node*, bool>
#endif
{
public:
    bool operator()(const json_node* p1, const json_node* p2) const
    {
        assert(p1 && p2);
        return p1->get_name() == p2->get_name();
    }
};

typedef unordered_set<json_node*, json_node_hash, json_node_equalto> json_node_map;
typedef vector<json_node*> json_node_list;

class json_node_array:
    public json_node
{
public:
    virtual ~json_node_array();
    virtual json_tag get_tag() const override { return jst_array; }
    virtual const string& get_name() const override { return _name; }
    virtual json_node* duplicate() const override;

protected:
    string              _name;
    json_node_list      _array;

public:
    void set_name(const gchar* str) { _name.assign(str); }
    void set_name(const gchar* str, int len) { _name.assign(str, len); }
    void set_name(const string& str) { _name = str; }
    bool is_empty() const { return _array.empty(); }
    int get_childs() const { return (int)_array.size(); }
    json_node* at(int i) const { return _array.at(i); }
    int parse(const gchar* src, int len);
    json_node_list& get_container() { return _array; }
};

class json_node_table:
    public json_node
{
public:
    virtual ~json_node_table();
    virtual json_tag get_tag() const override { return jst_table; }
    virtual const string& get_name() const override { return _name; }
    virtual json_node* duplicate() const override;

protected:
    string              _name;
    json_node_map       _table;

public:
    void set_name(const gchar* str) { _name.assign(str); }
    void set_name(const gchar* str, int len) { _name.assign(str, len); }
    void set_name(const string& str) { _name = str; }
    bool is_empty() const { return _table.empty(); }
    int get_childs() const { return (int)_table.size(); }
    json_node* find(const string& name) const;
    int parse(const gchar* src, int len);
    json_node_map& get_container() { return _table; }
};

class json_node_value:
    public json_node
{
public:
    virtual json_tag get_tag() const override { return jst_value; }
    virtual const string& get_name() const override { return _strval; }
    virtual json_node* duplicate() const override;

protected:
    string              _strval;

public:
    void set_value_string(const gchar* str) { _strval.assign(str); }
    void set_value_string(const gchar* str, int len) { _strval.assign(str, len); }
    void set_value_string(const string& str) { _strval = str; }
    int parse(const gchar* src, int len);
};

class json_node_pair:
    public json_node
{
public:
    virtual json_tag get_tag() const override { return jst_pair; }
    virtual const string& get_name() const override { return _name; }
    virtual json_node* duplicate() const override;

protected:
    string              _name;
    json_node_value     _value;

public:
    void set_name(const gchar* str) { _name.assign(str); }
    void set_name(const gchar* str, int len) { _name.assign(str, len); }
    void set_name(const string& str) { _name = str; }
    void set_value(const gchar* str) { _value.set_value_string(str); }
    void set_value(const gchar* str, int len) { _value.set_value_string(str, len); }
    void set_value(const string& str) { _value.set_value_string(str); }
    const json_node_value& get_value() const { return _value; }
    const string& get_value_string() const { return _value.get_name(); }
    int parse(const gchar* src, int len);
};

class json_key:
    public json_node_pair
{
public:
    json_key(const string& name) { _name = name; }
    json_key(const gchar* str) { _name.assign(str); }
    json_key(const gchar* str, int len) { _name.assign(str, len); }
};

class json_parser
{
public:
    json_parser() { _root = nullptr; }
    ~json_parser() { destroy(); }
    bool parse(const gchar* src, int len);
    bool parse(const gchar* filename);
    void destroy();
    json_node* get_root() const { return _root; }

protected:
    json_node*          _root;
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#ifndef MCMC_AWAKE
#define MCMC_AWAKE

#include <stdarg.h>
#include <gbpRNG.h>
#include <gsl/gsl_linalg.h>

#define MCMC_COVERAGE_LINEAR 0
#define MCMC_COVERAGE_LOG 1
#define MCMC_NAME_SIZE 256

#define MCMC_MODE_SERIAL 1
#define MCMC_MODE_PARALLEL 2
#define MCMC_MODE_AUTOTUNE 8
#define MCMC_MODE_NO_MAP_WRITE 16
#define MCMC_MODE_REPORT_PROPS 32
#define MCMC_MODE_REFLECT_PRIORS 64
#define MCMC_MODE_MINIMIZE_IO 128
#define MCMC_MODE_AUTOTUNE_REINIT_TEMP 256
#define MCMC_MODE_ANALYZE_ALL_RUNS 512
#define MCMC_MODE_DEFAULT MCMC_MODE_AUTOTUNE | MCMC_MODE_SERIAL

#define MCMC_MAP_RETURN_GOOD 0
#define MCMC_MAP_RETURN_BAD 1

#define MCMC_DEFAULT_SUCCESS_TARGET 35.
#define MCMC_DEFAULT_SUCCESS_THRESH 5.
#define MCMC_DEFAULT_COVARIANCE_THRESH 10.
#define MCMC_DEFAULT_N_AUTOTUNE 1
#define MCMC_DEFAULT_N_AUTOTUNE_RANDOMIZE 0
#define MCMC_DEFAULT_N_AUTOTUNE_TEMPERATURE 100
#define MCMC_DEFAULT_N_AUTOTUNE_COVMTX (10 * MCMC->n_P * MCMC->n_P * (1 + (int)(100. / MCMC->success_target)))
#define MCMC_AUTOTUNE_CONVERGENCE_THRESH 1e-6

typedef struct MCMC_DS_info MCMC_DS_info;
struct MCMC_DS_info {
    char          name[MCMC_NAME_SIZE];
    int           n_D;
    int *         D;
    int           n_M;
    double *      M_target;
    double *      dM_target;
    double **     array;
    char **       array_name;
    void *        params;
    int           n_arrays;
    double *      M_best;
    double *      M_best_parameters;
    double *      M_peak_parameters;
    double *      M_min;
    double *      M_max;
    double *      M_lo_68;
    double *      M_hi_68;
    double *      M_lo_95;
    double *      M_hi_95;
    double *      M_avg;
    double *      dM_avg;
    MCMC_DS_info *next;
};

typedef struct MCMC_info MCMC_info;
struct MCMC_info {
    SID_Comm *comm;
    char      filename_output_root[SID_MAX_FILENAME_LENGTH];
    char      filename_output_dir[SID_MAX_FILENAME_LENGTH];
    char      problem_name[MCMC_NAME_SIZE];
    int (*map_P_to_M)(double *, MCMC_info *, double **);
    void (*compute_MCMC_ln_likelihood)(MCMC_info *MCMC,
                                       double **  M,
                                       double *   P,
                                       double *   ln_likeliood_DS,
                                       int *      n_DoF_DS,
                                       double *   ln_likeliood_all,
                                       int *      n_DoF_all);
    int           my_chain;
    int           n_chains;
    void *        params;
    char **       P_names;
    double *      P_init;
    double *      P_new;
    double *      P_last;
    double *      P_chain;
    double *      P_limit_min;
    double *      P_limit_max;
    double *      P_min;
    double *      P_max;
    double *      P_avg;
    double *      dP_avg;
    double *      P_best;
    double *      P_peak;
    double *      P_lo_68;
    double *      P_hi_68;
    double *      P_lo_95;
    double *      P_hi_95;
    double        ln_Pr_new;
    double        ln_Pr_last;
    double        ln_Pr_chain;
    int *         n_M;
    double **     M_new;
    double **     M_last;
    int           n_arrays;
    char **       array_name;
    double **     array;
    double *      V;
    gsl_matrix *  m;
    gsl_vector *  b;
    int           n_P;
    int           n_DS;
    int           n_M_total;
    double        temperature;
    int           n_avg;
    int           n_avg_covariance;
    int           n_iterations;
    int           n_iterations_burn;
    int           n_thin;
    int           coverage_size;
    int           flag_autocor_on;
    int           flag_no_map_write;
    int           flag_integrate_on;
    int           flag_analysis_on;
    int           flag_init_chain;
    int           first_map_call;
    int           first_link_call;
    size_t        n_map_calls;
    int           mode;
    int           n_success;
    int           n_fail;
    int           n_propositions;
    double        success_target;
    double        success_threshold;
    double        covariance_threshold;
    int           n_autotune;
    int           n_autotune_randomize;
    int           n_autotune_temperature;
    int           n_autotune_covariance;
    int           first_parameter_call;
    int           first_chain_call;
    int           first_likelihood_call;
    double        ln_likelihood_last;
    double        ln_likelihood_new;
    double        ln_likelihood_best;
    double        ln_likelihood_peak;
    double        ln_likelihood_chain;
    double *      ln_likelihood_DS;
    double *      ln_likelihood_DS_best;
    double *      ln_likelihood_DS_peak;
    int *         n_DoF_DS;
    int *         n_DoF_DS_best;
    int *         n_DoF_DS_peak;
    int           n_DoF;
    int           n_DoF_best;
    int           n_DoF_peak;
    int           P_name_length;
    char          P_name_format[8];
    int           seed;
    RNG_info *    RNG;
    double *      P;
    MCMC_DS_info *DS;
    MCMC_DS_info *last;
    // This stuff is used when MCMC_MODE_MINIMIZE_IO is on
    char *  flag_success_buffer;
    double *ln_likelihood_new_buffer;
    double *P_new_buffer;
    double *M_new_buffer;
};

#ifdef __cplusplus
extern "C" {
#endif
void init_MCMC(MCMC_info * MCMC,
               const char *problem_name,
               void *      params,
               int (*f)(double *, MCMC_info *, double **),
               int     n_P,
               double *P_init,
               char ** P_names,
               double *P_limit_min,
               double *P_limit_max,
               int     n_arrays,
               ...);
void init_MCMC_DS(MCMC_DS_info **new_DS, const char *name, int n_D, int *D, double *DS, double *dDS, void *params, int n_arrays, va_list vargs);
void start_loop_MCMC(MCMC_info *MCMC);
void end_loop_MCMC(MCMC_info *MCMC);
void restart_MCMC(MCMC_info *MCMC);
void init_MCMC_arrays(MCMC_info *MCMC);
void free_MCMC_DS(MCMC_info *MCMC);
void free_MCMC_arrays(MCMC_info *MCMC);
void free_MCMC_covariance(MCMC_info *MCMC);
void free_MCMC(MCMC_info *MCMC);
void add_MCMC_DS(MCMC_info *MCMC, const char *name, int n_D, int *D, double *DS, double *dDS, void *params, int n_arrays, ...);
void autotune_MCMC(MCMC_info *MCMC);
void autotune_MCMC_randomize(MCMC_info *MCMC);
void autotune_MCMC_temperature(MCMC_info *MCMC);
void autotune_MCMC_covariance(MCMC_info *MCMC);
void compute_MCMC_ln_likelihood_default(MCMC_info *MCMC,
                                        double **  M,
                                        double *   P,
                                        double *   ln_likeliood_DS,
                                        int *      n_DoF_DS,
                                        double *   ln_likeliood_all,
                                        int *      n_DoF_all);
void compute_MCMC_chain_stats(double **x,
                              int      n_x,
                              int      n_avg,
                              double * x_min,
                              double * x_bar_in,
                              double * x_max,
                              double * x_sigma,
                              double **auto_cor,
                              double * slopes,
                              double * dP_sub,
                              double * ln_likelihood_in,
                              double * ln_likelihood_min,
                              double * ln_likelihood_avg,
                              double * ln_likelihood_max);
void compute_MCMC(MCMC_info *MCMC);
void analyze_MCMC(MCMC_info *MCMC);
void set_MCMC_covariance(MCMC_info *MCMC, double *V);
void set_MCMC_temperature(MCMC_info *MCMC, double temperature);
void set_MCMC_iterations(MCMC_info *MCMC, int n_avg, int n_thin, int n_burn, int n_integrate);
void set_MCMC_coverage_size(MCMC_info *MCMC, int coverage_size);
void set_MCMC_likelihood_function(MCMC_info *MCMC, void (*likelihood_function)(MCMC_info *, double **, double *, double *, int *, double *, int *));
void set_MCMC_directory(MCMC_info *MCMC, const char *directory);
void set_MCMC_mode(MCMC_info *MCMC, int mode);
void set_MCMC_autotune(MCMC_info *MCMC,
                       double     success_target,
                       double     success_threshold,
                       double     covariance_threshold,
                       int        n_autotune,
                       int        n_autotune_randomize,
                       int        n_autotune_temperature,
                       int        n_autotune_covariance);
void read_MCMC_covariance(MCMC_info *MCMC, char *filename);
void read_MCMC_state(MCMC_info *MCMC);
void read_MCMC_configuration(MCMC_info *MCMC, char *filename_output_dir, int chain);
void rm_MCMC_directory(MCMC_info *MCMC);
int  generate_MCMC_chain(MCMC_info *MCMC);
void generate_MCMC_proposition(MCMC_info *MCMC, int flag_chain_init);
void generate_MCMC_parameters(MCMC_info *MCMC);
#ifdef __cplusplus
}
#endif

#endif

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/*
  C code for actionAngle calculations
*/
#ifndef __GALPY_ACTIONANGLE_H__
#define __GALPY_ACTIONANGLE_H__
#include <gsl/gsl_roots.h>
#include <gsl/gsl_spline.h>
#include "interp_2d.h"
/*
  Macro for dealing with potentially unused variables due to OpenMP
 */
/* If we're not using GNU C, elide __attribute__ if it doesn't exist*/
#ifndef __has_attribute      // Compatibility with non-clang compilers.
#define __has_attribute(x) 0
#endif
#if defined(__GNUC__) || __has_attribute(unused)
#  define UNUSED __attribute__((unused))
#else
#  define UNUSED /*NOTHING*/
#endif
/*
  Structure declarations
*/
struct actionAngleArg{ //I think this isn't used JB 06/24/14
  double (*potentialEval)(double R, double Z, double phi, double t,
			  int nargs, double * args);
  int nargs;
  double * args;
  interp_2d * i2d;
  gsl_interp_accel * acc;
};
struct pragmasolver{
  gsl_root_fsolver *s;
};
/*
  Function declarations
*/
double evaluatePotentials(double,double,int, struct potentialArg *);
void parse_actionAngleArgs(int,struct potentialArg *,int *,double *);
#endif /* actionAngle.h */

max_stars_repo_name: fardal/galpy

max_issues_repo_name: fardal/galpy

max_forks_repo_name: fardal/galpy


text:
#pragma once

#include "ShaderCompiler.h"
#include "BgfxCallback.h"

#include <Babylon/JsRuntime.h>
#include <Babylon/JsRuntimeScheduler.h>

#include <GraphicsImpl.h>

#include <NativeWindow.h>

#include <napi/napi.h>

#include <bgfx/bgfx.h>
#include <bgfx/platform.h>
#include <bimg/bimg.h>
#include <bx/allocator.h>

#include <gsl/gsl>

#include <assert.h>

#include <arcana/containers/weak_table.h>
#include <arcana/threading/cancellation.h>
#include <unordered_map>

namespace Babylon
{
    class ClearState final
    {
    public:
        void UpdateColor(float r, float g, float b, float a)
        {
            const bool needToUpdate = r != Red || g != Green || b != Blue || a != Alpha;
            if (needToUpdate)
            {
                Red = r;
                Green = g;
                Blue = b;
                Alpha = a;
                Update();
            }
        }

        void UpdateFlags(const Napi::CallbackInfo& info)
        {
            const auto flags = static_cast<uint16_t>(info[0].As<Napi::Number>().Uint32Value());
            Flags = flags;
            Update();
        }

        void UpdateDepth(const Napi::CallbackInfo& info)
        {
            const auto depth = info[0].As<Napi::Number>().FloatValue();
            const bool needToUpdate = Depth != depth;
            if (needToUpdate)
            {
                Depth = depth;
                Update();
            }
        }

        void UpdateStencil(const Napi::CallbackInfo& info)
        {
            const auto stencil = static_cast<uint8_t>(info[0].As<Napi::Number>().Int32Value());
            const bool needToUpdate = Stencil != stencil;
            if (needToUpdate)
            {
                Stencil = stencil;
                Update();
            }
        }

        arcana::weak_table<std::function<void()>>::ticket AddUpdateCallback(std::function<void()> callback)
        {
            return m_callbacks.insert(std::move(callback));
        }

        void Update()
        {
            m_callbacks.apply_to_all([](std::function<void()>& callback) {
                callback();
            });
        }

        uint32_t Color() const
        {
            uint32_t color = 0x0;
            color += static_cast<uint8_t>(Red * std::numeric_limits<uint8_t>::max());
            color = color << 8;
            color += static_cast<uint8_t>(Green * std::numeric_limits<uint8_t>::max());
            color = color << 8;
            color += static_cast<uint8_t>(Blue * std::numeric_limits<uint8_t>::max());
            color = color << 8;
            color += static_cast<uint8_t>(Alpha * std::numeric_limits<uint8_t>::max());
            return color;
        }

        float Red{68.f / 255.f};
        float Green{51.f / 255.f};
        float Blue{85.f / 255.f};
        float Alpha{1.f};
        float Depth{1.f};
        uint16_t Flags{BGFX_CLEAR_COLOR | BGFX_CLEAR_DEPTH};
        uint8_t Stencil{0};

    private:
        arcana::weak_table<std::function<void()>> m_callbacks{};
    };

    class ViewClearState final
    {
    public:
        ViewClearState(uint16_t viewId, ClearState& clearState)
            : m_viewId{viewId}
            , m_clearState{clearState}
            , m_callbackTicket{m_clearState.AddUpdateCallback([this]() { Update(); })}
        {
        }

        void UpdateColor(float r, float g, float b, float a = 1.f)
        {
            m_clearState.UpdateColor(r, g, b, a);
        }

        void UpdateColor(const Napi::CallbackInfo& info)
        {
            const auto r = info[0].As<Napi::Number>().FloatValue();
            const auto g = info[1].As<Napi::Number>().FloatValue();
            const auto b = info[2].As<Napi::Number>().FloatValue();
            const auto a = info[3].IsUndefined() ? 1.f : info[3].As<Napi::Number>().FloatValue();
            m_clearState.UpdateColor(r, g, b, a);
        }

        void UpdateFlags(const Napi::CallbackInfo& info)
        {
            m_clearState.UpdateFlags(info);
        }

        void UpdateDepth(const Napi::CallbackInfo& info)
        {
            m_clearState.UpdateDepth(info);
        }

        void UpdateStencil(const Napi::CallbackInfo& info)
        {
            m_clearState.UpdateStencil(info);
        }

        void UpdateViewId(uint16_t viewId)
        {
            m_viewId = viewId;
            Update();
        }

    private:

        void Update() const
        {
            bgfx::setViewClear(m_viewId, m_clearState.Flags, m_clearState.Color(), m_clearState.Depth, m_clearState.Stencil);
            // discard any previous set state
            bgfx::discard();
        }

        uint16_t m_viewId{};
        ClearState& m_clearState;
        arcana::weak_table<std::function<void()>>::ticket m_callbackTicket;
    };

    struct FrameBufferData final
    {
    private:
        std::unique_ptr<ClearState> m_clearState{};

    public:
        FrameBufferData(bgfx::FrameBufferHandle frameBuffer, uint16_t viewId, uint16_t width, uint16_t height, bool actAsBackBuffer = false)
            : m_clearState{std::make_unique<ClearState>()}
            , FrameBuffer{frameBuffer}
            , ViewId{viewId}
            , ViewClearState{ViewId, *m_clearState}
            , Width{width}
            , Height{height}
            , ActAsBackBuffer{actAsBackBuffer}
        {
            assert(ViewId < bgfx::getCaps()->limits.maxViews);
        }

        FrameBufferData(bgfx::FrameBufferHandle frameBuffer, uint16_t viewId, ClearState& clearState, uint16_t width, uint16_t height, bool actAsBackBuffer = false)
            : m_clearState{}
            , FrameBuffer{frameBuffer}
            , ViewId{viewId}
            , ViewClearState{ViewId, clearState}
            , Width{width}
            , Height{height}
            , ActAsBackBuffer{actAsBackBuffer}
        {
            assert(ViewId < bgfx::getCaps()->limits.maxViews);
        }

        FrameBufferData(FrameBufferData&) = delete;

        ~FrameBufferData()
        {
            bgfx::destroy(FrameBuffer);
        }

        void UseViewId(uint16_t viewId)
        {
            ViewId = viewId;
            ViewClearState.UpdateViewId(ViewId);
        }

        void SetUpView(uint16_t viewId)
        {
            bgfx::setViewFrameBuffer(viewId, FrameBuffer);
            UseViewId(viewId);
            bgfx::setViewRect(ViewId, 0, 0, Width, Height);
        }

        bgfx::FrameBufferHandle FrameBuffer{bgfx::kInvalidHandle};
        bgfx::ViewId ViewId{};
        Babylon::ViewClearState ViewClearState;
        uint16_t Width{};
        uint16_t Height{};
        // When a FrameBuffer acts as a back buffer, it means it will not be used as a texture in a shader.
        // For example as a post process. It will be used as-is in a swapchain or for direct rendering (XR)
        // When this flag is true, projection matrix will not be flipped for API that would normaly need it.
        // Namely Direct3D and Metal.
        bool ActAsBackBuffer{false};
    };

    struct FrameBufferManager final
    {
        FrameBufferManager()
        {
            m_boundFrameBuffer = m_backBuffer = new FrameBufferData(BGFX_INVALID_HANDLE, GetNewViewId(), bgfx::getStats()->width, bgfx::getStats()->height);
        }

        FrameBufferData* CreateNew(bgfx::FrameBufferHandle frameBufferHandle, uint16_t width, uint16_t height)
        {
            return new FrameBufferData(frameBufferHandle, GetNewViewId(), width, height);
        }

        FrameBufferData* CreateNew(bgfx::FrameBufferHandle frameBufferHandle, ClearState& clearState, uint16_t width, uint16_t height, bool actAsBackBuffer)
        {
            return new FrameBufferData(frameBufferHandle, GetNewViewId(), clearState, width, height, actAsBackBuffer);
        }

        void Bind(FrameBufferData* data)
        {
            m_boundFrameBuffer = data;

            // TODO: Consider doing this only on bgfx::reset(); the effects of this call don't survive reset, but as
            // long as there's no reset this doesn't technically need to be called every time the frame buffer is bound.
            m_boundFrameBuffer->SetUpView(GetNewViewId());

            // bgfx::setTexture()? Why?
            // TODO: View order?
            m_renderingToTarget = !m_boundFrameBuffer->ActAsBackBuffer;
        }

        FrameBufferData& GetBound() const
        {
            return *m_boundFrameBuffer;
        }

        void Unbind(FrameBufferData* data)
        {
            // this assert is commented because of an issue with XR described here : https://github.com/BabylonJS/BabylonNative/issues/344
            //assert(m_boundFrameBuffer == data);
            (void)data;
            m_boundFrameBuffer = m_backBuffer;
            m_renderingToTarget = false;
        }

        uint16_t GetNewViewId()
        {
            m_nextId++;
            assert(m_nextId < bgfx::getCaps()->limits.maxViews);
            return m_nextId;
        }

        void Reset()
        {
            m_nextId = 0;
        }

        bool IsRenderingToTarget() const
        {
            return m_renderingToTarget;
        }

    private:
        FrameBufferData* m_boundFrameBuffer{nullptr};
        FrameBufferData* m_backBuffer{nullptr};
        uint16_t m_nextId{0};
        bool m_renderingToTarget{false};
    };

    struct TextureData final
    {
        ~TextureData()
        {
            if (bgfx::isValid(Handle))
            {
                bgfx::destroy(Handle);
            }
        }

        bgfx::TextureHandle Handle{bgfx::kInvalidHandle};
        uint32_t Width{0};
        uint32_t Height{0};
        uint32_t Flags{0};
        uint8_t AnisotropicLevel{0};
    };

    struct ImageData final
    {
        ~ImageData()
        {
            if (Image)
            {
                bimg::imageFree(Image.get());
            }
        }
        std::unique_ptr<bimg::ImageContainer> Image;
    };

    struct UniformInfo final
    {
        uint8_t Stage{};
        bgfx::UniformHandle Handle{bgfx::kInvalidHandle};
        bool YFlip{false};
    };

    struct ProgramData final
    {
        ProgramData() = default;
        ProgramData(const ProgramData&) = delete;
        ProgramData(ProgramData&&) = delete;

        ~ProgramData()
        {
            bgfx::destroy(Program);
        }

        std::unordered_map<std::string, uint32_t> VertexAttributeLocations{};
        std::unordered_map<std::string, UniformInfo> VertexUniformInfos{};
        std::unordered_map<std::string, UniformInfo> FragmentUniformInfos{};

        bgfx::ProgramHandle Program{};

        struct UniformValue
        {
            std::vector<float> Data{};
            uint16_t ElementLength{};
            bool YFlip{false};
        };

        std::unordered_map<uint16_t, UniformValue> Uniforms{};

        void SetUniform(bgfx::UniformHandle handle, gsl::span<const float> data, bool YFlip, size_t elementLength = 1)
        {
            UniformValue& value = Uniforms[handle.idx];
            value.Data.assign(data.begin(), data.end());
            value.ElementLength = static_cast<uint16_t>(elementLength);
            value.YFlip = YFlip;
        }
    };

    class IndexBufferData;
    class VertexBufferData;

    struct VertexArray final
    {
        ~VertexArray()
        {
            for (auto& vertexBuffer : vertexBuffers)
            {
                bgfx::destroy(vertexBuffer.vertexLayoutHandle);
            }
        }

        struct IndexBuffer
        {
            const IndexBufferData* data{};
        };

        IndexBuffer indexBuffer{};

        struct VertexBuffer
        {
            const VertexBufferData* data{};
            uint32_t startVertex{};
            bgfx::VertexLayoutHandle vertexLayoutHandle{};
        };

        std::vector<VertexBuffer> vertexBuffers{};
    };

    class NativeEngine final : public Napi::ObjectWrap<NativeEngine>
    {
        static constexpr auto JS_CLASS_NAME = "_NativeEngine";
        static constexpr auto JS_ENGINE_CONSTRUCTOR_NAME = "Engine";
        static constexpr auto JS_AUTO_RENDER_PROPERTY_NAME = "_AUTO_RENDER";

    public:
        NativeEngine(const Napi::CallbackInfo& info);
        NativeEngine(const Napi::CallbackInfo& info, JsRuntime& runtime, Plugins::Internal::NativeWindow& nativeWindow);
        ~NativeEngine();

        static void Initialize(Napi::Env, bool autoRender);

        FrameBufferManager& GetFrameBufferManager();
        void Dispatch(std::function<void()>);

        void ScheduleRender();

        const bool AutomaticRenderingEnabled{};
        JsRuntimeScheduler RuntimeScheduler;

    private:
        void Dispose();

        void Dispose(const Napi::CallbackInfo& info);
        Napi::Value GetEngine(const Napi::CallbackInfo& info); // TODO: Hack, temporary method. Remove as part of the change to get rid of NapiBridge.
        void RequestAnimationFrame(const Napi::CallbackInfo& info);
        Napi::Value CreateVertexArray(const Napi::CallbackInfo& info);
        void DeleteVertexArray(const Napi::CallbackInfo& info);
        void BindVertexArray(const Napi::CallbackInfo& info);
        Napi::Value CreateIndexBuffer(const Napi::CallbackInfo& info);
        void DeleteIndexBuffer(const Napi::CallbackInfo& info);
        void RecordIndexBuffer(const Napi::CallbackInfo& info);
        void UpdateDynamicIndexBuffer(const Napi::CallbackInfo& info);
        Napi::Value CreateVertexBuffer(const Napi::CallbackInfo& info);
        void DeleteVertexBuffer(const Napi::CallbackInfo& info);
        void RecordVertexBuffer(const Napi::CallbackInfo& info);
        void UpdateDynamicVertexBuffer(const Napi::CallbackInfo& info);
        Napi::Value CreateProgram(const Napi::CallbackInfo& info);
        Napi::Value GetUniforms(const Napi::CallbackInfo& info);
        Napi::Value GetAttributes(const Napi::CallbackInfo& info);
        void SetProgram(const Napi::CallbackInfo& info);
        void SetState(const Napi::CallbackInfo& info);
        void SetZOffset(const Napi::CallbackInfo& info);
        Napi::Value GetZOffset(const Napi::CallbackInfo& info);
        void SetDepthTest(const Napi::CallbackInfo& info);
        Napi::Value GetDepthWrite(const Napi::CallbackInfo& info);
        void SetDepthWrite(const Napi::CallbackInfo& info);
        void SetColorWrite(const Napi::CallbackInfo& info);
        void SetBlendMode(const Napi::CallbackInfo& info);
        void SetMatrix(const Napi::CallbackInfo& info);
        void SetInt(const Napi::CallbackInfo& info);
        void SetIntArray(const Napi::CallbackInfo& info);
        void SetIntArray2(const Napi::CallbackInfo& info);
        void SetIntArray3(const Napi::CallbackInfo& info);
        void SetIntArray4(const Napi::CallbackInfo& info);
        void SetFloatArray(const Napi::CallbackInfo& info);
        void SetFloatArray2(const Napi::CallbackInfo& info);
        void SetFloatArray3(const Napi::CallbackInfo& info);
        void SetFloatArray4(const Napi::CallbackInfo& info);
        void SetMatrices(const Napi::CallbackInfo& info);
        void SetMatrix3x3(const Napi::CallbackInfo& info);
        void SetMatrix2x2(const Napi::CallbackInfo& info);
        void SetFloat(const Napi::CallbackInfo& info);
        void SetFloat2(const Napi::CallbackInfo& info);
        void SetFloat3(const Napi::CallbackInfo& info);
        void SetFloat4(const Napi::CallbackInfo& info);
        Napi::Value CreateTexture(const Napi::CallbackInfo& info);
        Napi::Value CreateDepthTexture(const Napi::CallbackInfo& info);
        void LoadTexture(const Napi::CallbackInfo& info);
        void LoadCubeTexture(const Napi::CallbackInfo& info);
        void LoadCubeTextureWithMips(const Napi::CallbackInfo& info);
        Napi::Value GetTextureWidth(const Napi::CallbackInfo& info);
        Napi::Value GetTextureHeight(const Napi::CallbackInfo& info);
        void SetTextureSampling(const Napi::CallbackInfo& info);
        void SetTextureWrapMode(const Napi::CallbackInfo& info);
        void SetTextureAnisotropicLevel(const Napi::CallbackInfo& info);
        void SetTexture(const Napi::CallbackInfo& info);
        void DeleteTexture(const Napi::CallbackInfo& info);
        Napi::Value CreateFrameBuffer(const Napi::CallbackInfo& info);
        void DeleteFrameBuffer(const Napi::CallbackInfo& info);
        void BindFrameBuffer(const Napi::CallbackInfo& info);
        void UnbindFrameBuffer(const Napi::CallbackInfo& info);
        void DrawIndexed(const Napi::CallbackInfo& info);
        void Draw(const Napi::CallbackInfo& info);
        void Clear(const Napi::CallbackInfo& info);
        void ClearColor(const Napi::CallbackInfo& info);
        void ClearStencil(const Napi::CallbackInfo& info);
        void ClearDepth(const Napi::CallbackInfo& info);
        Napi::Value GetRenderWidth(const Napi::CallbackInfo& info);
        Napi::Value GetRenderHeight(const Napi::CallbackInfo& info);
        void SetViewPort(const Napi::CallbackInfo& info);
        void GetFramebufferData(const Napi::CallbackInfo& info);
        Napi::Value GetRenderAPI(const Napi::CallbackInfo& info);

        void UpdateSize(size_t width, size_t height);

        template<typename SchedulerT>
        arcana::task<void, std::exception_ptr> GetRequestAnimationFrameTask(SchedulerT&);

        bool m_isRenderScheduled{false};

        arcana::cancellation_source m_cancelSource{};

        ShaderCompiler m_shaderCompiler;

        ProgramData* m_currentProgram{nullptr};
        arcana::weak_table<std::unique_ptr<ProgramData>> m_programDataCollection{};

        JsRuntime& m_runtime;
        Graphics::Impl& m_graphicsImpl;

        bx::DefaultAllocator m_allocator;
        uint64_t m_engineState;

        FrameBufferManager m_frameBufferManager{};

        Plugins::Internal::NativeWindow::NativeWindow::OnResizeCallbackTicket m_resizeCallbackTicket;

        template<int size, typename arrayType>
        void SetTypeArrayN(const Napi::CallbackInfo& info);

        template<int size>
        void SetFloatN(const Napi::CallbackInfo& info);

        template<int size>
        void SetMatrixN(const Napi::CallbackInfo& info);

        // Scratch vector used for data alignment.
        std::vector<float> m_scratch{};

        Napi::FunctionReference m_requestAnimationFrameCallback{};

        // webgl/opengl draw call parameters allow to set first index and number of indices used for that call
        // but with bgfx, those parameters must be set when binding the index buffer
        // at the time of webgl binding, we don't know those values yet
        // so a pointer to the to-bind buffer is kept and the buffer is bound to bgfx at the time of the drawcall
        const IndexBufferData* m_currentBoundIndexBuffer{};
    };
}

max_stars_repo_name: syntheticmagus/BabylonNative

max_issues_repo_name: syntheticmagus/BabylonNative

max_forks_repo_name: syntheticmagus/BabylonNative