text:
/* StdVarTrajs_Sampling.c

Simulation for calculating allele trajectories, adaptation from standing variation
To be fed into MSMS

Simulation uses routines found with the GNU Scientific Library (GSL)
(http://www.gnu.org/software/gsl/)
Since GSL is distributed under the GNU General Public License
(http://www.gnu.org/copyleft/gpl.html), you must download it
separately from this file.

This program can be compiled with e.g. GCC using a command like:
gcc StdVarTrajs_Sampling -lm -lgsl -lgslcblas -I/usr/local/include -L/usr/local/lib StdVarTrajs_Sampling.c

Then run by executing:
./StdVarTrajs_Sampling N s h x0 reps cutoff
Where:
- N is the (haploid) size
- s is the selection coefficient of the beneficial allele
- h is dominance of the beneficial allele
- x0 is frequency of the allele when it started to be selected for
- reps is how many times the second allele should FIX before simulation stops
(the number of actual runs is greater due to stochastic loss)
- cutoff denotes how high the derived allele can be in the neutral phase ('0' for fixation, '1' for x0 + 1/2N)
*/

/* Preprocessor statements */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stddef.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <sys/stat.h>
#include <sys/types.h>

#define INITTS 310000

/* Function prototypes */
void geninit(double *geninit, double pee);
void selection(double *geninit, double h, double s);
void reproduction(double *geninit);
void bentraj(double *XFreq, unsigned int *btime, double *genotype, double x0, unsigned int rep, unsigned int N, double s, double h, const gsl_rng *r);
void neutraj(double *NFreq, unsigned int *ntime, double *genotype, double x0, unsigned int rep, unsigned int N, unsigned int cutoff, const gsl_rng *r);
void trajprint(double **AFreq, unsigned int *ttime, unsigned int rep);
void Wait();

/* Initialising genotypes */
void geninit(double *geninit, double pee){

	*(geninit + 0) = (1-pee)*(1-pee);
	*(geninit + 1) = 2*pee*(1-pee);
	*(geninit + 2) = pee*pee;

}	/* End of gen initiation routine */

/* Selection routine (single locus) */
void selection(double *geninit, double h, double s){
	/* Fitness of each genotype */
	double WAA, WAa, Waa;
	double Wmean;				/* Mean fitness */

	Waa = 1;
	WAa = 1 + h*s;
	WAA = 1 + s;

	/* Mean fitness calculation */
	Wmean = ((*(geninit + 0))*Waa) + ((*(geninit + 1))*WAa) + ((*(geninit + 2))*WAA);

	/* Changing frequencies by selection */
	*(geninit + 0) = ((*(geninit + 0))*Waa)/Wmean;
	*(geninit + 1) = ((*(geninit + 1))*WAa)/Wmean;
	*(geninit + 2) = ((*(geninit + 2))*WAA)/Wmean;

}	/* End of selection routine */

/* Reproduction routine */
void reproduction(double *geninit){

	/* Fed-in genotype frequencies (for ease of programming) */
	double gaas, gAas, gAAs;
	/* Haplotypes */
	double xA, xa;

	/* Initial definition of genotypes */
	gaas = *(geninit + 0);
	gAas = *(geninit + 1);
	gAAs = *(geninit + 2);

	/* Baseline change in haplotype frequencies with obligate outcrossing */
	xA = gAAs + (gAas)/2.0;
	xa = gaas + (gAas)/2.0;

	/* Change in frequencies (HWE) */
	*(geninit + 0) = pow(xa,2);
	*(geninit + 1) = 2.0*xA*xa;
	*(geninit + 2) = pow(xA,2);

}	/* End of reproduction routine */

/* Routine to simulate beneficial allele trajectory */
void bentraj(double *XFreq, unsigned int *btime, double *genotype, double x0, unsigned int rep, unsigned int N, double s, double h, const gsl_rng *r){

	unsigned int gt = 1;		/* Generation time */
	unsigned int done = 0;		/* Counter if routine done or not */
	unsigned int exitl = 0;		/* Decide if exit loop or not */
	unsigned int i;				/* Genotype counter */
	double x = 0;				/* Beneficial allele frequency */

	unsigned int *gensamp = calloc(3,sizeof(unsigned int));		/* New population samples */

	while(done == 0){
		exitl = 0;
		gt = 1;
		x = x0;
		*(XFreq + 0) = x0;
		/* Initiating genotypes */
		geninit(genotype,x0);
		while(exitl == 0){

			/* Selection routine */
    		selection(genotype,h,s);

			/* Reproduction routine */
    		reproduction(genotype);

    		/* Sampling based on new frequencies */
	       	gsl_ran_multinomial(r,3,N,genotype,gensamp);
       		for(i = 0; i < 3; i++){
    	   		*(genotype + i) = (*(gensamp + i))/(1.0*N);
	       	}

	       	x = 0.5*(*(genotype + 1)) + (*(genotype + 2));
			*(XFreq + gt) = x;
			gt += 1;
			if(gt > INITTS){
				fprintf(stderr,"Number of generations in selection trajectory exceed vector length (INITTS).\n");
				exit(1);
			}

			if( (x == 1) || (x == 0) ){
				exitl = 1;
			}
		}

		if(x == 1){
			done = 1;
			*(btime + rep) = (gt - 1);
		}
	}

	free(gensamp);

}

/* Routine to simulate neutral allele trajectory */
void neutraj(double *NFreq, unsigned int *ntime, double *genotype, double x0, unsigned int rep, unsigned int N, unsigned int cutoff, const gsl_rng *r){

	unsigned int gt = 1;		/* Generation time */
	unsigned int done = 0;		/* Counter if routine done or not */
	unsigned int exitl = 0;		/* Decide if exit loop or not */
	unsigned int i;				/* Genotype counter */
	double x = 0;				/* Neutral allele frequency */
	double xmax = 1;			/* Value of x to use as maximum cutoff */

	unsigned int *gensamp = calloc(3,sizeof(unsigned int));		/* New population samples */

	if(cutoff == 0){
		xmax = 1;
	}else if(cutoff == 1){
		xmax = x0 + (1.0/(2.0*N));
		if(xmax > 1){
			xmax = 1;
		}
	}
/*	printf("xmax is %lf\n",xmax);	*/

	while(done == 0){
		exitl = 0;
		gt = 0;		/* Note start from zero here since x0 case defined in selection trajectory */
		x = x0;

		/* Initiating genotypes */
		geninit(genotype,x0);

		while(exitl==0){

			/* Reproduction routine */
    		reproduction(genotype);

			/* Random sampling of allele frequencies (drift) */
	       	gsl_ran_multinomial(r,3,N,genotype,gensamp);
       		for(i = 0; i < 3; i++){
    	   		*(genotype + i) = (*(gensamp + i))/(1.0*N);
	       	}
	       	x = 0.5*(*(genotype + 1)) + (*(genotype + 2));
			*(NFreq + gt) = x;

			gt++;
			if(gt > INITTS){
				printf("Number of generations in neutral trajectory exceed vector length (INITTS).\n");
				exit(1);
			}

			if( (x == 0) || (x >= xmax) ){
				exitl = 1;
			}
		}

		if(x == 0){
			done = 1;
			*(ntime + rep) = (gt - 1);
		}
	}

	free(gensamp);

}

/* Print out trajectories to file */
void trajprint(double **AFreq, unsigned int *ttime, unsigned int rep){
	unsigned int j;
	char filename[32];
	FILE *ofp_tr;				/* Pointer for file output */

	sprintf(filename,"Traj/ATraj%d.dat",(rep + 1));
	ofp_tr = fopen(filename,"w");

	for(j = 0; j < *(ttime + rep); j++){
		fprintf(ofp_tr,"%lf %lf %lf\n",*((*(AFreq + 0)) + j),*((*(AFreq + 1)) + j),*((*(AFreq + 2)) + j));
	}

	fclose(ofp_tr);

}

void Wait(){
	printf("Press Enter to Continue");
	while( getchar() != '\n' );
	printf("\n");
}

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

	/* Declare variables here */
	unsigned int i, a;			/* Rep counter, memory counter */
	unsigned int N = 0;			/* Population Size */
	unsigned int Nreps = 0;		/* Number of repetitions */
	unsigned int acc = 0;		/* Accumulator (when merging two trajectories) */
	unsigned int es = 10;		/* How many extra steps to put in the trajectory files */
	unsigned int cutoff = 0;	/* What kind of back-in-time cutoff to use */
	int j;						/* Timestep counter */
	double s = 0;				/* Strength of selection */
	double h = 0;				/* Dominance level */
	double x0 = 0;				/* Initial allele frequency */
	double lim = 0;				/* Limit of acceptable allele frequencies */
	double dx = 0;				/* Rescaled timestep so can be used by MSMS */
	FILE *ofp_sd;				/* Pointer for seed output */

	/* GSL random number definitions */
	const gsl_rng_type * T;
	gsl_rng * r;

	/* Reading in data from command line */
	/* <Program> N s h x0 NReps cutoff */
	if(argc != 7){
		fprintf(stderr,"Six inputs are needed (N s h x0 Reps cutoff).\n");
		exit(1);
	}

	N = atoi(argv[1]);
	if(N <= 0){
		fprintf(stderr,"Total Population size N is zero or negative, not allowed.\n");
		exit(1);
	}

	/* Defining min frequency limit ** based on unscaled N ** */
	lim = (1.0/(1.0*N));

	/* Defining timestep, dx = 1/4N since that's what MSMS uses */
	dx = (1.0/(4.0*N));

	s = strtod(argv[2],NULL);
	if(s < 0){
		fprintf(stderr,"Allele strength s is negative, not allowed.\n");
		exit(1);
	}

	h = strtod(argv[3],NULL);
	if(h < 0 || h > 1){
		fprintf(stderr,"Dominance value must lie between 0 and 1.\n");
		exit(1);
	}

	x0 = strtod(argv[4],NULL);
	if(x0 < lim || x0 > (1-lim)){
		fprintf(stderr,"Initial mutant frequency must lie between %0.5lf and %0.5lf.\n",lim,1-lim);
		exit(1);
	}

	/* Number of samples/reps to take */
	Nreps = atoi(argv[5]);
	if(Nreps <= 0){
		fprintf(stderr,"Must set positive number of repetitions.\n");
		exit(1);
	}

	/* Determining cutoff */
	cutoff = atoi(argv[6]);
	if(!(cutoff == 0 || cutoff == 1)){
		fprintf(stderr,"'Cutoff' switch must equal 0 or 1.\n");
		exit(1);
	}

	unsigned int *btime = calloc(Nreps,sizeof(unsigned int));			/* Number of generations needed for ben allele fixation */
	unsigned int *ntime = calloc(Nreps,sizeof(unsigned int));			/* Number of generations needed for neut allele loss */
	unsigned int *ttime = calloc(Nreps,sizeof(unsigned int));			/* Number of generations needed for whole processes */
	double *genotype = calloc(3,sizeof(double));						/* Genotype frequencies */

	/* create a generator chosen by the
    environment variable GSL_RNG_TYPE */

	gsl_rng_env_setup();
	if (!getenv("GSL_RNG_SEED")) gsl_rng_default_seed = time(0);
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	ofp_sd = fopen("Seed.dat","w");
	fprintf(ofp_sd,"%lu\n",gsl_rng_default_seed);
	fclose(ofp_sd);

	mkdir("Traj/", 0777);

	/* Executing simulation */
	for(i = 0; i < Nreps; i++){

		if(Nreps > 100){
			if(i%100 == 0){
				printf("Starting run %d\n",i);
			}
		}

		double *XFreq = calloc(INITTS,sizeof(double));			/* Beneficial allele frequency */
		double *NFreq = calloc(INITTS,sizeof(double));			/* Neutral allele frequency */

		/* Simulating beneficial trajectories */
		bentraj(XFreq, btime, genotype, x0, i, N, s, h, r);

		/* Simulation neutral trajectory */
		neutraj(NFreq, ntime, genotype, x0, i, N, cutoff, r);

		*(ttime + i) = (*(btime+i)) + (*(ntime+i) + 2 + es);
/*		printf("%d %d %d\n",*(btime+i),*(ntime+i),*(ttime+i));	*/

		/* Defining memory for AFreq */
		double **AFreq = calloc(3,sizeof(double *));					/* Allele frequency over time */
		for(a = 0; a < 3; a++){											/* Assigning space for each run */
			AFreq[a] = calloc(*(ttime + i),sizeof(double));
		}

		/* Merging trajectories
		Timing is set up so tfix = 0
		And discrete generations is scaled by timestep, dx
		*/

		/* First add ten timesteps for x = 0
		So trajectory fully covered in MSMS */
		acc = 0;
		for(j = 0; j < es; j++){
			*((*(AFreq + 0)) + j) = ((*(ttime + i) - es - 1)*dx) + (es - j);
			*((*(AFreq + 1)) + j) = 1.0;
			*((*(AFreq + 2)) + j) = 0.0;
		}
		for(j = *(ntime + i); j >= 0; j--){
			*((*(AFreq + 0)) + acc + es) = (*(ttime + i) - acc - es - 1)*dx;
			*((*(AFreq + 1)) + acc + es) = 1.0-(*(NFreq + j));
			*((*(AFreq + 2)) + acc + es) = *(NFreq + j);
			acc++;
		}
		for(j = 0; j <= *(btime + i); j++){
			*((*(AFreq + 0)) + acc + es) = (*(ttime + i) - acc - es - 1)*dx;
			*((*(AFreq + 1)) + acc + es) = 1.0-(*(XFreq + j));
			*((*(AFreq + 2)) + acc + es) = *(XFreq + j);
			acc++;
		}
		/*
		*(ttime + i) = (acc-1);
		printf("%d\n",*(ttime+i));
		*/

		free(NFreq);
		free(XFreq);

		/* Printout of trajectories */
		trajprint(AFreq, ttime, i);

		for(a = 0; a < 3; a++){
	 		free(AFreq[a]);
		}
		free(AFreq);

	}

	gsl_rng_free(r);
	free(genotype);
	free(ttime);
	free(ntime);
	free(btime);

	return 0;

}	/* End of main program */

/* End of File */

max_stars_repo_name: MattHartfield/SweepDomSelf

max_issues_repo_name: MattHartfield/SweepDomSelf

max_forks_repo_name: MattHartfield/SweepDomSelf


text:
/*
  #################################################################
  ##   Private MATRIX Library ver 0.3                            ##
  ##   Made by "Gun Woo Park", Email : gunwoo.park@unina.it      ##
  ##   Fellow of SUPOLEN, Marie Curie Actions                    ##
  ##   DICMaPI, University of Naples Federico II                 ##
  #################################################################
*/

// all the integer is converted to MKL_LONG

// @ denoted for tested function
// Lots of functions are based on Call by Reference.
// Notice that for matrix operation, we have to put &.
// Vector can use by rows = 1.

// Note that I seperated library for linear algebra and matrix.
// If some error occurs during compiling, please check to include linear algebra code.

#ifndef MATRIX_H
#define MATRIX_H
#define TRUE 1
#define FALSE 0
#define BIT 64
#include <iostream>
#include <fstream>
#include <mkl.h>
#include <math.h>
//#include <omp.h>

extern "C" {
#include <gsl/gsl_sort_double.h>
}

class MATRIX
{
 public:
  MKL_LONG rows;
  MKL_LONG cols;
  MKL_LONG size;
  double *data;  // 1-dimensional array
  bool VECTOR;
  bool INITIALIZATION;
  bool DIAGONALIZATION;
  double *eigen_value;
  MKL_LONG DIAGONALIZATION_INITIAL()                     //@
  {
    /* eigen_value = new double [rows]; */
    /* eigen_value = (double*) mkl_malloc(rows*sizeof(double), BIT); */
    eigen_value = new double [rows];
    DIAGONALIZATION = TRUE;
    return 0;
  }
  MKL_LONG
    print_eigen();
  MKL_LONG
    print_eigen
    (MKL_LONG n_ele);                                //@
  // Public Member Function
  MKL_LONG
    p_self();                                     //@
  MKL_LONG
    initial();
  MKL_LONG
    initial
    (MKL_LONG N_r, MKL_LONG N_c);                    //@
  MKL_LONG
    initial
    (MKL_LONG N_r, MKL_LONG N_c, double x);
  MKL_LONG
    print();                                      //@
  MKL_LONG
    print
    (MKL_LONG n_row, MKL_LONG n_col);                  //@

  MKL_LONG
    fprint_skip
    (std::ofstream& file, MKL_LONG N_skip);
  MKL_LONG
    fprint_skip_transpose
    (std::ofstream& file, MKL_LONG N_skip);
  MKL_LONG
    fprint_LONG_skip
    (std::ofstream& file, MKL_LONG N_skip);
  MKL_LONG
    fprint_LONG_skip_transpose
    (std::ofstream& file, MKL_LONG N_skip);
  MKL_LONG
    fprint_LONG_skip_transpose_LIMROWS
    (std::ofstream& file, MKL_LONG N_skip, MKL_LONG N_lim_rows);
  MKL_LONG
    fprint
    (std::ofstream& file);                             // ???
  MKL_LONG
    fprint_transpose
    (std::ofstream& file);                             // ???
  MKL_LONG
    fprint_LONG
    (std::ofstream& file);                             // ???
  MKL_LONG
    fprint_LONG_transpose
    (std::ofstream& file);                             // ???
  MKL_LONG
    fprint_row
    (std::ofstream& file, MKL_LONG given_row);
  MKL_LONG
    fprint_out_skip
    (std::ofstream& file, MKL_LONG N_skip);
  MKL_LONG
    fprint_out
    (std::ofstream& file);                             // ???

  MKL_LONG
    test_arr();                                   //@
  MKL_LONG
    test_arr_symm();                              //@
  MKL_LONG
    nonzero();
  MKL_LONG
    ABS_cond
    (double x);
  double
    ABS_sum();
  MKL_LONG
    force_clean(double TOLERANCE);
  double
    sum();
  double
    average();
  MKL_LONG
    add(MATRIX& given_MAT);

  MKL_LONG
  read_exist_data(const char* fn_given_array);

  // Constructor
  MATRIX();                                         //@
  MATRIX
    (MKL_LONG N_r, MKL_LONG N_c);                         //@
  MATRIX
    (MKL_LONG N_r, MKL_LONG N_c, double x);               //@
  // Copy Constructor
  MKL_LONG
    copy_obj
    (const MATRIX& Mat);
  MATRIX
    (const MATRIX& Mat);                        //@
  /* int connect_copy_constructor(const MATRIX& Mat); */
  // Destructor
  virtual
    ~MATRIX();                                        //@

  MKL_LONG
    data_delete();                                //@
  MKL_LONG
    undefined_error();                            //@
  /* MKL_LONG index(MKL_LONG i, MKL_LONG j);                          //@ */ // inlined
  MKL_LONG
    sort();
  MKL_LONG
    sort2
    (MATRIX& given_index);
  /* MKL_LONG sort2(size_t* index); */
  // inlined functions
  // note that function inside class definition is inlined
  // these are frequently called function that have small amount of computation

  double
    norm()
  {
    double result = 0.0;
    for(MKL_LONG  i=0; i<size; i++)
      {
        result += pow(data[i],2.0);
      }
    return sqrt(result);
  }

  MKL_LONG
    set_value
    (double x)
  {
    if(INITIALIZATION)
      {
        for(MKL_LONG  i=0; i<size; i++)
          {
            data[i] = x;
          } // i
      }
    else
      {
        undefined_error();
      }
    return 0;
  }

  MKL_LONG
    index
    (MKL_LONG i, MKL_LONG j)
  {
    return i*cols+j;
  }

  // operator
  // inlined operator overloading
  MATRIX&
    operator=
    (const MATRIX &Mat)
    {
      copy_obj(Mat);
      return *this;
    }

  MATRIX&
    operator+=
    (const MATRIX &Mat)
    {
      // std::cout << "the operator += for MATRIX class is in buggy status\n";
      // std::cout << "if this operator is called, to check the internal problems\n";
      if (size != Mat.size)
        std::cout << "The += operator must happen when two size are the same\n";
      // copy_obj(Mat);
      // for(MKL_LONG i=0; i<size; i++)
      //   {
      //     data[i] += Mat.data[i];
      //   }
      cblas_daxpy(size, 1.0, Mat.data, 1, data, 1);
      return *this;
    }

  double&
    operator()
    (MKL_LONG  i, MKL_LONG  j)
  {
    return data[i*cols+j];
  }

  double&
    operator()
    (MKL_LONG  i)
  {
    return data[i];
  }

};

inline MKL_LONG
make_unit_vector
(MATRIX& given_vec)
{
  double norm = 0.;
  for(MKL_LONG  i=0; i<given_vec.size; i++)
    {
      norm += pow(given_vec.data[i], 2.0);
    }
  norm = sqrt(norm);
  for(MKL_LONG  i=0; i<given_vec.size; i++)
    {
      given_vec.data[i] /= norm;
    }
  return norm;
}

inline MKL_LONG
matrix_mul
(MATRIX& given_vec, double val)
{
  cblas_dscal(given_vec.size, val, given_vec.data, 1);
  return 0;
}

MKL_LONG
nonzero
(const MATRIX &A);

#endif

ext: h


text:
#ifndef PITCH_H
#define PITCH_H

#include <cstring>

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

#define DECL_PITCH(name)  bool pitch_ ## name (gsl_vector *frame, double *f0res, double *T0res)

DECL_PITCH(AMDF);

#endif // PITCH_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


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

/* steffenson.c -- steffenson root finding algorithm

   This is Newton's method with an Aitken "delta-squared"
   acceleration of the iterates. This can improve the convergence on
   multiple roots where the ordinary Newton algorithm is slow.

   x[i+1] = x[i] - f(x[i]) / f'(x[i])

   x_accelerated[i] = x[i] - (x[i+1] - x[i])**2 / (x[i+2] - 2*x[i+1] + x[i])

   We can only use the accelerated estimate after three iterations,
   and use the unaccelerated value until then.

 */

#include <config.h>

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

#include "roots.h"

typedef struct
  {
    double f, df;
    double x;
    double x_1;
    double x_2;
    int count;
  }
steffenson_state_t;

static int steffenson_init (void * vstate, gsl_function_fdf * fdf, double * root);
static int steffenson_iterate (void * vstate, gsl_function_fdf * fdf, double * root);

static int
steffenson_init (void * vstate, gsl_function_fdf * fdf, double * root)
{
  steffenson_state_t * state = (steffenson_state_t *) vstate;

  const double x = *root ;

  state->f = GSL_FN_FDF_EVAL_F (fdf, x);
  state->df = GSL_FN_FDF_EVAL_DF (fdf, x) ;

  state->x = x;
  state->x_1 = 0.0;
  state->x_2 = 0.0;

  state->count = 1;

  return GSL_SUCCESS;

}

static int
steffenson_iterate (void * vstate, gsl_function_fdf * fdf, double * root)
{
  steffenson_state_t * state = (steffenson_state_t *) vstate;

  double x_new, f_new, df_new;

  double x_1 = state->x_1 ;
  double x = state->x ;

  if (state->df == 0.0)
    {
      GSL_ERROR("derivative is zero", GSL_EZERODIV);
    }

  x_new = x - (state->f / state->df);

  GSL_FN_FDF_EVAL_F_DF(fdf, x_new, &f_new, &df_new);

  state->x_2 = x_1 ;
  state->x_1 = x ;
  state->x = x_new;

  state->f = f_new ;
  state->df = df_new ;

  if (!finite (f_new))
    {
      GSL_ERROR ("function value is not finite", GSL_EBADFUNC);
    }

  if (state->count < 3)
    {
      *root = x_new ;
      state->count++ ;
    }
  else
    {
      double u = (x - x_1) ;
      double v = (x_new - 2 * x + x_1);

      if (v == 0)
        *root = x_new;  /* avoid division by zero */
      else
        *root = x_1 - u * u / v ;  /* accelerated value */
    }

  if (!finite (df_new))
    {
      GSL_ERROR ("derivative value is not finite", GSL_EBADFUNC);
    }

  return GSL_SUCCESS;
}

static const gsl_root_fdfsolver_type steffenson_type =
{"steffenson",                          /* name */
 sizeof (steffenson_state_t),
 &steffenson_init,
 &steffenson_iterate};

const gsl_root_fdfsolver_type  * gsl_root_fdfsolver_steffenson = &steffenson_type;

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
//
// C++ Interface: %{MODULE}
//
// Description:
//
//
// Author: %{AUTHOR} <%{EMAIL}>, (C) %{YEAR}
//
// Copyright: See COPYING file that comes with this distribution
//
//

#ifndef TNet_Vector_h
#define TNet_Vector_h

#include <cstddef>
#include <cstdlib>
#include <stdexcept>
#include <iostream>

#ifdef HAVE_ATLAS
extern "C"{
  #include <cblas.h>
  #include <clapack.h>
}
#endif

#include "Common.h"
#include "MathAux.h"
#include "Types.h"
#include "Error.h"

namespace TNet
{
  template<typename _ElemT> class Vector;
  template<typename _ElemT> class SubVector;
  template<typename _ElemT> class Matrix;
  template<typename _ElemT> class SpMatrix;

  // we need to declare the friend functions here
  template<typename _ElemT>
    std::ostream & operator << (std::ostream & rOut, const Vector<_ElemT> & rV);

  template<typename _ElemT>
    std::istream & operator >> (std::istream & rIn, Vector<_ElemT> & rV);

  template<typename _ElemT>
    _ElemT
    BlasDot(const Vector<_ElemT>& rA, const Vector<_ElemT>& rB);

  /** **************************************************************************
   ** **************************************************************************
   *  @brief Provides a matrix abstraction class
   *
   *  This class provides a way to work with matrices in TNet.
   *  It encapsulates basic operations and memory optimizations.
   *
   */
  template<typename _ElemT>
    class Vector
    {
    public:

    /// defines a type of this
    typedef Vector<_ElemT> ThisType;

    Vector(): mpData(NULL)
#ifdef STK_MEMALIGN_MANUAL
    ,mpFreeData(NULL)
#endif
    , mDim(0)
      {}

      /**
       * @brief Copy constructor
       * @param rV
       */
      Vector(const Vector<_ElemT>& rV)
  	  { mpData=NULL; Init(rV.Dim()); Copy(rV); }

      /* Type conversion constructor. */
      template<typename _ElemU>
      explicit Vector(const Vector<_ElemU>& rV)
  	  { mpData=NULL; Init(rV.Dim()); Copy(rV); }

      Vector(const _ElemT* ppData, const size_t s)
      { mpData=NULL; Init(s); Copy(ppData); }

      explicit Vector(const size_t s, bool clear=true)
      { mpData=NULL; Init(s,clear); }

      ~Vector()
      { Destroy(); }

       Vector<_ElemT> &operator = (const Vector <_ElemT> &other)
       { Init(other.Dim()); Copy(other); return *this; } // Needed for inclusion in std::vector

      Vector<_ElemT>&
      Init(size_t length, bool clear=true);

      /**
       * @brief Dealocates the window from memory and resets the dimensions to (0)
       */
      void
      Destroy();

      /**
       * @brief Returns @c true if vector is initialized
       */
      bool
      IsInitialized() const
      { return mpData != NULL; }

      /**
       * @brief Sets all elements to 0
       */
      void
      Zero();

      void
      Set(_ElemT f);

      inline size_t
      Dim() const
      { return mDim; }

      /**
       * @brief Returns size of matrix in memory (in bytes)
       */
      inline size_t
      MSize() const
      {
        return (mDim + (((16 / sizeof(_ElemT)) - mDim%(16 / sizeof(_ElemT)))
                          % (16 / sizeof(_ElemT)))) * sizeof(_ElemT);
      }

      /**
       *  @brief Gives access to the vector memory area
       *  @return pointer to the first field
       */
      inline _ElemT*
      pData()
      { return mpData; }

      /**
       *  @brief Gives access to the vector memory area
       *  @return pointer to the first field (const version)
       */
      inline const _ElemT*
      pData() const
      { return mpData; }

      /**
       *  @brief Gives access to a specified vector element (const).
       */
      inline _ElemT
      operator [] (size_t i) const
      {
#ifdef PARANOID
		assert(i<mDim);
#endif
		return *(mpData + i);
	  }

      /**
       *  @brief Gives access to a specified vector element (non-const).
       */
      inline _ElemT &
      operator [] (size_t i)
      {
#ifdef PARANOID
		assert(i<mDim);
#endif
		return *(mpData + i);
	  }

      /**
       *  @brief Gives access to a specified vector element (const).
       */
      inline _ElemT
		operator () (size_t i) const
      {
#ifdef PARANOID
		assert(i<mDim);
#endif
		return *(mpData + i);
	  }

      /**
       *  @brief Gives access to a specified vector element (non-const).
       */
      inline _ElemT &
		operator () (size_t i)
      {
#ifdef PARANOID
		assert(i<mDim);
#endif
		return *(mpData + i);
	  }

      /**
       * @brief Returns a matrix sub-range
       * @param o Origin
       * @param l Length
       * See @c SubVector class for details
       */
      SubVector<_ElemT>
      Range(const size_t o, const size_t l)
      { return SubVector<_ElemT>(*this, o, l); }

      /**
       * @brief Returns a matrix sub-range
       * @param o Origin
       * @param l Length
       * See @c SubVector class for details
       */
      const SubVector<_ElemT>
      Range(const size_t o, const size_t l) const
      { return SubVector<_ElemT>(*this, o, l); }

      //########################################################################
      //########################################################################

      /// Copy data from another vector
      Vector<_ElemT>&
      Copy(const Vector<_ElemT>& rV);

      /// Copy data from another vector of a different type.
      template<typename _ElemU> Vector<_ElemT>&
      Copy(const Vector<_ElemU>& rV);

      /// Load data into the vector
      Vector<_ElemT>&
      Copy(const _ElemT* ppData);

      Vector<_ElemT>&
      CopyVectorizedMatrixRows(const Matrix<_ElemT> &rM);

      Vector<_ElemT>&
      RemoveElement(size_t i);

      Vector<_ElemT>&
      ApplyLog();

      Vector<_ElemT>&
      ApplyLog(const Vector<_ElemT>& rV);//ApplyLog to rV and put the result in (*this)

      Vector<_ElemT>&
      ApplyExp();

      Vector<_ElemT>&
      ApplySoftMax();

      Vector<_ElemT>&
      Invert();

      Vector<_ElemT>&
      DotMul(const Vector<_ElemT>& rV); // Multiplies each element (*this)(i) by rV(i).

      Vector<_ElemT>&
      BlasAxpy(const _ElemT alpha, const Vector<_ElemT>& rV);

      Vector<_ElemT>&
      BlasGemv(const _ElemT alpha, const Matrix<_ElemT>& rM, const MatrixTrasposeType trans, const Vector<_ElemT>& rV, const _ElemT beta = 0.0);

      //########################################################################
      //########################################################################

      Vector<_ElemT>&
      Add(const Vector<_ElemT>& rV)
      { return BlasAxpy(1.0, rV); }

      Vector<_ElemT>&
      Subtract(const Vector<_ElemT>& rV)
      { return BlasAxpy(-1.0, rV); }

      Vector<_ElemT>&
      AddScaled(_ElemT alpha, const Vector<_ElemT>& rV)
      { return BlasAxpy(alpha, rV); }

      Vector<_ElemT>&
      Add(_ElemT c);

      Vector<_ElemT>&
      MultiplyElements(const Vector<_ElemT>& rV);

      // @brief elementwise : rV.*rR+beta*this --> this
      Vector<_ElemT>&
      MultiplyElements(_ElemT alpha, const Vector<_ElemT>& rV, const Vector<_ElemT>& rR,_ElemT beta);

      Vector<_ElemT>&
      DivideElements(const Vector<_ElemT>& rV);

      /// @brief elementwise : rV./rR+beta*this --> this
      Vector<_ElemT>&
      DivideElements(_ElemT alpha, const Vector<_ElemT>& rV, const Vector<_ElemT>& rR,_ElemT beta);

      Vector<_ElemT>&
      Subtract(_ElemT c);

      Vector<_ElemT>&
      Scale(_ElemT c);

      //########################################################################
      //########################################################################

      /// Performs a row stack of the matrix rMa
      Vector<_ElemT>&
      MatrixRowStack(const Matrix<_ElemT>& rMa);

      // Extracts a row of the matrix rMa.  .. could also do this with vector.Copy(rMa[row]).
      Vector<_ElemT>&
      Row(const Matrix<_ElemT>& rMa, size_t row);

      // Extracts a column of the matrix rMa.
      Vector<_ElemT>&
      Col(const Matrix<_ElemT>& rMa, size_t col);

      // Takes all elements to a power.
      Vector<_ElemT>&
      Power(_ElemT power);

      _ElemT
      Max() const;

      _ElemT
      Min() const;

      /// Returns sum of the elements
      _ElemT
      Sum() const;

      /// Returns sum of the elements
      Vector<_ElemT>&
      AddRowSum(const Matrix<_ElemT>& rM);

      /// Returns sum of the elements
      Vector<_ElemT>&
      AddColSum(const Matrix<_ElemT>& rM);

      /// Returns log(sum(exp())) without exp overflow
      _ElemT
      LogSumExp() const;

      //########################################################################
      //########################################################################

      friend std::ostream &
      operator << <> (std::ostream& rOut, const Vector<_ElemT>& rV);

      friend _ElemT
      BlasDot<>(const Vector<_ElemT>& rA, const Vector<_ElemT>& rB);

      /**
       * Computes v1^T * M * v2.
       * Not as efficient as it could be where v1==v2 (but no suitable blas
       * routines available).
       */
      _ElemT
      InnerProduct(const Vector<_ElemT> &v1, const Matrix<_ElemT> &M, const Vector<_ElemT> &v2) const;

    //##########################################################################
    //##########################################################################
    //protected:
    public:
      /// data memory area
      _ElemT*   mpData;
#ifdef STK_MEMALIGN_MANUAL
      /// data to be freed (in case of manual memalignment use, see common.h)
      _ElemT*   mpFreeData;
#endif
      size_t  mDim;      ///< Number of elements
    }; // class Vector

  /**
   * @brief Represents a non-allocating general vector which can be defined
   * as a sub-vector of higher-level vector
   */
  template<typename _ElemT>
    class SubVector : public Vector<_ElemT>
    {
    protected:
      /// Constructor
      SubVector(const Vector<_ElemT>& rT,
                const size_t  origin,
                const size_t  length)
      {
        assert(origin+length <= rT.mDim);
        Vector<_ElemT>::mpData = rT.mpData+origin;
        Vector<_ElemT>::mDim   = length;
      }
      //only Vector class can call this protected constructor
      friend class Vector<_ElemT>;

    public:
      /// Constructor
      SubVector(Vector<_ElemT>& rT,
                const size_t  origin,
                const size_t  length)
      {
        assert(origin+length <= rT.mDim);
        Vector<_ElemT>::mpData = rT.mpData+origin;
        Vector<_ElemT>::mDim   = length;
      }

      /**
       * @brief Constructs a vector representation out of a standard array
       *
       * @param pData pointer to data array to associate with this vector
       * @param length length of this vector
       */
      inline
      SubVector(_ElemT *ppData,
                size_t length)
      {
        Vector<_ElemT>::mpData = ppData;
        Vector<_ElemT>::mDim   = length;
      }

      /**
       * @brief Destructor
       */
      ~SubVector()
      {
        Vector<_ElemT>::mpData = NULL;
      }
    };

    // Useful shortcuts
    typedef Vector<BaseFloat> BfVector;
    typedef SubVector<BaseFloat> BfSubVector;

    //Adding two vectors of different types
    template <typename _ElemT, typename _ElemU>
    void Add(Vector<_ElemT>& rDst, const Vector<_ElemU>& rSrc)
    {
      assert(rDst.Dim() == rSrc.Dim());
      const _ElemU* p_src = rSrc.pData();
      _ElemT* p_dst = rDst.pData();

      for(size_t i=0; i<rSrc.Dim(); i++) {
        *p_dst++ += (_ElemT)*p_src++;
      }
    }

    //Scales adding two vectors of different types
    template <typename _ElemT, typename _ElemU>
    void AddScaled(Vector<_ElemT>& rDst, const Vector<_ElemU>& rSrc, _ElemT scale)
    {
      assert(rDst.Dim() == rSrc.Dim());

      Vector<_ElemT> tmp(rSrc);
      rDst.BlasAxpy(scale, tmp);

/*
      const _ElemU* p_src = rSrc.pData();
      _ElemT* p_dst = rDst.pData();

      for(size_t i=0; i<rDst.Dim(); i++) {
        *p_dst++ += *p_src++ * scale;
      }
*/
    }

} // namespace TNet

//*****************************************************************************
//*****************************************************************************
// we need to include the implementation
#include "Vector.tcc"

/******************************************************************************
 ******************************************************************************
 * The following section contains specialized template definitions
 * whose implementation is in Vector.cc
 */

#endif // #ifndef TNet_Vector_h

max_stars_repo_name: troylee/nnet-asr

max_issues_repo_name: troylee/nnet-asr

max_forks_repo_name: troylee/nnet-asr


text:
#define PY_SSIZE_T_CLEAN
#include <Python.h>
#include <gsl_vector.h>

#include "svm_main.h"
#include "test_main.h"

PyObject* __opt(PyObject* self, PyObject* args) {
    PyObject* elements;
    int k_type;    // 1 if RBF kernel is used 0 for my kernel
	double C;	// hyperparameter for optimization
	if (!PyArg_ParseTuple(args, "Oid", &elements, &k_type, &C))
		return NULL;

	return __get_w_b(elements, k_type, C);
}

PyObject* __test(PyObject* self, PyObject* args) {
	PyObject* training;
	PyObject* testing;
	PyObject* alphas;
	double b;
	int rbf;
	if (!PyArg_ParseTuple(args, "OOOdi", &training, &testing, &alphas, &b, &rbf))
		return NULL;
	return __test_t(training, testing, alphas, b, rbf);
}

static PyMethodDef methods[] = {
	{"opt", __opt, METH_VARARGS, ""},
	{"test", __test, METH_VARARGS, ""},
	{NULL, NULL, 0, NULL}
};

static struct PyModuleDef c_opt = {
	PyModuleDef_HEAD_INIT,
	"C OPT",
	"C interop to optimize a custom implementation of an SVM",
	-1,
	methods
};

PyMODINIT_FUNC PyInit_c_opt() {
	return PyModule_Create(&c_opt);
}


text:
#include <stdlib.h>
#include <cblas.h>

#include "tasks.h"

void syrk_task_par_reconfigure(int nth)
{
    // empty
}

void syrk_task_par_finalize(void)
{
    // empty
}

void syrk_task_par(void *ptr, int nth, int me)
{
    struct syrk_task_arg *arg = (struct syrk_task_arg*) ptr;

    int n       = arg->n;
    int k       = arg->k;
    double *A21 = arg->A21;
    double *A22 = arg->A22;
    int ldA     = arg->ldA;

    // Balance the load by flops.
    int part[nth + 1];
    const int total_work = n * (n + 1) / 2;
    const int ideal_part_work = total_work / nth;
    part[0] = 0;
    part[nth] = n;
    for (int k = 1; k < nth; ++k) {
        part[k] = part[k - 1];
        int work = 0;
        while (work < ideal_part_work && part[k] < n) {
            work += n - part[k];
            part[k] += 1;
        }
    }

    const int my_first_col = part[me];
    const int my_num_cols  = part[me + 1] - part[me];
    const int i1 = my_first_col;
    const int i2 = i1 + my_num_cols;
    const int m2 = n - i2;

    // TODO Insert picture (see krnl_syrk in task_chol_par).

    if (my_num_cols > 0) {
        cblas_dsyrk(CblasColMajor, CblasLower, CblasNoTrans,
                    my_num_cols, k,
                    -1.0, A21 + i1,                           ldA,
                     1.0, A22 + i1 + my_first_col * ldA,      ldA);
    }

    if (m2 > 0) {
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
                    m2, my_num_cols, k,
                    -1.0, A21 + i2,                      ldA,
                          A21 + my_first_col,            ldA,
                     1.0, A22 + i2 + my_first_col * ldA, ldA);
    }
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_system_schema_h_
#define SQ_INCLUDE_GUARD_system_schema_h_

#include "core/Primitive.h"
#include "core/typeutil.h"

#include <cstddef>
#include <gsl/gsl>
#include <optional>
#include <string_view>

namespace sq::system {

class SchemaImpl;

/**
 * Represents the schema for a primitive type.
 */
class PrimitiveTypeSchema {
public:
  constexpr explicit PrimitiveTypeSchema(std::string_view name,
                                         std::string_view doc)
      : name_{name}, doc_{doc} {}

  SQ_ND std::string_view name() const;
  SQ_ND std::string_view doc() const;

private:
  std::string_view name_;
  std::string_view doc_;
};

/**
 * Represents the schema for a parameter of field of a system object.
 */
class ParamSchema {
public:
  constexpr ParamSchema(std::string_view name, std::string_view doc,
                        std::size_t index, std::size_t type_index,
                        bool required, std::string_view default_value_str,
                        std::string_view default_value_doc)
      : name_{name}, doc_{doc}, index_{index}, type_index_{type_index},
        required_{required}, default_value_str_{default_value_str},
        default_value_doc_{default_value_doc} {}

  SQ_ND std::string_view name() const;
  SQ_ND std::string_view doc() const;
  SQ_ND std::size_t index() const;
  SQ_ND const PrimitiveTypeSchema &type() const;
  SQ_ND bool required() const;
  SQ_ND std::optional<Primitive> default_value() const;
  SQ_ND std::string_view default_value_doc() const;

private:
  std::string_view name_;
  std::string_view doc_;
  std::size_t index_;
  std::size_t type_index_;
  bool required_;
  std::string_view default_value_str_;
  std::string_view default_value_doc_;
};

class TypeSchema;

/**
 * Represents the schema for a field of a system object.
 */
class FieldSchema {
public:
  constexpr FieldSchema(std::string_view name, std::string_view doc,
                        std::size_t params_begin_index,
                        std::size_t params_end_index,
                        std::size_t return_type_index, bool return_list,
                        bool null)
      : name_{name}, doc_{doc}, params_begin_index_{params_begin_index},
        params_end_index_{params_end_index},
        return_type_index_{return_type_index},
        return_list_{return_list}, null_{null} {}

  SQ_ND std::string_view name() const;
  SQ_ND std::string_view doc() const;
  SQ_ND gsl::span<const ParamSchema> params() const;
  SQ_ND const TypeSchema &return_type() const;
  SQ_ND bool return_list() const;
  SQ_ND bool null() const;

private:
  std::string_view name_;
  std::string_view doc_;
  std::size_t params_begin_index_;
  std::size_t params_end_index_;
  std::size_t return_type_index_;
  bool return_list_;
  bool null_;
};

/**
 * Represents the schema for a system object.
 */
class TypeSchema {
public:
  constexpr TypeSchema(std::string_view name, std::string_view doc,
                       std::size_t fields_begin_index,
                       std::size_t fields_end_index)
      : name_{name}, doc_{doc}, fields_begin_index_{fields_begin_index},
        fields_end_index_{fields_end_index} {}

  SQ_ND std::string_view name() const;
  SQ_ND std::string_view doc() const;
  SQ_ND gsl::span<const FieldSchema> fields() const;

private:
  std::string_view name_;
  std::string_view doc_;
  std::size_t fields_begin_index_;
  std::size_t fields_end_index_;
};

/**
 * Represents the whole SQ schema.
 */
struct Schema {
public:
  SQ_ND gsl::span<const TypeSchema> types() const;
  SQ_ND gsl::span<const PrimitiveTypeSchema> primitive_types() const;
  SQ_ND const TypeSchema &root_type() const;
};

/**
 * Get the whole SQ schema.
 */
const Schema &schema();

} // namespace sq::system

#endif // SQ_INCLUDE_GUARD_system_schema_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>

/* slow/dumb rmedian which constructs actual window for each sample, sorts
 * it and finds median */
static int
slow_rmedian(const gsl_filter_end_t endtype, const gsl_vector * x, gsl_vector * y, const int K)
{
  const int n = (int) x->size;
  const int H = K / 2;
  double *window = malloc(K * sizeof(double));
  int i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(endtype, x, i, H, H, window);
      double yi;
      int j;

      /* fill first half of window with previous filter output values */
      for (j = i - H; j < i; ++j)
        {
          if (j < 0)
            {
              if (endtype == GSL_FILTER_END_PADVALUE)
                window[j - i + H] = gsl_vector_get(x, 0);
              else if (endtype == GSL_FILTER_END_PADZERO)
                window[j - i + H] = 0.0;
            }
          else
            {
              window[j - i + H] = gsl_vector_get(y, j);
            }
        }

      yi = gsl_stats_median(window, 1, wsize);
      gsl_vector_set(y, i, yi);
    }

  free(window);

  return GSL_SUCCESS;
}

/* test square wave input (root signal) */
static void
test_rmedian_root(const gsl_filter_end_t etype, const size_t n, const size_t k)
{
  const double tol = 1.0e-12;
  gsl_filter_rmedian_workspace *w = gsl_filter_rmedian_alloc(k);
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  char buf[2048];
  size_t i;

  /* test a root sequence (square input): x = [zero one zero] */
  gsl_vector_set_all(x, 0.0);

  for (i = n / 3; i <= n / 2; ++i)
    gsl_vector_set(x, i, 1.0);

  /* compute y = rmedian(x) and test y = x */
  gsl_filter_rmedian(etype, x, y, w);

  sprintf(buf, "n=%zu k=%zu RMF square wave root sequence", n, k);
  compare_vectors(tol, y, x, buf);

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_filter_rmedian_free(w);
}

/* test random input and in-place */
static void
test_rmedian_random(const gsl_filter_end_t etype, const size_t n, const int K, gsl_rng * r)
{
  const double tol = 1.0e-12;
  gsl_filter_rmedian_workspace *w = gsl_filter_rmedian_alloc(K);
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *z = gsl_vector_alloc(n);
  char buf[2048];

  /* test filter with random input against slow algorithm */
  random_vector(x, r);

  /* y = rmedian(x) */
  gsl_filter_rmedian(etype, x, y, w);

  /* y = rmedian(x) with slow algorithm */
  slow_rmedian(etype, x, z, K);

  /* test y = z */
  sprintf(buf, "n=%zu K=%d RMF symmetric random slow test", n, K);
  compare_vectors(tol, y, z, buf);

  /* test in-place filter */

  /* z = rmedian(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_filter_rmedian(etype, z, z, w);

  sprintf(buf, "n=%zu K=%d RMF symmetric random in-place", n, K);
  compare_vectors(tol, z, y, buf);

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
  gsl_filter_rmedian_free(w);
}

void
test_rmedian(gsl_rng * rng_p)
{
  /* test root sequences */

  test_rmedian_root(GSL_FILTER_END_PADZERO, 1000, 3);
  test_rmedian_root(GSL_FILTER_END_PADZERO, 2000, 101);

  test_rmedian_root(GSL_FILTER_END_PADVALUE, 1000, 3);
  test_rmedian_root(GSL_FILTER_END_PADVALUE, 2000, 101);

  /* test random input */

  test_rmedian_random(GSL_FILTER_END_PADZERO, 10, 1, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADZERO, 100, 3, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADZERO, 1000, 3, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADZERO, 100, 1001, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADZERO, 5, 7, rng_p);

  test_rmedian_random(GSL_FILTER_END_PADVALUE, 10, 1, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADVALUE, 100, 3, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADVALUE, 1000, 3, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADVALUE, 100, 1001, rng_p);
  test_rmedian_random(GSL_FILTER_END_PADVALUE, 5, 7, rng_p);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#ifndef QDM_GSL_H
#define QDM_GSL_H 1

#include <stdio.h>

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

gsl_vector *
qdm_vector_seq(double from, double to, double by);

void
qdm_vector_set_seq(gsl_vector *v, double from, double to);

size_t
qdm_vector_search(const gsl_vector *v, double needle);

void
qdm_vector_csv_fwrite(FILE *f, const gsl_vector *v);

void
qdm_matrix_csv_fwrite(FILE *f, const gsl_matrix *m);

int
qdm_matrix_tmm(gsl_matrix *m, gsl_matrix *result);

int
qdm_matrix_det_tmm(gsl_matrix *m, double *det);

gsl_vector *
qdm_vector_sorted(const gsl_vector *v);

gsl_vector *
qdm_vector_quantile(gsl_vector *data, gsl_vector *probs);

double
qdm_vector_rss(const gsl_vector *y, const gsl_vector *fx);

double
qdm_vector_sum(gsl_vector *v);

size_t
qdm_vector_greater_than(const gsl_vector *v, double value);

void
qdm_matrix_select_upper_triangle(gsl_matrix *m);

int
qdm_matrix_to_csc_matrix(csc **result, gsl_matrix *m);

gsl_vector *
qdm_vector_copy(const gsl_vector *src);

gsl_matrix *
qdm_matrix_copy(const gsl_matrix *src);

int
qdm_vector_hd5_read(
    hid_t id,
    const char *name,
    gsl_vector **v
);

int
qdm_vector_hd5_write(
    hid_t id,
    const char *name,
    const gsl_vector *v
);

int
qdm_matrix_hd5_read(
    hid_t id,
    const char *name,
    gsl_matrix **m
);

int
qdm_matrix_hd5_write(
    hid_t id,
    const char *name,
    const gsl_matrix *m
);

gsl_vector *
qdm_matrix_filter(
    const gsl_matrix *m,
    size_t needle_column,
    double needle_value,
    size_t select_column
);

#endif /* QDM_GSL_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
/*
 *  specialfunctionsmodule.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef SPECIALFUNCTIONSMODULE_H
#define SPECIALFUNCTIONSMODULE_H
/*
    SLI Module implementing functions from the GNU Science Library.
    The GSL is available from sources.redhat.com/gsl.
*/

/*
    NOTE: Special functions are available only if the GSL is installed.
          If no GSL is available, calling special functions will result
          in a SLI error message.  HEP 2002-09-19.
*/

// Generated includes:
#include "config.h"

// Includes from sli:
#include "slifunction.h"
#include "slimodule.h"

#ifdef HAVE_GSL
// External include:
#include <gsl/gsl_integration.h>
#endif

// NOTE: all gsl headers are included in specialfunctionsmodule.cc

class SpecialFunctionsModule : public SLIModule
{

  // Part 1: Methods pertaining to the module ----------------------

public:
  SpecialFunctionsModule( void ){};
  // ~SpecialFunctionsModule(void);

  // The Module is registered by a call to this Function:
  void init( SLIInterpreter* );

  // This function will return the name of our module:
  const std::string name( void ) const;

  // Part 2: Classes for the implemented functions -----------------

public:
  /**
   * Classes which implement the GSL Funktions.
   * These must be public, since we want to export
   * objects of these.
   */
  class GammaIncFunction : public SLIFunction
  {
  public:
    GammaIncFunction()
    {
    }
    void execute( SLIInterpreter* ) const;
  };
  class LambertW0Function : public SLIFunction
  {
  public:
    LambertW0Function()
    {
    }
    void execute( SLIInterpreter* ) const;
  };
  class LambertWm1Function : public SLIFunction
  {
  public:
    LambertWm1Function()
    {
    }
    void execute( SLIInterpreter* ) const;
  };

  class ErfFunction : public SLIFunction
  {
  public:
    ErfFunction()
    {
    }
    void execute( SLIInterpreter* ) const;
  };

  class ErfcFunction : public SLIFunction
  {
  public:
    ErfcFunction()
    {
    }
    void execute( SLIInterpreter* ) const;
  };

  class GaussDiskConvFunction : public SLIFunction
  {
  public:
    void execute( SLIInterpreter* ) const;

    // need constructor and destructor to set up integration workspace
    GaussDiskConvFunction( void );
    ~GaussDiskConvFunction( void );

  private:
    // quadrature parameters, see GSL Reference
    static const int MAX_QUAD_SIZE;
    static const double QUAD_ERR_LIM;
    static const double QUAD_ERR_SCALE;

// integration workspace
#ifdef HAVE_GSL
    gsl_integration_workspace* w_;

    /**
     * Integrand function.
     * @note This function must be static with C linkage so that it can
     *       be passed to the GSL. Alternatively, one could define it
     *       outside the class.
     */
    static double f_( double, void* );
    static gsl_function F_; // GSL wrapper struct for it
#endif
  };

  // Part 3: One instatiation of each new function class -----------

public:
  const GammaIncFunction gammaincfunction;
  const LambertW0Function lambertw0function;
  const LambertWm1Function lambertwm1function;
  const ErfFunction erffunction;
  const ErfcFunction erfcfunction;
  const GaussDiskConvFunction gaussdiskconvfunction;

  // Part 3b: Internal variables
private:
};

// Part 4: Documentation for all functions -------------------------

/* BeginDocumentation

Name: Gammainc - incomplete gamma function

Synopsis: x a Gammainc -> result

Description: Computes the incomplete Gamma function
             int(t^(a-1)*exp(-t), t=0..x) / Gamma(a)

Parameters:  x      (double): upper limit of integration
             a      (double): order of Gamma function

Examples: 2.2 1.5 Gammainc -> 0.778615

Author: H E Plesser

FirstVersion: 2001-07-26

Remarks: This is the incomplete Gamma function P(a,x) defined as no. 6.5.1
         in Abramowitz&Stegun.  Requires the GSL.

References: http://sources.redhat.com/gsl/ref
*/

/* BeginDocumentation

Name: Erf - error function

Synopsis: x Erf -> result

Description: Computes the error function
             erf(x) = 2/sqrt(pi) int_0^x dt exp(-t^2)

Parameters:  x (double): error function argument

Examples: 0.5 erf -> 0.5205

Author: H E Plesser

FirstVersion: 2001-07-30

Remarks: Requires the GSL.

References: http://sources.redhat.com/gsl/ref

SeeAlso: Erfc
*/

/* BeginDocumentation

Name: Erfc - complementary error function

Synopsis: x Erfc -> result

Description: Computes the error function
             erfc(x) = 1 - erf(x) = 2/sqrt(pi) int_x^inf dt exp(-t^2)

Parameters:  x (double): error function argument

Examples: 0.5 erfc -> 0.4795

Author: H E Plesser

FirstVersion: 2001-07-30

Remarks: Requires the GSL.

References: http://sources.redhat.com/gsl/ref

SeeAlso: Erf
*/

/* BeginDocumentation

Name:GaussDiskConv - Convolution of a Gaussian with an excentric disk

Synopsis:R r0 GaussDiskConv -> result

Description:Computes the convolution of an excentric normalized Gaussian
with a disk

       C[R, r0] = IInt[ disk(rvec; R) * Gauss(rvec - r_0vec) d^2rvec ]
                = 2 Int[ r Exp[-r0^2-r^2] I_0[2 r r_0] dr, r=0..R]

Parameters:R   radius of the disk, centered at origin
r0  distance of Gaussian center from origin

Examples:SLI ] 3.2 2.3 GaussDiskConv =
0.873191

Author:H E Plesser

FirstVersion: 2002-07-12

Remarks:This integral is needed to compute the response of a DOG model to
 excentric light spots, see [1].  For technicalities, see [2].  Requires GSL.

References: [1] G. T. Einevoll and P. Heggelund, Vis Neurosci 17:871-885 (2000).
 [2] Hans E. Plesser, Convolution of an Excentric Gaussian with a Disk,
     Technical Report, arken.nlh.no/~itfhep, 2002

*/

#endif

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


text:
/* interpolation/gsl_spline2d.h
 *
 * Copyright 2012 David Zaslavsky
 *
 * 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.
 */

#ifndef __GSL_SPLINE2D_H__
#define __GSL_SPLINE2D_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/*
 * A 2D interpolation object which stores the arrays defining the function.
 * In all other respects, this is just like a gsl_interp2d object.
 */
typedef struct
{
  gsl_interp2d interp_object; /* low-level interpolation object */
  double * xarr;              /* x data array */
  double * yarr;              /* y data array */
  double * zarr;              /* z data array */
} gsl_spline2d;

GSL_FUN gsl_spline2d * gsl_spline2d_alloc(const gsl_interp2d_type * T, size_t xsize, size_t ysize);

GSL_FUN int gsl_spline2d_init(gsl_spline2d * interp, const double xa[],
                      const double ya[], const double za[],
                      size_t xsize, size_t ysize);

GSL_FUN void gsl_spline2d_free(gsl_spline2d * interp);

GSL_FUN double gsl_spline2d_eval(const gsl_spline2d * interp, const double x,
                         const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

GSL_FUN int gsl_spline2d_eval_e(const gsl_spline2d * interp, const double x,
                        const double y, gsl_interp_accel* xa, gsl_interp_accel* ya,
                        double * z);

GSL_FUN double gsl_spline2d_eval_deriv_x(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

GSL_FUN int gsl_spline2d_eval_deriv_x_e(const gsl_spline2d * interp, const double x,
                                const double y, gsl_interp_accel* xa,
                                gsl_interp_accel* ya, double * z);

GSL_FUN double gsl_spline2d_eval_deriv_y(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya);

GSL_FUN int gsl_spline2d_eval_deriv_y_e(const gsl_spline2d * interp, const double x,
                                const double y, gsl_interp_accel* xa,
                                gsl_interp_accel* ya, double * z);

GSL_FUN double gsl_spline2d_eval_deriv_xx(const gsl_spline2d * interp, const double x,
                                  const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

GSL_FUN int gsl_spline2d_eval_deriv_xx_e(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya, double * z);

GSL_FUN double gsl_spline2d_eval_deriv_yy(const gsl_spline2d * interp, const double x,
                                  const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

GSL_FUN int gsl_spline2d_eval_deriv_yy_e(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya, double * z);

GSL_FUN double gsl_spline2d_eval_deriv_xy(const gsl_spline2d * interp, const double x,
                                  const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

GSL_FUN int gsl_spline2d_eval_deriv_xy_e(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya, double * z);

GSL_FUN size_t gsl_spline2d_min_size(const gsl_spline2d * interp);

GSL_FUN const char * gsl_spline2d_name(const gsl_spline2d * interp);

GSL_FUN int gsl_spline2d_set(const gsl_spline2d * interp, double zarr[],
                     const size_t i, const size_t j, const double z);
GSL_FUN double gsl_spline2d_get(const gsl_spline2d * interp, const double zarr[],
                        const size_t i, const size_t j);

__END_DECLS

#endif /* __GSL_SPLINE2D_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#pragma once

#include <gsl/gsl>

#include <vector>
#include <memory>
#include <mutex>

struct AVFrame;

namespace video_streamer
{

class StreamSession;
class FrameEncoder;

class GroupStreamSession
{
public:
    explicit GroupStreamSession(std::shared_ptr<FrameEncoder> frameEncoder);

    void addStream(std::unique_ptr<StreamSession> streamSession);
    void pushFrame(gsl::not_null<const AVFrame*> frame);

private:
    std::shared_ptr<FrameEncoder> m_frameEncoder;
    std::vector<std::unique_ptr<StreamSession>> m_streamSession;
    mutable std::mutex m_mtx;
};

} // video_streamer

max_stars_repo_name: nicledomaS/VideoStream

max_issues_repo_name: nicledomaS/VideoStream

max_forks_repo_name: nicledomaS/VideoStream


text:
// Copyright (c) 2021 SmartPolarBear
//
// 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.

//
// Created by cleve on 9/1/2021.
//
#pragma once

#include <base/iterable_stack.h>

#include <interpreter/vm/opcode.h>
#include <interpreter/vm/chunk.h>
#include <interpreter/vm/value.h>
#include <interpreter/vm/exceptions.h>

#include <interpreter/vm/heap.h>
#include <interpreter/vm/heap_allocator.h>
#include <interpreter/vm/garbage_collector.h>

#include <interpreter/vm/string_object.h>
#include <interpreter/vm/closure_object.h>
#include <interpreter/vm/instance_object.h>

#include <memory>
#include <map>
#include <ranges>

#include <gsl/gsl>
#include "class_object.h"

namespace clox::interpreting::vm
{

enum class virtual_machine_status
{
	OK,
	RUNTIME_ERROR,
};

template<typename F, typename R, class ...Args>
concept BinaryOperator=
requires(F&& f, Args&& ... args) {
	{ std::invoke(std::forward<F>(f), std::forward<Args>(args)...) }->std::same_as<R>;
};

class virtual_machine final
{
public:
	friend class garbage_collector;

	static inline constexpr size_t CALL_STACK_RESERVED_SIZE = 64;
	static inline constexpr size_t STACK_RESERVED_SIZE = 16384;

	using value_list_type = std::vector<value>;
	using global_table_type = std::unordered_map<std::string, value>;
	using function_table_type = std::unordered_map<full_opcode_type, value>;
	using ip_type = chunk::iterator_type;

	using index_type = gsl::index;

	class call_frame final
	{
	public:

	public:

		explicit call_frame(closure_object_raw_pointer closure, ip_type ip, size_t offset)
				: closure_(closure), ip_(ip), stack_offset_(offset)
		{
		}

		[[nodiscard]] size_t stack_offset() const
		{
			return stack_offset_;
		}

		[[nodiscard]] function_object_raw_pointer function() const
		{
			return closure_->function();
		}

		[[nodiscard]] closure_object_raw_pointer closure() const
		{
			return closure_;
		}

		[[nodiscard]] ip_type& ip()
		{
			return ip_;
		}

	private:
		closure_object_raw_pointer closure_{};
		ip_type ip_{};
		size_t stack_offset_{};
	};

	using call_frame_list_type = std::vector<call_frame>;

public:
	virtual_machine() = delete;

	~virtual_machine();

	explicit virtual_machine(helper::console& cons,
			std::shared_ptr<object_heap> heap);

	virtual_machine_status run(clox::interpreting::vm::closure_object* closure);

private:
	virtual_machine_status run();

	// {return status, exit}
	std::tuple<std::optional<virtual_machine_status>, bool> run_code(chunk::code_type instruction, call_frame& frame);

	template<class ...TArgs>
	void runtime_error(std::string_view fmt, const TArgs& ...args)
	{

		cons_->error() << std::format("[Line {}] in file {}:",
				top_call_frame().function()->body()->line_of(top_call_frame().ip()),
				top_call_frame().function()->body()->filename()) << std::endl;

		cons_->error() << std::format(fmt, args...) << "\n";

		cons_->error() << "Call stack:" << std::endl;
		for (auto& frm: call_frames_ | std::ranges::views::reverse)
		{
			auto func = frm.function();
			auto line = func->body()->line_of(frm.ip() - 1);
			cons_->error() << std::format("[Line {}] in ", line);
			if (func->name().empty())
			{
				cons_->error() << "script" << std::endl;
			}
			else
			{
				cons_->error() << func->name() << std::endl;
			}
		}

		reset_stack();
	}

	template<typename TOp>
	requires BinaryOperator<TOp, floating_value_type, floating_value_type, floating_value_type>
	inline void binary_op(TOp op)
	{
		try
		{
			auto l = peek(0), r = peek(1);

			auto right = get_number_promoted(l);

			auto left = get_number_promoted(r);

			auto ret = op(left, right);

			if (std::holds_alternative<floating_value_type>(l) ||
				std::holds_alternative<floating_value_type>(r))
			{
				pop_two_and_push(ret);
			}
			else if (std::holds_alternative<integer_value_type>(l) ||
					 std::holds_alternative<integer_value_type>(r))
			{
				pop_two_and_push(static_cast<integer_value_type>(ret));
			}
			else if (std::holds_alternative<boolean_value_type>(l) ||
					 std::holds_alternative<boolean_value_type>(r))
			{
				pop_two_and_push(static_cast<scanning::boolean_literal_type>(ret));
			}
			else
			{
				pop_two_and_push(
						static_cast<integer_value_type>(ret)); // cannot combine for the sake of the rules of type promoting
			}
		}
		catch (const std::exception& e)
		{
			this->runtime_error("Invalid operands for binary operator: {}", e.what());
			throw e;
		}
	}

	template<typename TOp>
	requires BinaryOperator<TOp, object_raw_pointer, string_object_raw_pointer, string_object_raw_pointer>
	inline void binary_op(TOp op)
	{
		try
		{
			auto right = get_string(peek(0));
			auto left = get_string(peek(1));

			pop_two_and_push(op(left, right));
		}
		catch (const std::exception& e)
		{
			this->runtime_error("Invalid operands for binary operator: {}", e.what());
			throw e;
		}
	}

	// for comparing, whose op returns bool value
	template<typename TOp>
	requires BinaryOperator<TOp, bool, floating_value_type, floating_value_type>
	inline void binary_op(TOp op)
	{
		try
		{
			auto right = get_number_promoted(peek(0));
			auto left = get_number_promoted(peek(1));

			pop_two_and_push(op(left, right));
		}
		catch (const std::exception& e)
		{
			this->runtime_error("Invalid operands for binary operator: {}", e.what());
			throw e;
		}
	}

	template<typename TOp>
	requires BinaryOperator<TOp, bool, string_object_raw_pointer, string_object_raw_pointer>
	inline void binary_op(TOp op)
	{
		try
		{
			auto right = get_string(peek(0));
			auto left = get_string(peek(1));

			pop_two_and_push(op(left, right));
		}
		catch (const std::exception& e)
		{
			this->runtime_error("Invalid operands for binary operator: {}", e.what());
			throw e;
		}
	}

	void call_value(const value& val, size_t arg_count);

	void call(closure_object_raw_pointer closure, size_t arg_count);

	bool is_false(const value& val);

	[[maybe_unused]] bool is_true(const value& val)
	{
		return !is_false(val);
	}

	// stack modification
	void reset_stack();

	template<object_pointer T>
	T peek_object(size_t offset = 0)
	{
		return std::visit([](auto&& val) -> T
		{
			using TV = std::decay_t<decltype(val)>;
			if constexpr (std::is_same_v<TV, object_raw_pointer>)
			{
				if (auto func = dynamic_cast<T>(val);func)
				{
					return func;
				}
				else
				{
					throw invalid_value{ val };
				}
			}
			else
			{
				throw invalid_value{ val };
			}
		}, peek(offset));
	}

	value peek(size_t offset = 0);

	value pop();

	void push(const value& val);

	inline void pop_two_and_push(const value& val);
	//

	// instruction reading
	value next_constant();

	std::string next_variable_name();

	chunk::code_type next_code();

	value& slot_at(const call_frame& frame, size_t slot);

	upvalue_object_raw_pointer capture_upvalue(value* val, index_type stack_index);

	void close_upvalues(index_type last);

	// call frame

	call_frame& top_call_frame()
	{
		return *call_frames_.rbegin();
	}

	void push_call_frame(closure_object_raw_pointer closure, chunk::iterator_type ip, size_t stack_offset)
	{
		call_frames_.emplace_back(closure, ip, stack_offset);
	}

	void pop_call_frame()
	{
		call_frames_.pop_back();
	}

	//

	// method
	bool bind_method(class_object_raw_pointer class_obj, resolving::function_id_type id);

	bool bind_method(instance_object_raw_pointer class_obj, resolving::function_id_type id);

	std::shared_ptr<object_heap> heap_{};

	value_list_type stack_{};

	global_table_type globals_{};

	function_table_type functions_{};

	call_frame_list_type call_frames_{};

	std::map<index_type, upvalue_object_raw_pointer> open_upvalues_{}; // it should be ordered

	mutable helper::console* cons_{ nullptr };
};

}

max_stars_repo_name: SmartPolarBear/clox

max_issues_repo_name: SmartPolarBear/clox

max_forks_repo_name: SmartPolarBear/clox


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

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

#include <cblas.h>
#include "morn_tensor.h"

struct TensorReshapePara
{
    MLayer *prev;

    int channel;
    int height;
    int width;

    int res_valid;
};

void *mTensorReshapePara(MList *ini,char *name)
{
    struct TensorReshapePara *para = (struct TensorReshapePara *)mMalloc(sizeof(struct TensorReshapePara));

    char *value = mINIRead(ini,name,"prev");
    para->prev = mNetworkLayer(ini,value);
    mException((para->prev == NULL),EXIT,"invalid prev");

    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    value = mINIRead(ini,name,"height");
    if(value != NULL) para->height= atoi(value);else para->height= DFLT;

    value = mINIRead(ini,name,"width");
    if(value != NULL) para->width = atoi(value);else para->width = DFLT;

    value = mINIRead(ini,name,"channel");
    if(value != NULL) para->channel = atoi(value);else para->channel = DFLT;

    return para;
}

void TensorReshapeSet(MLayer *layer)
{
    if(layer->state != DFLT) return;
    struct TensorReshapePara *para = (struct TensorReshapePara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out= layer->tns;

    int size = in->height*in->width*in->channel;

    if(para->height <=0) para->height = size/(para->width *para->channel);
    if(para->width  <=0) para->width  = size/(para->height*para->channel);
    if(para->channel<=0) para->channel= size/(para->height*para->width  );
    mException((para->height<0)||(para->width<0)||(para->channel<0)||(para->height*para->width*para->channel!=size),EXIT,"invalid input shape, with shape input is %d,%d,%d",in->channel,in->height,in->width);

    mTensorRedefine(out,in->batch,para->channel,para->height,para->width,in->data);
    if(morn_network_flag==MORN_TRAIN)
    {
        if(INVALID_TENSOR(res)) mTensorRedefine(res,in->batch,in->channel,in->height,in->width,in->data);
        else                    mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);
    }
}

void mTensorReshapeForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Reshape",mLayerType(layer)),EXIT,"invalid layer type");

    TensorReshapeSet(layer);

    layer->state = MORN_FORWARD;
}

void mTensorReshapeBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Reshape",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorReshapePara *para = (struct TensorReshapePara *)(layer->para);
    if(para->res_valid==0) return;

    MTensor *res= para->prev->res;
    MTensor *out= layer->res;
    int size = res->height*res->width*res->channel;
    for(int b=0;b<res->batch;b++)if(res->data[b]!=out->data[b]) memcpy(res->data[b],out->data[b],size*sizeof(float));

    para->prev->state = MORN_BACKWARD;
}

max_stars_repo_name: ishine/Morn

max_issues_repo_name: ishine/Morn

max_forks_repo_name: ishine/Morn


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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "ccl.h"

/*
 * convert_concentration_single finds the concentration c'
 * for a mass definition with overdensity Delta' given
 * the concentration c for a mass definition with overdensity
 * Delta assuming an NFW density profile.
 *
 * To do so, it solves the following equation numerically:
 *    c^3 * f(c) = Delta' f(c')
 * where f(x) = x^3 * (x + 1) / ((1+x) * log(x + 1) - x).
 *
 * The equation is solved using a Newton-Raphson approach.
 * The functions nfw_fx, nfw_f, nfw_df and nfw_fdf implement
 * the function whose zero we try to find and its derivative.
 *
 * ccl_convert_concentration does the same thing for an
 * array of input concentrations.
 */

static double nfw_fx(double x)
{
  if(x>0.01) {
    double xp1=1+x;
    double lxp1=log(xp1);
    double den=1./(xp1*lxp1-x);
    return xp1*x*x*x*den;
  }
  else {
    return 2*x;
  }
}

static double nfw_f(double x,void *params)
{
  double offset = *((double *)params);
  return nfw_fx(x)-offset;
}

static double nfw_df(double x,void *params)
{
  if(x>0.01) {
    double xp1=1+x;
    double lxp1=log(xp1);
    double den=1./(xp1*lxp1-x);
    return x*x*(3*xp1*xp1*lxp1-x*(4*x+3))*den*den;
  }
  else {
    return 2;
  }
}

static void nfw_fdf(double x,void *params,
		    double *y, double *dy)
{
  double offset = *((double *)params);
  if(x>0.01) {
    double xp1=1+x;
    double lxp1=log(xp1);
    double den=1./(xp1*lxp1-x);
    *y=xp1*x*x*x*den-offset;
    *dy=x*x*(3*xp1*xp1*lxp1-x*(4*x+3))*den*den;
  }
  else {
    *y=2*x-offset;
    *dy=2;
  }
}

static int convert_concentration_single(double d_factor, double c_old,
					double *c_new, double c_start)
{
  double c0, offset = d_factor * nfw_fx(c_old);
  int status, iter=0, max_iter=100;
  gsl_function_fdf FDF;
  gsl_root_fdfsolver *s=NULL;
  FDF.f = &nfw_f;
  FDF.df = &nfw_df;
  FDF.fdf = &nfw_fdf;
  FDF.params = &offset;

  s=gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
  if (s==NULL)
    return CCL_ERROR_MEMORY;

  gsl_root_fdfsolver_set (s, &FDF, c_start);
  *c_new = c_start;
  do
    {
      iter++;
      c0 = *c_new;
      status = gsl_root_fdfsolver_iterate (s);
      *c_new = gsl_root_fdfsolver_root (s);
      status = gsl_root_test_delta (*c_new, c0, 0, 1e-4);

    }
  while (status == GSL_CONTINUE && iter < max_iter);
  gsl_root_fdfsolver_free (s);

  return status;
}

void ccl_convert_concentration(ccl_cosmology *cosmo,
			       double delta_old, int nc, double c_old[],
			       double delta_new, double c_new[],int *status)
{
  if(nc<=0)
    return;

  int ii,st=0;
  double d_factor = delta_old/delta_new;
  for(ii=0;ii<nc;ii++) {
    double c_new_h;
    st=convert_concentration_single(d_factor, c_old[ii], &(c_new[ii]), c_old[ii]);
    if(st!=GSL_SUCCESS) {
      *status=CCL_ERROR_ROOT;
      ccl_cosmology_set_status_message(cosmo,
				       "ccl_mass_conversion.c: NR solver failed to find a root\n");
      return;
    }
  }
}

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
#pragma once

#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <DirectXMath.h>
#include <gsl\gsl>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"

namespace Library
{
	class SkyboxMaterial;

	class Skybox final : public DrawableGameComponent
	{
		RTTI_DECLARATIONS(Skybox, DrawableGameComponent)

	public:
		Skybox(Game& game, const std::shared_ptr<Camera>& camera, const std::wstring& cubeMapFileName, float scale);
		Skybox(const Skybox&) = delete;
		Skybox(Skybox&&) = default;
		Skybox& operator=(const Skybox&) = delete;
		Skybox& operator=(Skybox&&) = default;
		~Skybox() = default;

		virtual void Initialize() override;
		virtual void Draw(const GameTime& gameTime) override;

	private:
		std::wstring mCubeMapFileName;
		DirectX::XMFLOAT4X4 mWorldMatrix{ MatrixHelper::Identity };
		float mScale;
		std::shared_ptr<SkyboxMaterial> mMaterial;
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		bool mUpdateMaterial{ true };
	};
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gbpInterpolate.h>

double interpolate(interp_info *interp, double x) {
    return (gsl_interp_eval(interp->interp, interp->x, interp->y, x, interp->accel));
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#pragma once

#include <gsl/span>

#include "halley/text/halleystring.h"
#include "halley/resources/resource.h"
#include "halley/maths/range.h"
#include "halley/maths/vector4.h"
#include "halley/data_structures/maybe.h"
#include <map>
#include <vector>

#include "halley/utils/type_traits.h"

#if defined(DEV_BUILD)
#define STORE_CONFIG_NODE_PARENTING
#endif

namespace Halley {
	class Serializer;
	class Deserializer;
	class ConfigNode;

	template<typename T>
	struct HasToConfigNode
	{
	private:
		typedef std::true_type yes;
		typedef std::false_type no;
		template<typename U> static auto test(int) -> decltype(std::declval<U>().toConfigNode(), yes());
		template<typename> static no test(...);

	public:
		static constexpr bool value = std::is_same<decltype(test<T>(0)),yes>::value;
	};

	template<typename T>
	struct HasConfigNodeConstructor
	{
	private:
		typedef std::true_type yes;
		typedef std::false_type no;
		template<typename U> static auto test(int) -> decltype(U(std::declval<ConfigNode>()), yes());
		template<typename> static no test(...);

	public:
		static constexpr bool value = std::is_same<decltype(test<T>(0)),yes>::value;
	};

    enum class ConfigNodeType
    {
        Undefined,
        String,
        Sequence,
        Map,
        Int,
        Float,
        Int2,
        Float2,
        Int3,
        Float3,
        Int4,
        Float4,
        Bytes,
        DeltaSequence, // For delta coding
        DeltaMap,      // For delta coding
        Noop,          // For delta coding
        Idx,           // For delta coding
        Del            // For delta coding
    };

    template <>
    struct EnumNames< ConfigNodeType >
    {
        constexpr std::array< const char*, 17 > operator()() const
        {
            return { { "undefined",
                "string",
                "sequence",
                "map",
                "int",
                "float",
                "int2",
                "float2",
                "int3",
                "float3",
                "int4",
                "float4",
                "bytes",
                "deltaSequence",
                "deltaMap",
                "noop",
                "idx"
                "del" } };
        }
    };

    class ConfigFile;

    class ConfigNode
    {
        friend class ConfigFile;

    public:
        using MapType = std::map< String, ConfigNode >;
        using SequenceType = std::vector< ConfigNode >;

        struct NoopType
        {
        };
        struct DelType
        {
        };
        struct IdxType
        {
            int start;
            int len;
            IdxType() = default;
            IdxType( int start, int len ) :
                start( start ), len( len ) { }
        };

        ConfigNode();
        explicit ConfigNode( const ConfigNode& other );
        ConfigNode( ConfigNode&& other ) noexcept;
        ConfigNode( MapType entryMap );
        ConfigNode( SequenceType entryList );
        explicit ConfigNode( String value );
        explicit ConfigNode( const char* value );
        explicit ConfigNode( bool value );
        explicit ConfigNode( int value );
        explicit ConfigNode( float value );
        explicit ConfigNode( Angle1f value );
        explicit ConfigNode( Vector2i value );
        explicit ConfigNode( Vector2f value );
        explicit ConfigNode( Vector3i value );
        explicit ConfigNode( Vector3f value );
        explicit ConfigNode( Vector4i value );
        explicit ConfigNode( Vector4f value );
        explicit ConfigNode( Bytes value );
        explicit ConfigNode( NoopType value );
        explicit ConfigNode( DelType value );
        explicit ConfigNode( IdxType value );

		template <typename T>
		explicit ConfigNode(const std::vector<T>& sequence)
		{
			SequenceType seq;
			seq.reserve(sequence.size());
			for (auto& e: sequence) {
				if constexpr (HasToConfigNode<T>::value) {
					seq.push_back(e.toConfigNode());
				} else {
					seq.push_back(ConfigNode(e));
				}
			}
			*this = seq;
		}

        ~ConfigNode();

        ConfigNode& operator=( const ConfigNode& other );
        ConfigNode& operator=( ConfigNode&& other ) noexcept;
        ConfigNode& operator=( bool value );
        ConfigNode& operator=( int value );
        ConfigNode& operator=( float value );
        ConfigNode& operator=( Angle1f value );
        ConfigNode& operator=( Vector2i value );
        ConfigNode& operator=( Vector2f value );
        ConfigNode& operator=( Vector3i value );
        ConfigNode& operator=( Vector3f value );
        ConfigNode& operator=( Vector4i value );
        ConfigNode& operator=( Vector4f value );

        ConfigNode& operator=( MapType entryMap );
        ConfigNode& operator=( SequenceType entryList );
        ConfigNode& operator=( String value );
        ConfigNode& operator=( Bytes value );
        ConfigNode& operator=( gsl::span< const gsl::byte > bytes );

        ConfigNode& operator=( const char* value );

        ConfigNode& operator=( NoopType value );
        ConfigNode& operator=( DelType value );
        ConfigNode& operator=( IdxType value );

		template <typename T>
		ConfigNode& operator=(const std::vector<T>& sequence)
		{
			SequenceType seq;
			seq.reserve(sequence.size());
			for (auto& e: sequence) {
				if constexpr (HasToConfigNode<T>::value) {
					seq.push_back(e.toConfigNode());
				} else {
					seq.push_back(ConfigNode(e));
				}
			}
			return *this = seq;
		}

        bool operator==( const ConfigNode& other ) const;
        bool operator!=( const ConfigNode& other ) const;

        ConfigNodeType getType() const;

        void serialize( Serializer& s ) const;
        void deserialize( Deserializer& s );

        int asInt() const;
        float asFloat() const;
        bool asBool() const;
        Angle1f asAngle1f() const;
        Vector2i asVector2i() const;
        Vector2f asVector2f() const;
        Vector3i asVector3i() const;
        Vector3f asVector3f() const;
        Vector4i asVector4i() const;
        Vector4f asVector4f() const;
        Range< float > asFloatRange() const;
        String asString() const;
        const Bytes& asBytes() const;

        int asInt( int defaultValue ) const;
        float asFloat( float defaultValue ) const;
        bool asBool( bool defaultValue ) const;
        String asString( const String& defaultValue ) const;
        Angle1f asAngle1f( Angle1f defaultValue ) const;
        Vector2i asVector2i( Vector2i defaultValue ) const;
        Vector2f asVector2f( Vector2f defaultValue ) const;
        Vector3i asVector3i( Vector3i defaultValue ) const;
        Vector3f asVector3f( Vector3f defaultValue ) const;
        Vector4i asVector4i( Vector4i defaultValue ) const;
        Vector4f asVector4f( Vector4f defaultValue ) const;

		template <typename T>
		std::vector<T> asVector() const
		{
			if (type == ConfigNodeType::Sequence) {
				std::vector<T> result;
				result.reserve(asSequence().size());
				for (const auto& e : asSequence()) {
					if constexpr (HasConfigNodeConstructor<T>::value) {
						result.emplace_back(T(e));
					} else {
						result.emplace_back(e.convertTo(Tag<T>()));
					}
				}
				return result;
			} else {
				throw Exception("Can't convert " + getNodeDebugId() + " from " + toString(getType()) + " to std::vector<T>.", HalleyExceptions::Resources);
			}
		}

		template <typename T>
		std::vector<T> asVector(const std::vector<T>& defaultValue) const
		{
			if (type == ConfigNodeType::Sequence) {
				return asVector<T>();
			} else {
				return defaultValue;
			}
		}

        const SequenceType& asSequence() const;
        const MapType& asMap() const;
        SequenceType& asSequence();
        MapType& asMap();

        void ensureType( ConfigNodeType type );

        bool hasKey( const String& key ) const;
        void removeKey( const String& key );

        ConfigNode& operator[]( const String& key );
        ConfigNode& operator[]( size_t idx );

        const ConfigNode& operator[]( const String& key ) const;
        const ConfigNode& operator[]( size_t idx ) const;

        SequenceType::iterator begin();
        SequenceType::iterator end();

        SequenceType::const_iterator begin() const;
        SequenceType::const_iterator end() const;

        void reset();
        void setOriginalPosition( int line, int column );
        void setParent( const ConfigNode* parent, int idx );
        void propagateParentingInformation( const ConfigFile* parentFile );

        inline void assertValid() const
        {
            Expects( intData != 0xCDCDCDCD );
            Expects( intData != 0xDDDDDDDD );
        }

        struct BreadCrumb
        {
            const BreadCrumb* prev = nullptr;
            String key;
            OptionalLite< int > idx;
            int depth = 0;

            BreadCrumb() = default;
            BreadCrumb( const BreadCrumb& prev, String key ) :
                prev( &prev ), key( std::move( key ) ), depth( prev.depth + 1 ) { }
            BreadCrumb( const BreadCrumb& prev, int index ) :
                prev( &prev ), idx( index ), depth( prev.depth + 1 ) { }

            bool hasKeyAt( const String& key, int depth ) const;
            bool hasIndexAt( int idx, int depth ) const;
        };

        class IDeltaCodeHints
        {
        public:
            virtual ~IDeltaCodeHints() = default;

            virtual std::optional< size_t > getSequenceMatch( const SequenceType& seq, const ConfigNode& newValue, size_t curIdx, const BreadCrumb& breadCrumb ) const = 0;
            virtual bool doesSequenceOrderMatter( const BreadCrumb& breadCrumb ) const { return true; }
            virtual bool canDeleteKey( const String& key, const BreadCrumb& breadCrumb ) const { return true; }
            virtual bool canDeleteAnyKey() const { return true; }
            virtual bool shouldBypass( const BreadCrumb& breadCrumb ) const { return false; }
            virtual bool areNullAndEmptyEquivalent( const BreadCrumb& breadCrumb ) const { return false; }
        };

        static ConfigNode createDelta( const ConfigNode& from, const ConfigNode& to, const IDeltaCodeHints* hints = nullptr );
        static ConfigNode applyDelta( const ConfigNode& from, const ConfigNode& delta );
        void applyDelta( const ConfigNode& delta );
        void decayDeltaArtifacts();

    private:
        template < typename T >
        class Tag
        {
        };

        union {
            String* strData;
            MapType* mapData;
            SequenceType* sequenceData;
            Bytes* bytesData;
            void* rawPtrData;
            int intData;
            float floatData;
            Angle1f ang1fData;
            Vector2i vec2iData;
            Vector2f vec2fData;
            Vector3i vec3iData;
            Vector3f vec3fData;
            Vector4i vec4iData;
            Vector4f vec4fData;
        };
        ConfigNodeType type = ConfigNodeType::Undefined;
        int auxData = 0; // Used by delta coding

#if defined(STORE_CONFIG_NODE_PARENTING)
        struct ParentingInfo
        {
            int line = 0;
            int column = 0;
            int idx = 0;
            const ConfigNode* node = nullptr;
            const ConfigFile* file = nullptr;
        };
        std::unique_ptr< ParentingInfo > parent;
#endif
        static ConfigNode undefinedConfigNode;
        static String undefinedConfigNodeName;

        template < typename T >
        void deserializeContents( Deserializer& s )
        {
            T v;
            s >> v;
            *this = std::move( v );
        }

        String getNodeDebugId() const;
        String backTrackFullNodeName() const;

        int convertTo( Tag< int > tag ) const;
        float convertTo( Tag< float > tag ) const;
        bool convertTo( Tag< bool > tag ) const;
        Angle1f convertTo( Tag< Angle1f > tag ) const;
        Vector2i convertTo( Tag< Vector2i > tag ) const;
        Vector2f convertTo( Tag< Vector2f > tag ) const;
        Vector3i convertTo( Tag< Vector3i > tag ) const;
        Vector3f convertTo( Tag< Vector3f > tag ) const;
        Vector4i convertTo( Tag< Vector4i > tag ) const;
        Vector4f convertTo( Tag< Vector4f > tag ) const;
        Range< float > convertTo( Tag< Range< float > > tag ) const;
        String convertTo( Tag< String > tag ) const;
        const Bytes& convertTo( Tag< Bytes& > tag ) const;

        bool isNullOrEmpty() const;

        static ConfigNode doCreateDelta( const ConfigNode& from, const ConfigNode& to, const BreadCrumb& breadCrumb, const IDeltaCodeHints* hints );
        static ConfigNode createMapDelta( const ConfigNode& from, const ConfigNode& to, const BreadCrumb& breadCrumb, const IDeltaCodeHints* hints );
        static ConfigNode createSequenceDelta( const ConfigNode& from, const ConfigNode& to, const BreadCrumb& breadCrumb, const IDeltaCodeHints* hints );
        void applyMapDelta( const ConfigNode& delta );
        void applySequenceDelta( const ConfigNode& delta );

        bool isEquivalent( const ConfigNode& other ) const;
        bool isEquivalentStrictOrder( const ConfigNode& other ) const;
    };
} // namespace Halley

max_stars_repo_name: chidddy/halley

max_issues_repo_name: chidddy/halley

max_forks_repo_name: chidddy/halley


text:
/* sum/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sum.h>

#include <gsl/gsl_ieee_utils.h>

#define N 50

void check_trunc (double * t, double expected, const char * desc);
void check_full (double * t, double expected, const char * desc);

int
main (void)
{
  gsl_ieee_env_setup ();

  {
    double t[N];
    int n;

    const double zeta_2 = M_PI * M_PI / 6.0;

    /* terms for zeta(2) */

    for (n = 0; n < N; n++)
      {
	double np1 = n + 1.0;
	t[n] = 1.0 / (np1 * np1);
      }

    check_trunc (t, zeta_2, "zeta(2)");
    check_full (t, zeta_2, "zeta(2)");
  }

  {
    double t[N];
    double x, y;
    int n;

    /* terms for exp(10.0) */
    x = 10.0;
    y = exp(x);

    t[0] = 1.0;
    for (n = 1; n < N; n++)
      {
	t[n] = t[n - 1] * (x / n);
      }

    check_trunc (t, y, "exp(10)");
    check_full (t, y, "exp(10)");
  }

  {
    double t[N];
    double x, y;
    int n;

    /* terms for exp(-10.0) */
    x = -10.0;
    y = exp(x);

    t[0] = 1.0;
    for (n = 1; n < N; n++)
      {
	t[n] = t[n - 1] * (x / n);
      }

    check_trunc (t, y, "exp(-10)");
    check_full (t, y, "exp(-10)");
  }

  {
    double t[N];
    double x, y;
    int n;

    /* terms for -log(1-x) */
    x = 0.5;
    y = -log(1-x);
    t[0] = x;
    for (n = 1; n < N; n++)
      {
	t[n] = t[n - 1] * (x * n) / (n + 1.0);
      }

    check_trunc (t, y, "-log(1/2)");
    check_full (t, y, "-log(1/2)");
  }

  {
    double t[N];
    double x, y;
    int n;

    /* terms for -log(1-x) */
    x = -1.0;
    y = -log(1-x);
    t[0] = x;
    for (n = 1; n < N; n++)
      {
	t[n] = t[n - 1] * (x * n) / (n + 1.0);
      }

    check_trunc (t, y, "-log(2)");
    check_full (t, y, "-log(2)");
  }

  {
    double t[N];
    int n;

    double result = 0.192594048773;

    /* terms for an alternating asymptotic series */

    t[0] = 3.0 / (M_PI * M_PI);

    for (n = 1; n < N; n++)
      {
	t[n] = -t[n - 1] * (4.0 * (n + 1.0) - 1.0) / (M_PI * M_PI);
      }

    check_trunc (t, result, "asymptotic series");
    check_full (t, result, "asymptotic series");
  }

  {
    double t[N];
    int n;

    /* Euler's gamma from GNU Calc (precision = 32) */

    double result = 0.5772156649015328606065120900824;

    /* terms for Euler's gamma */

    t[0] = 1.0;

    for (n = 1; n < N; n++)
      {
	t[n] = 1/(n+1.0) + log(n/(n+1.0));
      }

    check_trunc (t, result, "Euler's constant");
    check_full (t, result, "Euler's constant");
  }

  {
    double t[N];
    int n;

    /* eta(1/2) = sum_{k=1}^{\infty} (-1)^(k+1) / sqrt(k)

       From Levin, Intern. J. Computer Math. B3:371--388, 1973.

       I=(1-sqrt(2))zeta(1/2)
        =(2/sqrt(pi))*integ(1/(exp(x^2)+1),x,0,inf) */

    double result = 0.6048986434216305;  /* approx */

    /* terms for eta(1/2) */

    for (n = 0; n < N; n++)
      {
	t[n] = (n%2 ? -1 : 1) * 1.0 /sqrt(n + 1.0);
      }

    check_trunc (t, result, "eta(1/2)");
    check_full (t, result, "eta(1/2)");
  }

  exit (gsl_test_summary ());
}

void
check_trunc (double * t, double expected, const char * desc)
{
  double sum_accel, prec;

  gsl_sum_levin_utrunc_workspace * w = gsl_sum_levin_utrunc_alloc (N);

  gsl_sum_levin_utrunc_accel (t, N, w, &sum_accel, &prec);
  gsl_test_rel (sum_accel, expected, 1e-8, "trunc result, %s", desc);

  /* No need to check precision for truncated result since this is not
     a meaningful number */

  gsl_sum_levin_utrunc_free (w);
}

void
check_full (double * t, double expected, const char * desc)
{
  double sum_accel, err_est, sd_actual, sd_est;

  gsl_sum_levin_u_workspace * w = gsl_sum_levin_u_alloc (N);

  gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err_est);
  gsl_test_rel (sum_accel, expected, 1e-8, "full result, %s", desc);

  sd_est = -log10 (err_est/fabs(sum_accel));
  sd_actual = -log10 (DBL_EPSILON + fabs ((sum_accel - expected)/expected));

  /* Allow one digit of slop */

  gsl_test (sd_est > sd_actual + 1.0, "full significant digits, %s (%g vs %g)", desc, sd_est, sd_actual);

  gsl_sum_levin_u_free (w);
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_cdf.h>

unsigned int k;
unsigned int n1;
unsigned int n2;
unsigned int t;

static void usage(int argc, char *argv[]) {
  printf (
  "Usage: %s k n1 n2 t\n"
  "\n"
  "If a population contains n1 elements of \"type 1\" and n2 elements\n"
  "of \"type 2\" then the hypergeometric distribution gives the probability\n"
  "of obtaining k elements of \"type 1\" in t samples from the population\n"
  "without replacement.\n"
  "This program computes the Pvalue of such a distribution\n", argv[0]);
  exit (-1);
}

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

  if (argc != 5)
    usage(argc, argv);

  k = atoi(argv[1]);
  n1 = atoi(argv[2]);
  n2 = atoi(argv[3]);
  t = atoi(argv[4]);

  double Q = gsl_cdf_hypergeometric_Q (k, n1, n2, t);
  double P = gsl_cdf_hypergeometric_P (k, n1, n2, t);

  printf("%g\n", Q);
  if (P < Q)
    fprintf(stderr, "Warning: P < Q, %g < %g\n", P, Q);

}

max_stars_repo_name: scovit/nust-helico

max_issues_repo_name: scovit/nust-helico

max_forks_repo_name: scovit/nust-helico