text:
/* population genetics tools- mainly for working with site frequency spectrum
/
/
/ Andrew Kern
*/

#include "stdio.h"
#include <math.h>
#include "snpFile.h"
#include "popGenTools.h"
#include "numerical.h"
#include "assert.h"
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf.h>

/* function for sfs generation -- stationary solution for snp frequency density*/
double my_f(double q, void * p){
  struct my_f_params * params = (struct my_f_params *) p;
  double beta = params->beta;
  int i = params->i;
  int n = params->n;

  return ((1.0 - exp(-2.0 * beta * (1.0 - q))) / (1.0 - exp(-2.0 * beta))) * (2.0 / (q * (1.0 - q))) * gsl_ran_binomial_pdf(i, q, n);
}

/* integral of f over allele freq */
double my_F(double beta, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_f_params fParams;
  gsl_function f;
  //gsl_integration_workspace * w = gsl_integration_workspace_alloc (100000);
  double result,error;
  size_t size= 1000;

  if (beta == 0){
    //gsl_integration_workspace_free(w);
    return 1.0 / params->i;
  }
  else{
    fParams.beta = beta;
    fParams.i = params->i;
    fParams.n = params->n;
    f.function = &(my_f);
    f.params = &fParams;
    //adaptive integration-- slow and steady
//	gsl_integration_qags(&f,0.0,1.0,1e-5,1e-5, 100000,w, &result,&error);
	//non-adaptive version is rougher and faster
	gsl_integration_qng(&f, 0.0, 1.0, 1e-2,1e-2,&result,&error,&size);
//    gsl_integration_workspace_free(w);
    return result;
  }
}

/* integral of f over allele freq using numerical recipes integration */
double my_F2(double beta, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_f_params fParams;
  double result;

  if (beta == 0){
    return 1.0 / params->i;
  }
  else{
    fParams.beta = beta;
    fParams.i = params->i;
    fParams.n = params->n;
    result = d_qromb2(&my_f, 0.0, 1.0, &fParams);
    return result;
  }
}

/* integral of f over allele freq-- using adaptive integration */
double my_FAdapt(double beta, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_f_params fParams;
  gsl_function f;
  double result,error;
  size_t size= 100;

  if (beta == 0){
    return 1.0 / params->i;
  }
  else{
    fParams.beta = beta;
    fParams.i = params->i;
    fParams.n = params->n;
    f.function = &(my_f);
    f.params = &fParams;
    //adaptive integration-- slow and steady
	gsl_integration_qags(&f,0.0,1.0,1e-3,1e-3,size,params->w, &result,&error);
    return result;
  }
}

/* this returns the prob of SNP freq i given n and beta */
double snpProb(double beta, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_F_params tempParams;
  long double tot;
  int i;
  gsl_function F;

  /* add up total prob space */
  tot = 0.0;
  tempParams.n = params->n;

  for(i = 1; i < params->n ; i++){
    tempParams.i = i;
    tempParams.w = params->w;
    F.function =  &my_F;
    F.params = &tempParams;
    tot += GSL_FN_EVAL(&F, beta);
  }

  /* reset parameters */
  F.function =  &my_FAdapt;
  F.params = params;

  return   GSL_FN_EVAL(&F, beta) / tot ;
}

/* trying to speedup the snpProb routine with a lookup of the denominators (in snpProb_params) */
double snpProbDenomLookup(double beta, void * p){
	struct my_snpProb_params * params = (struct my_snpProb_params *) p;
	struct my_F_params tempParams;
	gsl_function F;

/* set parameters */
	tempParams.n = params->n;
	tempParams.i = params->i;
	tempParams.w = params->w;
	F.function =  &my_FAdapt;
	F.params = &tempParams;

	return   GSL_FN_EVAL(&F, beta) / params->denom ;
}

/* sampleSize vector - returns a vector of 1s or 0s depending on sample sizes in the data */
gsl_vector *sampleSizeVector(struct snp data[], int snpNumber, int maxSampleSize){
  gsl_vector *bools;
  int i;

  //alloc and initialize bools vector
  bools = gsl_vector_alloc(maxSampleSize + 1);
  gsl_vector_set_zero(bools);

  //go through data, set each sampleSize value to 1
  for(i = 0; i < snpNumber; i++){
    gsl_vector_set(bools,data[i].n, 1);
  }
  return(bools);
}

/* sampleSize vector - returns a vector of 1s or 0s depending on sample sizes in the data */
void sampleSizeVectorFill(struct snp data[], int snpNumber, int maxSampleSize, gsl_vector *bools){
  int i;

  //initialize bools vector
  gsl_vector_set_zero(bools);

  //go through data, set each sampleSize value to 1
  for(i = 0; i < snpNumber; i++){
    gsl_vector_set(bools,data[i].n, 1);
  }
}
/* returns a vector of snpProbDenominators from i = 2 to n, 1 indexed */
gsl_vector *makeSnpProbDenomVector(double beta, int maxSampleSize, gsl_vector *sampleSizeVector){
	gsl_vector *probs;
	int i, j, test;
	double tot;
	struct my_F_params tempParams;
	gsl_function F;
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
	probs = gsl_vector_alloc(maxSampleSize + 1);

	tempParams.w = w;
/* go from n = 2 to maxSampleSize, calc total prob space for each, put in vector probs */
	for(j = 2; j < maxSampleSize + 1; j++){
	//check if sampleSize is included
		test = gsl_vector_get(sampleSizeVector,j);
		if (test){
			tot = 0.0;
			tempParams.n = j;
			for(i = 1; i < j ; i++){
				tempParams.i = i;
				F.function =  &my_FAdapt;
				F.params = &tempParams;
				tot += GSL_FN_EVAL(&F, beta);
			}
			gsl_vector_set(probs, j, tot);
		}
	}
	gsl_integration_workspace_free(w);
	return probs;
}

/* snpProbDenom - returns denom for snpProb */
double snpProbDenom(double beta, int sampleSize){
  int i;
  double tot;
  struct my_F_params tempParams;
  gsl_function F;

  tot = 0.0;
  tempParams.n = sampleSize;
  for(i = 1; i < sampleSize ; i++){
	tempParams.i = i;
	F.function =  &my_F;
	F.params = &tempParams;
	tot += GSL_FN_EVAL(&F, beta);
      }
  return(tot);
}

/* snpProbMatrix- returns a matrix of log snpProbs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1) */
gsl_matrix *snpProbMatrix(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools){
  gsl_vector *denoms;
  gsl_matrix *probs;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //get denoms
  denoms = makeSnpProbDenomVector(beta,maxSampleSize, sampleSizeVector);
  FParams.w = gsl_integration_workspace_alloc(1000);
  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//have freq j?
	if(gsl_matrix_get(sfsBools, i, j)){
	  //calc prob
	  FParams.i = j;
	  FParams.n = i;
	  FParams.denom = gsl_vector_get(denoms,i);
	  f.function = &snpProbDenomLookup;
	  f.params = &FParams;
	  gsl_matrix_set(probs, i, j, log(GSL_FN_EVAL(&f, beta)));
	}
      }
    }
  }
  gsl_vector_free(denoms);
  gsl_integration_workspace_free(FParams.w);
  return(probs);
}
/* snpProbMatrix- returns a matrix of  snpProbs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
variable sampleSizes (2 to max), columns represent freqs (1 to max -1) */
gsl_matrix *snpProbMatrixNotLog(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools){
	gsl_vector *denoms;
	gsl_matrix *probs;
	struct my_snpProb_params FParams;
	gsl_function f;
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (100000);
	int i, j;

	FParams.w = w;
//alloc probs
	probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
	gsl_matrix_set_zero(probs);
//get denoms
	denoms = makeSnpProbDenomVector(beta,maxSampleSize, sampleSizeVector);

//go through sampleSizes
	for(i = 2; i <= maxSampleSize; i++){
	//have sample size i?
		if (gsl_vector_get(sampleSizeVector, i)){
	//go through freqs
			for(j = 1; j < maxSampleSize; j++){
	//have freq j?
				if(gsl_matrix_get(sfsBools, i, j)){
	//calc prob
					FParams.i = j;
					FParams.n = i;
					FParams.denom = gsl_vector_get(denoms,i);
					f.function = &snpProbDenomLookup;
					f.params = &FParams;
					gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, beta));
				}
			}
		}
	}
	gsl_vector_free(denoms);
	return(probs);
}
/* snpProbMatrixNotLogFull- returns a matrix of  snpProbs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1). All freqs/samplesizes included */
gsl_matrix *snpProbMatrixNotLogFull(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector){
  gsl_vector *denoms;
  gsl_matrix *probs;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);

  //get denoms

  denoms = makeSnpProbDenomVector(beta,maxSampleSize, sampleSizeVector);

  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //check for sampleSize
    if (gsl_vector_get(sampleSizeVector,i)){
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//calc prob
	FParams.i = j;
	FParams.n = i;
	FParams.denom = gsl_vector_get(denoms,i);
	f.function = &snpProbDenomLookup;
	f.params = &FParams;
	gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, beta));
      }
    }
  }
  gsl_vector_free(denoms);
  return(probs);
}

/* snpProbVectorNotLog- returns a vector of log snpProbs at specified sampleSize columns represent freqs (1 to max -1) */
gsl_vector *snpProbVectorNotLog(double beta,  int sampleSize){
  double denom;
  gsl_vector *probs;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i;

  //alloc probs
  probs = gsl_vector_alloc(sampleSize);
  gsl_vector_set_zero(probs);
  //get denoms
  denom = snpProbDenom(beta, sampleSize);

  //go through freqs
  for(i = 1; i < sampleSize; i++){
    FParams.i = i;
    FParams.n = sampleSize;
    FParams.denom = denom;
    f.function = &snpProbDenomLookup;
    f.params = &FParams;
    gsl_vector_set(probs, i, GSL_FN_EVAL(&f, beta));
  }
  return(probs);
}

gsl_vector *snpProbVectorNotLogPre(double beta,  int sampleSize, gsl_vector *probs){
  double denom;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i;

   gsl_vector_set_zero(probs);
  //get denoms
  denom = snpProbDenom(beta, sampleSize);

  //go through freqs
  for(i = 1; i < sampleSize; i++){
    FParams.i = i;
    FParams.n = sampleSize;
    FParams.denom = denom;
    f.function = &snpProbDenomLookup;
    f.params = &FParams;
    gsl_vector_set(probs, i, GSL_FN_EVAL(&f, beta));
  }
  return(probs);
}

/* likelihood function for SFS given the data; p here contains the data and the weights */
double my_lik(double beta, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	struct my_F_params FParams;
	gsl_function f;
	double lik;
	int i;

	lik = 0;
	for(i = 0; i < params->snpNumber; i++){
		FParams.i = params->data[i].i;
		FParams.n = params->data[i].n;
		f.function = &snpProb;
		f.params = &FParams;
		lik += log(GSL_FN_EVAL(&f, beta));
		}
	return -lik;
}

/* likelihood function for SFS where each SNP has independent selection coeff */
double sfsLikBetaVector(gsl_vector *betas, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	struct my_F_params FParams;
	gsl_function f;
	double lik;
	int i;

	lik = 0;
	for(i = 0; i < params->snpNumber; i++){
		FParams.i = params->data[i].i;
		FParams.n = params->data[i].n;
		FParams.w = params->w;
		f.function = &snpProb;
		f.params = &FParams;
		//	printf("here lik:%f beta:%f i: %d n: %d\n",lik,gsl_vector_get(betas, i),FParams.i,FParams.n);

		lik += log(GSL_FN_EVAL(&f, gsl_vector_get(betas, i)));
		}
	return -lik;
}

/*likelihood from Bustamante et al. 2001 */
double bustaLik(double beta, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	struct my_F_params FParams;
	gsl_function f;
	double lik;
	int i;

	lik = 0;
	for(i = 0; i < params->snpNumber; i++){
		FParams.i = params->data[i].i;
		FParams.n = params->data[i].n;
		f.function = &snpProb;
		f.params = &FParams;
		lik += log(GSL_FN_EVAL(&f, beta));
		}
	return -lik;
}

/* weighted likelihood function for SFS given the data; here contains the data and the weights */
double weightedLik(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  struct my_F_params FParams;
  gsl_function f;
  double lik;
  int i;

  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    FParams.i = params->data[i].i;
    FParams.n = params->data[i].n;
    f.function = &snpProb;
    f.params = &FParams;
    lik += gsl_vector_get(params->weights, i) * log(GSL_FN_EVAL(&f, beta));
  }
  return -lik;
}

/* weighted likelihood function for SFS given the data; here contains the data and the weights and uses the snpProbLookup routine */
double weightedLikLook(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  double lik;
  int i;
  gsl_matrix *probs;

  //make prob matrix (note this is in logs)
  probs = snpProbMatrix(beta, params->maxSampleSize, params->sampleSizeVector, params->sfsBools);
  //now go and tally likelihood
  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)  * gsl_vector_get(params->weights,i);
  }
  gsl_matrix_free(probs);
  return -lik;
}

/* Folded SFS version of above; weighted likelihood function for SFS given the data; here contains the data and the weights and uses the snpProbLookup routine */
double weightedLikLookFolded(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  double lik;
  int i;
  gsl_matrix *probs;

  //make prob matrix (note this is in logs)
  probs = snpProbMatrix(beta, params->maxSampleSize, params->sampleSizeVector, params->sfsBools);
  //now go and tally likelihood
  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)  * gsl_vector_get(params->weights,i);
	lik += gsl_matrix_get(probs,params->data[i].n, params->data[i].n - params->data[i].i)  * gsl_vector_get(params->weights,i);
  }
  gsl_matrix_free(probs);
  return -lik;
}

/* PRF log Likelihood  */
double prfLik(double beta, void * p){
	 struct prfLik_params * params = (struct prfLik_params *) p;
	double lik;
	double term1, term2,theta;
  	int i;
	struct my_F_params tempParams;
	gsl_function F;
	double fTerm = 0.;

	theta = gsl_vector_get(params->betas,1);
	term1 = term2 = 0;
	for(i = 1; i < params->maxSampleSize;i++){
		tempParams.i = i;
		tempParams.n = params->maxSampleSize;
		F.function = &my_F;
		F.params = &tempParams;
		fTerm = GSL_FN_EVAL(&F,beta);
		term1 += gsl_matrix_get(params->sfsCounts,params->maxSampleSize,i) * log( theta* fTerm) ;
		term2 += fTerm;
	}
	return((term1 - (theta * term2)) * -1);
}

double weightedLikLookCI(double beta, void * p){
  struct my_likCI_params *params = (struct my_likCI_params *) p;
  struct my_lik_params likParams;
  gsl_function l;
  double result;

  likParams.data = params->data;
  likParams.snpNumber = params->snpNumber;
  likParams.sampleSizeVector = params->sampleSizeVector;
  likParams.maxSampleSize = params->maxSampleSize;
  likParams.weights = params->weights;
  likParams.sfsBools = params->sfsBools;
  l.function = &weightedLikLook;
  l.params = &likParams;
  result = GSL_FN_EVAL(&l, beta) - (params->lMax) - (params->logUnits) ;
  //  printf("%f\t%f\t%f\n",params->lMax,GSL_FN_EVAL(&l, beta),result);
  return result;
}

double weightedLikLookCIFolded(double beta, void * p){
  struct my_likCI_params *params = (struct my_likCI_params *) p;
  struct my_lik_params likParams;
  gsl_function l;
  double result;

  likParams.data = params->data;
  likParams.snpNumber = params->snpNumber;
  likParams.sampleSizeVector = params->sampleSizeVector;
  likParams.maxSampleSize = params->maxSampleSize;
  likParams.weights = params->weights;
  likParams.sfsBools = params->sfsBools;
  l.function = &weightedLikLookFolded;
  l.params = &likParams;
  result = GSL_FN_EVAL(&l, beta) - (params->lMax) - (params->logUnits) ;
  //  printf("%f\t%f\t%f\n",params->lMax,GSL_FN_EVAL(&l, beta),result);
  return result;
}

/* wrapper for my_lik for use in mnbrak */
double likWrap(double beta){
	gsl_function f;

	f.function = &my_lik;
	f.params = NULL;
	return GSL_FN_EVAL(&f, beta);
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrap(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLik;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrapLook(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLikLook;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

double  wlikWrapLookFolded(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLikLookFolded;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

double  wlikCIWrapLook(double beta, void * p){
	gsl_function f;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	f.function = &weightedLikLookCI;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/* wrapper for prfLik for use in mnbrak */
double  prfLikWrap(double beta, void * p){
	gsl_function f;
	struct prfLik_params * params = (struct prfLik_params *) p;

	f.function = &prfLik;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/*error checking function for likelihood function */
void printProfile(double min, double max, void * p){
  gsl_function f;
  struct my_lik_params *params =  (struct my_lik_params *) p;
  int i;

  f.function = &weightedLikLook;
  f.params = params;
  for(i = min; i < max; i++){
    printf("%d\t%f\n",i, GSL_FN_EVAL(&f, i));
  }

}

/* get ML estimate of beta using likelihood function; BRENT method for gsl */
double ml_est(double * lik_beta_hat){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc;
	gsl_function L;

	/* get minimum bracket */
	ax = -10.0;
	bx = 0.0;
	cx = 10.0;

	mnbrak(&ax, &bx, &cx, &fa, &fb, &fc, &my_lik);

	/* check if there is a minimum, if not return inf */
	if (fa == fb || fc == fb || ax > cx){
		*lik_beta_hat = fa;
		return 666.0;
		}
	else{

		/* do mimization */

		L.function = &my_lik;
		//should have param line here but no params....

		T = gsl_min_fminimizer_brent;
		s = gsl_min_fminimizer_alloc(T);
		gsl_min_fminimizer_set(s, &L, bx, ax, cx);

		do{
			iter++;
			status = gsl_min_fminimizer_iterate(s);
			 m = gsl_min_fminimizer_x_minimum (s);
			  a = gsl_min_fminimizer_x_lower (s);
			  b = gsl_min_fminimizer_x_upper (s);

			  status = gsl_min_test_interval (a, b, 0.001, 0.0);

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

	  return m;
	}
}

/* get weighted ML estimate of beta using likelihood function; BRENT method for gsl */
double weighted_ml_est(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, min;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -10.0;
	bx = 1.0;
	cx = 10;
	L.function = &weightedLik;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrap, bx, params);

	//check if there is a minimum, if not return 666
	if (fa == fb || fc == fb || ax > cx){
	  min = fa;
	  if(fb < min){
	    min = fb;
	  }
	  if(fc < min){
	    min = fc;
	  }
	  *lik_beta_hat = min;

	  return 666.0;
	}
	else{
	  // do mimization

	  //initialize lik function
	  L.function = &weightedLik;
	  L.params = params;

	  //min routine
	  T = gsl_min_fminimizer_brent;
	  s = gsl_min_fminimizer_alloc(T);
	  gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	  do{
	    iter++;
	    status = gsl_min_fminimizer_iterate(s);
	    m = gsl_min_fminimizer_x_minimum (s);
	    a = gsl_min_fminimizer_x_lower (s);
	    b = gsl_min_fminimizer_x_upper (s);
	    status = gsl_min_test_interval (a, b, 0.001, 0.0);
	  }
	  while (status == GSL_CONTINUE && iter < max_iter);
	  *lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	  gsl_min_fminimizer_free(s);
	  return m;
	}
}

double weighted_ml_est_lookup(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 1000;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20000;
	L.function = &weightedLikLook;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLook, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}

	if (fb == fc){
		return HUGE_VAL;
	}
	// do mimization

	//initialize lik function
	L.function = &weightedLikLook;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.01, 0.01);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* same as above but outputs 666 on minimization error */
double weighted_ml_est_lookup_errReport(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20;
	L.function = &weightedLikLook;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLook, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}
	//check for error
	if (fb >= fc || fb >= fa){
	  return(666.0);
	}

	// do mimization

	//initialize lik function
	L.function = &weightedLikLook;
	L.params = params;

	//min routine
	s = params->min;
	gsl_set_error_handler_off ();
	status = gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	if (status == GSL_EINVAL){
	  return(666);
	}
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.01, 0.01);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	//gsl_min_fminimizer_free(s);
	return m;
}

/* get weighted ML lower CI beta using weightedLikLookCI method; BRENT method for gsl */
double weighted_ml_CILower_lookup(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	double m, a, b, ax, cx;
	gsl_function L;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	// get minimum bracket
	ax = -50.0;
	cx = params->beta_hat;

	// do root finding
	//initialize lik function
	L.function = &weightedLikLookCI;
	L.params = params;

	//min routine
	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &L, ax, cx);
	do{
	  iter++;
	  status = gsl_root_fsolver_iterate(s);
	  m = gsl_root_fsolver_root(s);
	  a = gsl_root_fsolver_x_lower(s);
	  b = gsl_root_fsolver_x_upper(s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	gsl_root_fsolver_free(s);
	return m;
}

/* get weighted ML upper CI beta using weightedLikLookCI method; BRENT method for gsl */
double weighted_ml_CIUpper_lookup(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	double m, a, b, ax, cx;
	gsl_function L;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	// get minimum bracket
	cx = 100.0;
	ax = params->beta_hat;

	// do root finding
	//initialize lik function
	L.function = &weightedLikLookCI;
	L.params = params;

	//min routine
	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &L, ax, cx);
	do{
	  iter++;
	  status = gsl_root_fsolver_iterate(s);
	  m = gsl_root_fsolver_root(s);
	  a = gsl_root_fsolver_x_lower(s);
	  b = gsl_root_fsolver_x_upper(s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	gsl_root_fsolver_free(s);
	return m;
}

//Folded Version
double weighted_ml_est_lookup_folded(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 1000;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20000;
	L.function = &weightedLikLook;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLookFolded, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}

	if (fb == fc){
		return HUGE_VAL;
	}
	// do mimization

	//initialize lik function
	L.function = &weightedLikLook;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.01, 0.01);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* get weighted ML lower CI beta using weightedLikLookCI method; BRENT method for gsl */
double weighted_ml_CILower_lookupFolded(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	double m, a, b, ax, cx;
	gsl_function L;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	// get minimum bracket
	ax = -50.0;
	cx = params->beta_hat;

	// do root finding
	//initialize lik function
	L.function = &weightedLikLookCIFolded;
	L.params = params;

	//min routine
	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &L, ax, cx);
	do{
	  iter++;
	  status = gsl_root_fsolver_iterate(s);
	  m = gsl_root_fsolver_root(s);
	  a = gsl_root_fsolver_x_lower(s);
	  b = gsl_root_fsolver_x_upper(s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	gsl_root_fsolver_free(s);
	return m;
}

/* get weighted ML upper CI beta using weightedLikLookCI method; BRENT method for gsl */
double weighted_ml_CIUpper_lookupFolded(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	double m, a, b, ax, cx;
	gsl_function L;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	// get minimum bracket
	cx = 100.0;
	ax = params->beta_hat;

	// do root finding
	//initialize lik function
	L.function = &weightedLikLookCIFolded;
	L.params = params;

	//min routine
	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &L, ax, cx);
	do{
	  iter++;
	  status = gsl_root_fsolver_iterate(s);
	  m = gsl_root_fsolver_root(s);
	  a = gsl_root_fsolver_x_lower(s);
	  b = gsl_root_fsolver_x_upper(s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	gsl_root_fsolver_free(s);
	return m;
}

//ML estimation using PRF likelihood model
double ml_prf(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct prfLik_params * params = (struct prfLik_params *) p;

	// get minimum bracket
	ax = -10.0;
	bx = 0.0;
	cx = 10;

	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &prfLikWrap, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}

	if (fb == fc){
		return HUGE_VAL;
	}
	// do mimization

	//initialize lik function
	L.function = &prfLik;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* summarizeSFS- returns a matrix of counts which summarize the SFS conditional
   upon sampleSize. rows = sampleSize, column = freqs */

gsl_matrix *summarizeSFS(int maxSampleSize, struct snp data[], int snpNumber){
  gsl_matrix *counts;
  int i;

  //alloc matrix of size maxSampleSize by maxSampleSize
  counts = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize + 1);
  gsl_matrix_set_zero(counts);
  assert(counts != NULL);

  //gsl_matrix_set_zero(counts);

  //go through data, fill up matrix
  for(i = 0; i < snpNumber; i++){
    gsl_matrix_set(counts, data[i].n, data[i].i, gsl_matrix_get(counts, data[i].n, data[i].i) + 1);
  }
  return(counts);

}

/* summarizeSFSFill- fills a matrix of counts which summarize the SFS conditional
   upon sampleSize. rows = sampleSize, column = freqs */

void summarizeSFSFill(gsl_matrix *counts, struct snp data[], int snpNumber){
  int i;

  //fill matrix of size maxSampleSize by maxSampleSize
  gsl_matrix_set_zero(counts);
  assert(counts != NULL);
  gsl_matrix_set_zero(counts);

  //go through data, fill up matrix
  for(i = 0; i < snpNumber; i++){
    gsl_matrix_set(counts, data[i].n, data[i].i, gsl_matrix_get(counts, data[i].n, data[i].i) + 1);
  }
}

/* summarizeSFSFill- fills a matrix of counts which summarize the SFS conditional
   upon sampleSize. rows = sampleSize, column = freqs */

int summarizeSFSFillFromTo(gsl_matrix *counts, struct snp data[], int snpNumber, int from, int to, int startPos){
  int i;

  //fill matrix of size maxSampleSize by maxSampleSize
  gsl_matrix_set_zero(counts);
  assert(counts != NULL);
  gsl_matrix_set_zero(counts);

  //go through data, fill up matrix
	i = startPos;
	while(data[i].pos < to){
		if(data[i].pos >= from)
    		gsl_matrix_set(counts, data[i].n, data[i].i, gsl_matrix_get(counts, data[i].n, data[i].i) + 1);
		i++;
  	}
	return(i);
}

/* summarizeSFS- returns a matrix of booleans which summarize the SFS conditional
   upon sampleSize. rows = sampleSize, column = freqs */

gsl_matrix *summarizeSFSBool(int maxSampleSize, struct snp data[], int snpNumber){
  gsl_matrix *counts;
  int i;

  //alloc matrix of size maxSampleSize by maxSampleSize
  counts = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize + 1);
  gsl_matrix_set_zero(counts);
  assert(counts != NULL);

  gsl_matrix_set_zero(counts);

  //go through data, fill up matrix
  for(i = 0; i < snpNumber; i++){
    gsl_matrix_set(counts, data[i].n, data[i].i, 1);
  }
  return(counts);

}

/* summarizeSFSBoolFill- fills a matrix of booleans which summarize the SFS conditional
   upon sampleSize. rows = sampleSize, column = freqs */

void summarizeSFSBoolFill(gsl_matrix *counts, struct snp data[], int snpNumber){
  int i;

  //fill matrix of size maxSampleSize by maxSampleSize
  //counts = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize + 1);
  gsl_matrix_set_zero(counts);
  assert(counts != NULL);

  gsl_matrix_set_zero(counts);

  //go through data, fill up matrix
  for(i = 0; i < snpNumber; i++){
    gsl_matrix_set(counts, data[i].n, data[i].i, 1);
  }
}

/* simulateSiteFreq- returns an int representing the frequency of
   of a site (unfolded), i, in sampleSize n, conditional on sel.
   coeff. beta , also needs a pointer to a gsl random number generator */

int simulateSiteFreq(int sampleSize, double alpha, void *r){
  int i;
  double sum, probs[sampleSize - 1], rand;
  gsl_function p;
  struct my_F_params simParams;
  gsl_rng * rn = (gsl_rng *) r;

  simParams.n = sampleSize;
  /* put probs in vector */
  for(i = 0; i < sampleSize - 1; i++){
    simParams.i = i + 1;
    p.function = &snpProb;
    p.params = &simParams;
    probs[i] = GSL_FN_EVAL(&p, alpha);
  }

  /* choose frequency */
  rand = gsl_rng_uniform(rn);
  sum = 0;
  for(i = 0; i < sampleSize - 1; i++){
    sum += probs[i];
    if (rand <= sum){
      return(i+1);
    }
  }
  return(666);
}

/* ascertainment stuff -- following Nielsen et al. 2004  */

/* simpleAscertain--  Nielsen et al. 2004 eqn 2 *, returns the prob of ascertainment
given a sample freq */
double simpleAscertain(double sampleFreq, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  double sampleSize, ascSize,tmpSampleFreq;
  double num1 = 0.0;
  double num2 = 0.0;

  sampleSize = params->n;
  ascSize = params->ascSize;
  tmpSampleFreq = (int) sampleFreq;

  //test for weirdness

  if (sampleFreq < ascSize){
    num1 = 0;
  }
  else{
    num1 = gsl_sf_choose(sampleFreq, ascSize);
  }

  if ((sampleSize -sampleFreq) < ascSize){
    num2 = 0;
  }
  else{
    num2 =  gsl_sf_choose(sampleSize - sampleFreq, ascSize);
  }

  return(1.0 - ((num1 + num2) / gsl_sf_choose(sampleSize, ascSize)));
}

/* outgroupAscertain-- equivalent to Nielsen et al. 2004 eqn 2,
 returns the prob of ascertainment given a sample freq, but
is for use when ascertainment was based on single outgroup comparisons
(e.g. divergence between reference genomes). specifically this
calculates the prob of sampling only ancestral alleles.
IMPORTANT- sampleFreq here is derived allele freq */
double outgroupAscertain(double sampleFreq, void * p){
	struct my_F_params * params = (struct my_F_params *) p;
	int sampleSize, ascSize, derived;
	double num = 0.0;

	sampleSize = params->n;
	ascSize = params->ascSize;
	derived = params->derivedFlag;
	if(derived){
	    num =  xchoosey( sampleFreq, ascSize);
	  }
	  else{
	    num =  xchoosey((sampleSize - sampleFreq), ascSize);
	  }
	  return(num /xchoosey(sampleSize, ascSize));
}

/*same as above for case where ascSamp is not subset of sample */
double outgroupAscertain2(double sampleFreq, void * p){
	struct my_F_params * params = (struct my_F_params *) p;
	int sampleSize, ascSize,ascFreq, derived, N, j;
	double num = 0.0, denom = 0.0;

	sampleSize = params->n;
	ascSize = params->ascSize;
	ascFreq = params->ascFreq;
	derived = params->derivedFlag;

	N = sampleSize + ascSize;
	j = sampleFreq + ascFreq;
  	//was ascertainment for derived?
  	if(derived == 1){
		num =  xchoosey(j,sampleFreq) * xchoosey((N-j), sampleSize - sampleFreq);
  	}
  	else{
		num =  xchoosey((N-j), (sampleSize - sampleFreq)) * xchoosey(j,  sampleFreq);
  	}
//	printf("num: %f\n", num );
//	printf("N: %d j: %d sampleSize: %d ascSize: %d derivedFlag: %d\n",N,j, sampleSize,ascSize,derived);
	denom = xchoosey(N, sampleSize);
	if (denom == 0 || num == 0){
		return(0);
	}
//	printf("oaVal: %f\n",num)/denom;
  	return(num / denom);
}

/* probAscertainmentGivenModel-- Nielsen et al. 2005 eqn 8
   returns the prob of ascertainment given the model */
double probAscertainmentGivenModel(double beta, void *p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_F_params tempParams;

  long double tot, tmpres, tmpres2;
  int i=1;
  gsl_function F, G;

  /* add up total prob space */
  tot = 0.0;
  tempParams.n = params->n;
  tempParams.ascSize = params->ascSize;

  for(i = 1; i < params->n ; i++){
    tempParams.i = i;

    F.function =  &(snpProb);
    F.params = &tempParams;
    tmpres = 0.0;
    tmpres = GSL_FN_EVAL(&F, beta);
    G.function = &simpleAscertain;
    G.params = &tempParams;
    tmpres2 = GSL_FN_EVAL(&G, i) ;
    tot += tmpres * tmpres2;

  }

  return(tot) ;
}

/* probAscertainmentGivenModel-- Nielsen et al. 2005 eqn 8
   returns the prob of ascertainment given the model */
double probOutgroupAscertainmentGivenModel(double beta, void *p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_F_params tempParams;

  long double tot, tmpres, tmpres2;
  int i=1;
  gsl_function F, G;

  /* add up total prob space */
  tot = 0.0;
  tempParams.n = params->n;
  tempParams.ascSize = params->ascSize;

  for(i = 1; i < params->n ; i++){
    tempParams.i = i;

    F.function =  &(snpProb);
    F.params = &tempParams;
    tmpres = 0.0;
    tmpres = GSL_FN_EVAL(&F, beta);
    G.function = &outgroupAscertain;
    G.params = &tempParams;
    tmpres2 = GSL_FN_EVAL(&G, i) ;
    tot += tmpres * tmpres2;

  }

  return(tot) ;
}

/* probAscertainmentGivenModelLookup-- Nielsen et al. 2005 eqn 8
   returns the prob of ascertainment given the model, lookup version
   requires snpProb matrix and ascProb matrix */
double probAscertainmentGivenModelLookup(double beta, int sampSize, gsl_matrix *snpProbs, gsl_matrix *ascProbs){

  long double tot;
  int i=1;

  /* add up total prob space */
  tot = 0.0;

  for(i = 1; i < sampSize ; i++){
    tot += gsl_matrix_get(snpProbs, sampSize, i) * gsl_matrix_get(ascProbs, sampSize, i);
  }
  return(tot) ;
}

/* probAscertainmentGivenModelHemiLookup-- same
   as above but expects vector not matrix of snpProbs */
double probAscertainmentGivenModelHemiLookup(double beta, int sampSize, gsl_vector *snpProbs, gsl_matrix *ascProbs){

  long double tot;
  int i=1;

  /* add up total prob space */
  tot = 0.0;

  for(i = 1; i < sampSize ; i++){
    tot += gsl_vector_get(snpProbs, i) * gsl_matrix_get(ascProbs, sampSize, i);
  }
  return(tot) ;
}

/* ascSize vector - returns a vector of 1s or 0s depending on ascertainment sample sizes in the data */
gsl_vector *ascSizeVector(struct snp data[], int snpNumber, int maxSampleSize){
  gsl_vector *bools;
  int i;

  //alloc and initialize bools vector
  bools = gsl_vector_alloc(maxSampleSize + 1);
  gsl_vector_set_zero(bools);

  //go through data, set each sampleSize value to 1
  for(i = 0; i < snpNumber; i++){
    gsl_vector_set(bools,data[i].ascSize, 1);
  }
  return(bools);
}

/* snpAscMatrix- returns a matrix of simpleAscertainment probs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1). Currently this only supports a single ascertainment size  */
gsl_matrix *snpAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools, \
		 int ascSize){
  gsl_matrix *probs;
  struct my_F_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//have freq j?
	if(gsl_matrix_get(sfsBools, i, j)){
	  //calc prob
	  FParams.i = j;
	  FParams.n = i;
	  FParams.ascSize = ascSize;
	  f.function = &simpleAscertain;
	  f.params = &FParams;
	  gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, j));
	}
      }
    }
  }
  return(probs);
}

/* snpOutgroupAscMatrix- returns a matrix of outgroupAscertainment probs.
matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1).
Currently this only supports a single ascertainment size. Also has inelegance
of calculating derived or anscertral matrix  */
gsl_matrix *snpOutgroupAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, int ascSize, int derived, \
int runMode){
	gsl_matrix *probs;
	struct my_F_params FParams;
	gsl_function f;
	int i, j;

//alloc probs
	probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
	gsl_matrix_set_zero(probs);
//go through sampleSizes
	for(i = 2; i <= maxSampleSize; i++){
	//have sample size i?
		if (gsl_vector_get(sampleSizeVector, i)){
	//go through freqs
			for(j = 1; j < i; j++){
	//calc prob
				FParams.i = j;
				FParams.n = i;
				FParams.ascSize = ascSize;
				FParams.ascFreq = derived;
				FParams.derivedFlag = derived;
				if(runMode == 0){
					f.function = &outgroupAscertain;
				}
				else{
					f.function = &outgroupAscertain2;
				}
				f.params = &FParams;
				gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, j));
			}
		}
	}
	return(probs);
}

gsl_matrix *snpOutgroupAscMatrixPre(int maxSampleSize, gsl_vector *sampleSizeVector, int ascSize, int derived, \
int runMode, gsl_matrix *probs){
	struct my_F_params FParams;
	gsl_function f;
	int i, j;

	gsl_matrix_set_zero(probs);
//go through sampleSizes
	for(i = 2; i <= maxSampleSize; i++){
	//have sample size i?
		if (gsl_vector_get(sampleSizeVector, i)){
	//go through freqs
			for(j = 1; j < i; j++){
	//calc prob
				FParams.i = j;
				FParams.n = i;
				FParams.ascSize = ascSize;
				FParams.ascFreq = derived;
				FParams.derivedFlag = derived;
				if(runMode == 0){
					f.function = &outgroupAscertain;
				}
				else{
					f.function = &outgroupAscertain2;
				}
				f.params = &FParams;
				gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, j));
			}
		}
	}
	return(probs);
}

/*estimateAscSFS-- ML estimates of SFS, given ascertainment. returns
 a matrix of probs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
variable sampleSizes (2 to max), columns represent freqs (1 to max -1). Currently this only supports a single ascertainment size */

gsl_matrix *estimateAscSFS(int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools, \
			   int ascSize, gsl_matrix *sfsSummary){
  gsl_matrix *probs, *ascProbs;
    int i, j;
  double ssTmp, p_hat;
  double xij;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //calculate ascProbs
  ascProbs = snpAscMatrix(maxSampleSize, sampleSizeVector, sfsBools, \
			  ascSize);

  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs to tally up denom at sampleSize
      ssTmp = 0.0;
      for(j = 1; j < maxSampleSize; j++){
	xij = gsl_matrix_get(sfsSummary, i, j);
	if (xij > 0){
	  ssTmp += xij / gsl_matrix_get(ascProbs, i, j);
	}
      }
      //go back through freqs to assign p_hats at sampleSize
      for(j = 1; j < maxSampleSize; j++){
	xij = gsl_matrix_get(sfsSummary, i, j);
	p_hat = (double) xij / gsl_matrix_get(ascProbs, i, j);
	gsl_matrix_set(probs, i, j, p_hat / ssTmp);
      }
    }
  }
  return(probs);
}

/* weighted likelihood function for SFS given the data and ascertainment; here contains the data and the weights and uses the snpProbLookup routine */
double weightedLikLookAsc(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  double lik, num;
  int i;
  gsl_matrix *probs, *ascProbs;
  gsl_vector *ascDenoms;

  //make prob matrix (note this is in logs)
  probs = snpProbMatrix(beta, params->maxSampleSize, params->sampleSizeVector, params->sfsBools);
  //make ascProb matrix (NOT in logs!)
  ascProbs = snpAscMatrix(params->maxSampleSize, params->sampleSizeVector, params->sfsBools, params->ascSize);

  //make vector of denoms
  ascDenoms = gsl_vector_alloc(params->maxSampleSize + 1);
  gsl_vector_set_zero(ascDenoms);
  for(i = params->ascSize; i <= params->maxSampleSize; i++){
     gsl_vector_set(ascDenoms,i,probAscertainmentGivenModelLookup(beta, i, probs, ascProbs));
  }

  //now go and tally likelihood
  lik = 0;
  num = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)	\
      + log(gsl_matrix_get(ascProbs,params->data[i].n,params->data[i].i)) \
      - log(gsl_vector_get(ascDenoms, params->data[i].n));
     }
  gsl_matrix_free(probs);
  gsl_matrix_free(ascProbs);
  gsl_vector_free(ascDenoms);

  return -lik;
}

/* weighted likelihood function for SFS given the data and outgroup ascertainment; here contains the data and the weights and uses the snpProbLookup routine */
double weightedLikLookOutgroupAsc(double beta, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	double lik, num;
	int i;
	gsl_matrix *probs, *ascProbsDerived, *ascProbsAncestral;
	gsl_vector *ascDenomsDerived, *ascDenomsAncestral;

	//make prob matrix
	probs = snpProbMatrixNotLogFull(beta, params->maxSampleSize, params->sampleSizeVector);
	//make ascProb matrices, one for derived one for ancestral
	ascProbsAncestral = snpOutgroupAscMatrix(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 0, params->runMode);

	ascProbsDerived = snpOutgroupAscMatrix(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 1, params->runMode);

	//make vectors of denoms- derived, ancestral
	ascDenomsDerived = gsl_vector_alloc(params->maxSampleSize + 1);
	ascDenomsAncestral = gsl_vector_alloc(params->maxSampleSize + 1);
	gsl_vector_set_zero(ascDenomsDerived);
	gsl_vector_set_zero(ascDenomsAncestral);
	for(i = params->ascSize; i <= params->maxSampleSize; i++){
		gsl_vector_set(ascDenomsDerived,i,probAscertainmentGivenModelLookup(beta, i, probs, ascProbsDerived));
		gsl_vector_set(ascDenomsAncestral,i,probAscertainmentGivenModelLookup(beta, i, probs, ascProbsAncestral));
	}

//now go and tally likelihood
	lik = 0;
	num = 0;
	for(i = 0; i < params->snpNumber; i++){
		if(params->data[i].derivedFlag == 1){
			num = (gsl_matrix_get(probs,params->data[i].n,params->data[i].i) * gsl_matrix_get(ascProbsDerived,params->data[i].n,params->data[i].i)) \
				/ gsl_vector_get(ascDenomsDerived, params->data[i].n);
			}
		else{
			num = (gsl_matrix_get(probs,params->data[i].n,params->data[i].i) * gsl_matrix_get(ascProbsAncestral,params->data[i].n,params->data[i].i)) \
				/ gsl_vector_get(ascDenomsAncestral, params->data[i].n);
		}
		lik += log(num);
	}

	gsl_matrix_free(probs);
	gsl_matrix_free(ascProbsDerived);
	gsl_matrix_free(ascProbsAncestral);
	gsl_vector_free(ascDenomsDerived);
	gsl_vector_free(ascDenomsAncestral);
	return -lik;
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrapLookAsc(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLikLookAsc;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrapLookOutgroupAsc(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLikLookOutgroupAsc;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/* get weighted ML estimate of beta using likelihood function (lookup) conditional on ascertainment
; BRENT method for gsl */
double weighted_ml_est_lookup_asc(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20;
	L.function = &weightedLikLookAsc;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLookAsc, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}
	// do mimization

	//initialize lik function
	L.function = &weightedLikLookAsc;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* get weighted ML estimate of beta using likelihood function (lookup) conditional on outgroup ascertainment
; BRENT method for gsl */
double weighted_ml_est_lookup_outgroup_asc(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20;
	L.function = &weightedLikLookOutgroupAsc;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLookOutgroupAsc, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}
	// do mimization

	//initialize lik function
	L.function = &weightedLikLookOutgroupAsc;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* likelihood function for SFS where each SNP has
 independent selection coeff and outgroup ascertainment */
double sfsLikBetaVectorOutgroupAsc(gsl_vector *betas, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	double pSNP, pAsc, pAscModel,lik, tmp;
	int i;
	gsl_matrix *ascProbsDerived, *ascProbsAncestral;
	gsl_vector *snpProbs;

	lik = 0;
//make ascProb matrix
	ascProbsAncestral = snpOutgroupAscMatrix(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 0, params->runMode);

	ascProbsDerived = snpOutgroupAscMatrix(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 1, params->runMode);
//go through snps
	for(i = 0; i < params->snpNumber; i++){
	//create vector of snpProbs, find pSNP_i
		snpProbs = snpProbVectorNotLog(gsl_vector_get(betas, i), params->data[i].n);
		pSNP = gsl_vector_get(snpProbs, params->data[i].i);

	//lookup pAsc_i
		if(params->data[i].derivedFlag){
			pAsc = gsl_matrix_get(ascProbsDerived, params->data[i].n, params->data[i].i);
			//get pAscModel
			pAscModel = probAscertainmentGivenModelHemiLookup(gsl_vector_get(betas, i),params->data[i].n, \
					snpProbs, ascProbsDerived);
		}
		else{
			pAsc = gsl_matrix_get(ascProbsAncestral, params->data[i].n, params->data[i].i);
			//get pAscModel
			pAscModel = probAscertainmentGivenModelHemiLookup(gsl_vector_get(betas, i),params->data[i].n, \
					snpProbs, ascProbsAncestral);
		}
		tmp = pSNP * pAsc / pAscModel;
		lik += log(tmp);
		gsl_vector_free(snpProbs);
	}
	gsl_matrix_free(ascProbsAncestral);
	gsl_matrix_free(ascProbsDerived);
	return -lik;
}

double sfsLikBetaVectorOutgroupAscPre(gsl_vector *betas, gsl_matrix *ascProbsAncestral,gsl_matrix *ascProbsDerived, gsl_vector *snpProbs, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	double pSNP, pAsc, pAscModel,lik, tmp;
	int i;

	lik = 0.0;
//make ascProb matrix
	ascProbsAncestral = snpOutgroupAscMatrixPre(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 0, params->runMode,ascProbsAncestral);

	ascProbsDerived = snpOutgroupAscMatrixPre(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 1, params->runMode,ascProbsDerived);
//go through snps
	for(i = 0; i < params->snpNumber; i++){
	//create vector of snpProbs, find pSNP_i
		snpProbs = snpProbVectorNotLogPre(gsl_vector_get(betas, i), params->data[i].n,snpProbs);
		pSNP = gsl_vector_get(snpProbs, params->data[i].i);

	//lookup pAsc_i
		if(params->data[i].derivedFlag){
			pAsc = gsl_matrix_get(ascProbsDerived, params->data[i].n, params->data[i].i);
			//get pAscModel
			pAscModel = probAscertainmentGivenModelHemiLookup(gsl_vector_get(betas, i),params->data[i].n, \
					snpProbs, ascProbsDerived);
		}
		else{
			pAsc = gsl_matrix_get(ascProbsAncestral, params->data[i].n, params->data[i].i);
			//get pAscModel
			pAscModel = probAscertainmentGivenModelHemiLookup(gsl_vector_get(betas, i),params->data[i].n, \
					snpProbs, ascProbsAncestral);
		}
		tmp = pSNP * pAsc / pAscModel;
		lik += log(tmp);
//		printf("%f\n",lik);
	}
	return(-lik);
}
//classical summary stats

//tajimas theta_pi
double theta_pi(struct snp *data, int sites){
	int i;
	double n,p, sum;

	sum = 0.0;
	for(i=0;i<sites;i++){
		n = (double) data[i].n;
		p = (double) data[i].i / n;
		sum += (n / (n-1.0)) * (2.0 * p * (1.0 - p));
	}

	return(sum);
}

//this simply counts size with 0 < i < sampleSize
int segSites(struct snp *data, int sites){
	int i, count = 0;
	//go thru sites
	for(i=0;i<sites;i++){
		//check if monomorphic ##note I'm doing zero and n
		if  (data[i].i != 0 && data[i].i != data[i].n)
			count++;
	}
	return(count);
}
/*
double theta_w(int segSites, int n){
	return((double)segSites / a1f(n));
}

//harmonic sum for Waterson's theta
double a1f(int n){
	int i;
	double sum = 0;

	for(i=1;i<n;i++){
		sum += 1.0 / (double) i;
	}
	return(sum);
}
/*
//tajd
double tajd(int nsam, int segsites, double sumk){
  double  a1, a2, b1, b2, c1, c2, e1, e2;

  if( segsites == 0 ) return( 0.0) ;

  a1 = a1f(nsam);
  a2 = a2f(nsam);
  b1 = b1f(nsam);
  b2 = b2f(nsam);
  c1 = c1f(a1, b1);
  c2 = c2f(nsam, a1, a2, b2);
  e1 = e1f(a1, c1);
  e2 = e2f(a1, a2, c2);

  return( (sumk - (segsites/a1))/sqrt((e1*segsites) + ((e2*segsites)*(segsites
								      -1))) ) ;
}

double a2f(int nsam) {
  double a2;
  int i;
  a2 = 0.0;
  for (i=1; i<=nsam-1; i++) a2 += 1.0/(i*i);
  return (a2);
}

double b1f(int nsam){
  double b1;
  b1 = (nsam + 1.0)/(3.0*(nsam-1.0));
  return (b1);
}

double b2f(int nsam){
  double b2;
  b2 = (2*(nsam*nsam + nsam + 3.0))/(9*nsam*(nsam - 1));
  return (b2);
}

double e1f(double a1f, double c1){
  double e1;
  e1 = c1/a1f;
  return (e1);
}

double e2f(double a1f, double a2, double c2){
  double e2;
  e2 = c2/((a1f*a1f)+a2);
  return (e2);
}

double c1f(double a1f, double b1){
  double c1;
  c1 = b1 - (1/a1f);
  return (c1);
}

double c2f(int nsam, double a1f, double a2, double b2){
  double c2;
  c2 = b2 - ((nsam+2)/(a1f*nsam)) + (a2/(a1f * a1f));
  return (c2);
}

*/

/////////
/////
////  MPI stuff

//MPI version of above for parallel processing
double sfsLikBetaVectorMPI(gsl_vector *betas, void * p, int myrank, int mprocs, MPI_Comm comm){
	struct my_lik_params * params = (struct my_lik_params *) p;
	struct my_F_params FParams;
	gsl_function f;
	double lik, final;
	int snpNumber = params->snpNumber;
	int i, start, end, window;

	//initialize; where am i?
	lik = 0.;
	window = ceil((float)snpNumber / mprocs);
	start = myrank * window;
	end = start + window;
	FParams.w = params->w;
	f.function = &snpProb;
	for(i = start; i < end && i < snpNumber; i++){
		FParams.i = params->data[i].i;
		FParams.n = params->data[i].n;

		f.params = &FParams;
	//	printf("here lik:%f beta:%f i: %d n: %d\n",lik,gsl_vector_get(betas, i),FParams.i,FParams.n);

		lik += log(GSL_FN_EVAL(&f, gsl_vector_get(betas, i)));
		}
	MPI_Reduce(&lik,&final,1, MPI_DOUBLE,MPI_SUM,0,comm);
	return -final;
}

double sfsLikBetaVectorOutgroupAscPreMPI(gsl_vector *betas, gsl_matrix *ascProbsAncestral,gsl_matrix *ascProbsDerived, \
	gsl_vector *snpProbs, void * p,int myrank, int mprocs, MPI_Comm comm){

	struct my_lik_params * params = (struct my_lik_params *) p;
	double pSNP, pAsc, pAscModel,lik, tmp, final;
	int i,start, end, window,snpNumber;

	//initialization; window for proc to process
	lik = 0;
	snpNumber = params->snpNumber;
	window = ceil((float)snpNumber / mprocs);
	start = myrank * window;
	end = start + window;

	//currently every process calculates the matrices below-- could be faster
	//make ascProb matrix
	ascProbsAncestral = snpOutgroupAscMatrixPre(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 0, params->runMode,ascProbsAncestral);

	ascProbsDerived = snpOutgroupAscMatrixPre(params->maxSampleSize, params->sampleSizeVector, \
		params->ascSize, 1, params->runMode,ascProbsDerived);

	//go through snps
	for(i = start; i < end && i < snpNumber; i++){
	//create vector of snpProbs, find pSNP_i
		snpProbs = snpProbVectorNotLogPre(gsl_vector_get(betas, i), params->data[i].n,snpProbs);
		pSNP = gsl_vector_get(snpProbs, params->data[i].i);

	//lookup pAsc_i
		if(params->data[i].derivedFlag){
			pAsc = gsl_matrix_get(ascProbsDerived, params->data[i].n, params->data[i].i);
			//get pAscModel
			pAscModel = probAscertainmentGivenModelHemiLookup(gsl_vector_get(betas, i),params->data[i].n, \
					snpProbs, ascProbsDerived);
		}
		else{
			pAsc = gsl_matrix_get(ascProbsAncestral, params->data[i].n, params->data[i].i);
			//get pAscModel
			pAscModel = probAscertainmentGivenModelHemiLookup(gsl_vector_get(betas, i),params->data[i].n, \
					snpProbs, ascProbsAncestral);
		}
		tmp = pSNP * pAsc / pAscModel;
		lik += log(tmp);
	}
	MPI_Reduce(&lik,&final,1, MPI_DOUBLE,MPI_SUM,0,comm);
	return -final;
}

max_stars_repo_name: andrewkern/pgLib

max_issues_repo_name: andrewkern/pgLib

max_forks_repo_name: andrewkern/pgLib


text:
#pragma once

#include "icode_generator.h"
#include <halley/data_structures/hash_map.h>
#include <gsl/gsl>
#include "halley/file/path.h"

namespace YAML
{
	class Node;
}

namespace Halley
{
	class ComponentSchema;
	class SystemSchema;
	class MessageSchema;
	class CustomTypeSchema;

	class Codegen
	{
		struct Stats
		{
			int written = 0;
			int skipped = 0;
			std::vector<Path> files;
		};

		static bool doNothing(float, String) { return true; }

	public:
		using ProgressReporter = std::function<bool(float, String)>;

		static void run(Path inDir, Path outDir);

		explicit Codegen(bool verbose = false);

		void loadSources(std::vector<std::pair<String, gsl::span<const gsl::byte>>> files, ProgressReporter progress = &doNothing);
		void validate(ProgressReporter progress = &doNothing);
		void process();
		bool writeFile(Path path, const char* data, size_t dataSize, bool stub) const;
		void writeFiles(Path directory, const CodeGenResult& files, Stats& stats) const;
		std::vector<Path> generateCode(Path directory, ProgressReporter progress = &doNothing);

	private:
		void addSource(String name, gsl::span<const gsl::byte> data);
		void addComponent(YAML::Node rootNode);
		void addSystem(YAML::Node rootNode);
		void addMessage(YAML::Node rootNode);
		void addType(YAML::Node rootNode);
		String getInclude(String typeName) const;

		bool verbose;
		HashMap<String, ComponentSchema> components;
		HashMap<String, SystemSchema> systems;
		HashMap<String, MessageSchema> messages;
		HashMap<String, CustomTypeSchema> types;
	};
}

ext: h

max_stars_repo_name: lye/halley

max_issues_repo_name: lye/halley

max_forks_repo_name: lye/halley


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

/* brent.c -- brent minimum finding algorithm */

#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_min.h>

#include "min.h"

typedef struct
  {
    double d, e, v, w;
    double f_v, f_w;
  }
brent_state_t;

static int brent_init (void *vstate, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper);
static int brent_iterate (void *vstate, gsl_function * f, double *x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper);

static int
brent_init (void *vstate, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper)
{
  brent_state_t *state = (brent_state_t *) vstate;

  const double golden = 0.3819660;      /* golden = (3 - sqrt(5))/2 */

  double v = x_lower + golden * (x_upper - x_lower);
  double w = v;

  double f_vw;

  x_minimum = 0 ;  /* avoid warnings about unused varibles */
  f_minimum = 0 ;
  f_lower = 0 ;
  f_upper = 0 ;

  state->v = v;
  state->w = w;

  state->d = 0;
  state->e = 0;

  SAFE_FUNC_CALL (f, v, &f_vw);

  state->f_v = f_vw;
  state->f_w = f_vw;

  return GSL_SUCCESS;
}

static int
brent_iterate (void *vstate, gsl_function * f, double *x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper)
{
  brent_state_t *state = (brent_state_t *) vstate;

  const double x_left = *x_lower;
  const double x_right = *x_upper;

  const double z = *x_minimum;
  double d = state->e;
  double e = state->d;
  double u, f_u;
  const double v = state->v;
  const double w = state->w;
  const double f_v = state->f_v;
  const double f_w = state->f_w;
  const double f_z = *f_minimum;

  const double golden = 0.3819660;      /* golden = (3 - sqrt(5))/2 */

  const double w_lower = (z - x_left);
  const double w_upper = (x_right - z);

  const double tolerance =  GSL_SQRT_DBL_EPSILON * fabs (z);

  double p = 0, q = 0, r = 0;

  const double midpoint = 0.5 * (x_left + x_right);

  if (fabs (e) > tolerance)
    {
      /* fit parabola */

      r = (z - w) * (f_z - f_v);
      q = (z - v) * (f_z - f_w);
      p = (z - v) * q - (z - w) * r;
      q = 2 * (q - r);

      if (q > 0)
        {
          p = -p;
        }
      else
        {
          q = -q;
        }

      r = e;
      e = d;
    }

  if (fabs (p) < fabs (0.5 * q * r) && p < q * w_lower && p < q * w_upper)
    {
      double t2 = 2 * tolerance ;

      d = p / q;
      u = z + d;

      if ((u - x_left) < t2 || (x_right - u) < t2)
        {
          d = (z < midpoint) ? tolerance : -tolerance ;
        }
    }
  else
    {
      e = (z < midpoint) ? x_right - z : -(z - x_left) ;
      d = golden * e;
    }

  if (fabs (d) >= tolerance)
    {
      u = z + d;
    }
  else
    {
      u = z + ((d > 0) ? tolerance : -tolerance) ;
    }

  state->e = e;
  state->d = d;

  SAFE_FUNC_CALL(f, u, &f_u);

  if (f_u <= f_z)
    {
      if (u < z)
        {
          *x_upper = z;
          *f_upper = f_z;
        }
      else
        {
          *x_lower = z;
          *f_lower = f_z;
        }

      state->v = w;
      state->f_v = f_w;
      state->w = z;
      state->f_w = f_z;
      *x_minimum = u;
      *f_minimum = f_u;
      return GSL_SUCCESS;
    }
  else
    {
      if (u < z)
        {
          *x_lower = u;
          *f_lower = f_u;
          return GSL_SUCCESS;
        }
      else
        {
          *x_upper = u;
          *f_upper = f_u;
          return GSL_SUCCESS;
        }

      if (f_u <= f_w || w == z)
        {
          state->v = w;
          state->f_v = f_w;
          state->w = u;
          state->f_w = f_u;
          return GSL_SUCCESS;
        }
      else if (f_u <= f_v || v == z || v == w)
        {
          state->v = u;
          state->f_v = f_u;
          return GSL_SUCCESS;
        }
    }

  return GSL_FAILURE;
}

static const gsl_min_fminimizer_type brent_type =
{"brent",                       /* name */
 sizeof (brent_state_t),
 &brent_init,
 &brent_iterate};

const gsl_min_fminimizer_type *gsl_min_fminimizer_brent = &brent_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:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_ger (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int M = 1;
   int N = 1;
   int lda = 1;
   float alpha = 1.0f;
   float A[] = { -0.515f };
   float X[] = { 0.611f };
   int incX = -1;
   float Y[] = { -0.082f };
   int incY = -1;
   float A_expected[] = { -0.565102f };
   cblas_sger(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "sger(case 1390)");
     }
   };
  };

  {
   int order = 102;
   int M = 1;
   int N = 1;
   int lda = 1;
   float alpha = 1.0f;
   float A[] = { -0.515f };
   float X[] = { 0.611f };
   int incX = -1;
   float Y[] = { -0.082f };
   int incY = -1;
   float A_expected[] = { -0.565102f };
   cblas_sger(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "sger(case 1391)");
     }
   };
  };

  {
   int order = 101;
   int M = 1;
   int N = 1;
   int lda = 1;
   double alpha = 1;
   double A[] = { -0.809 };
   double X[] = { -0.652 };
   int incX = -1;
   double Y[] = { 0.712 };
   int incY = -1;
   double A_expected[] = { -1.273224 };
   cblas_dger(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dger(case 1392)");
     }
   };
  };

  {
   int order = 102;
   int M = 1;
   int N = 1;
   int lda = 1;
   double alpha = 1;
   double A[] = { -0.809 };
   double X[] = { -0.652 };
   int incX = -1;
   double Y[] = { 0.712 };
   int incY = -1;
   double A_expected[] = { -1.273224 };
   cblas_dger(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dger(case 1393)");
     }
   };
  };

  {
   int order = 101;
   int M = 1;
   int N = 1;
   int lda = 1;
   float alpha[2] = {0.0f, 0.0f};
   float A[] = { -0.651f, 0.856f };
   float X[] = { -0.38f, -0.235f };
   int incX = -1;
   float Y[] = { -0.627f, 0.757f };
   int incY = -1;
   float A_expected[] = { -0.651f, 0.856f };
   cblas_cgeru(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cgeru(case 1394) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cgeru(case 1394) imag");
     };
   };
  };

  {
   int order = 101;
   int M = 1;
   int N = 1;
   int lda = 1;
   float alpha[2] = {0.0f, 0.0f};
   float A[] = { -0.651f, 0.856f };
   float X[] = { -0.38f, -0.235f };
   int incX = -1;
   float Y[] = { -0.627f, 0.757f };
   int incY = -1;
   float A_expected[] = { -0.651f, 0.856f };
   cblas_cgerc(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cgerc(case 1395) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cgerc(case 1395) imag");
     };
   };
  };

  {
   int order = 102;
   int M = 1;
   int N = 1;
   int lda = 1;
   float alpha[2] = {0.0f, 0.0f};
   float A[] = { -0.651f, 0.856f };
   float X[] = { -0.38f, -0.235f };
   int incX = -1;
   float Y[] = { -0.627f, 0.757f };
   int incY = -1;
   float A_expected[] = { -0.651f, 0.856f };
   cblas_cgeru(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cgeru(case 1396) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cgeru(case 1396) imag");
     };
   };
  };

  {
   int order = 102;
   int M = 1;
   int N = 1;
   int lda = 1;
   float alpha[2] = {0.0f, 0.0f};
   float A[] = { -0.651f, 0.856f };
   float X[] = { -0.38f, -0.235f };
   int incX = -1;
   float Y[] = { -0.627f, 0.757f };
   int incY = -1;
   float A_expected[] = { -0.651f, 0.856f };
   cblas_cgerc(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cgerc(case 1397) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cgerc(case 1397) imag");
     };
   };
  };

  {
   int order = 101;
   int M = 1;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-1, 0};
   double A[] = { -0.426, 0.757 };
   double X[] = { -0.579, -0.155 };
   int incX = -1;
   double Y[] = { 0.831, 0.035 };
   int incY = -1;
   double A_expected[] = { 0.049724, 0.90607 };
   cblas_zgeru(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zgeru(case 1398) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zgeru(case 1398) imag");
     };
   };
  };

  {
   int order = 101;
   int M = 1;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-1, 0};
   double A[] = { -0.426, 0.757 };
   double X[] = { -0.579, -0.155 };
   int incX = -1;
   double Y[] = { 0.831, 0.035 };
   int incY = -1;
   double A_expected[] = { 0.060574, 0.86554 };
   cblas_zgerc(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zgerc(case 1399) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zgerc(case 1399) imag");
     };
   };
  };

  {
   int order = 102;
   int M = 1;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-1, 0};
   double A[] = { -0.426, 0.757 };
   double X[] = { -0.579, -0.155 };
   int incX = -1;
   double Y[] = { 0.831, 0.035 };
   int incY = -1;
   double A_expected[] = { 0.049724, 0.90607 };
   cblas_zgeru(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zgeru(case 1400) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zgeru(case 1400) imag");
     };
   };
  };

  {
   int order = 102;
   int M = 1;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-1, 0};
   double A[] = { -0.426, 0.757 };
   double X[] = { -0.579, -0.155 };
   int incX = -1;
   double Y[] = { 0.831, 0.035 };
   int incY = -1;
   double A_expected[] = { 0.060574, 0.86554 };
   cblas_zgerc(order, M, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zgerc(case 1401) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zgerc(case 1401) imag");
     };
   };
  };

}


text:
/**
 *
 * @file qwrapper_zgemm_tile.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 * Version of zgemm for tile storage, to avoid dependency problem when
 * computations are done within the tile. alpha and beta are passed as
 * pointers so they can depend on runtime values.
 *
 * @param[in] Alock
 *          Pointer to tile owning submatrix A.
 *
 * @param[in] Block
 *          Pointer to tile owning submatrix B.
 *
 * @param[in] Clock
 *          Pointer to tile owning submatrix C.
 *
 **/
void QUARK_CORE_zgemm_tile(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, PLASMA_enum transB,
                           int m, int n, int k, int nb,
                           const PLASMA_Complex64_t *alpha, const PLASMA_Complex64_t *A, int lda,
                                                            const PLASMA_Complex64_t *B, int ldb,
                           const PLASMA_Complex64_t *beta,        PLASMA_Complex64_t *C, int ldc,
                           const PLASMA_Complex64_t *Alock,
                           const PLASMA_Complex64_t *Block,
                           const PLASMA_Complex64_t *Clock)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_zgemm_tile_quark, task_flags,
        sizeof(PLASMA_enum),              &transA, VALUE,
        sizeof(PLASMA_enum),              &transB, VALUE,
        sizeof(int),                      &m,      VALUE,
        sizeof(int),                      &n,      VALUE,
        sizeof(int),                      &k,      VALUE,
        sizeof(PLASMA_Complex64_t),       alpha,           INPUT,
        sizeof(PLASMA_Complex64_t)*nb*nb, A,               NODEP,          /* input; see Alock */
        sizeof(int),                      &lda,    VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb, B,               NODEP,          /* input; see Block */
        sizeof(int),                      &ldb,    VALUE,
        sizeof(PLASMA_Complex64_t),       beta,            INPUT,
        sizeof(PLASMA_Complex64_t)*nb*nb, C,                       NODEP,  /* inout; see Clock */
        sizeof(int),                      &ldc,    VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb, Alock,           INPUT,
        sizeof(PLASMA_Complex64_t)*nb,    Block,           INPUT,
        sizeof(PLASMA_Complex64_t)*nb,    Clock,                   INOUT,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgemm_tile_quark = PCORE_zgemm_tile_quark
#define CORE_zgemm_tile_quark PCORE_zgemm_tile_quark
#endif
void CORE_zgemm_tile_quark(Quark *quark)
{
    PLASMA_enum transA, transB;
    int m, n, k, lda, ldb, ldc;
    const PLASMA_Complex64_t *alpha, *beta;
    const PLASMA_Complex64_t *A, *B;
    PLASMA_Complex64_t *C;

    quark_unpack_args_13( quark, transA, transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc );
    cblas_zgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        m, n, k,
        CBLAS_SADDR(*alpha), A, lda,
                             B, ldb,
        CBLAS_SADDR(*beta),  C, ldc );
}


text:
/*
 file includeDef.h
 list of headers and definitions used by the program OmegaMaxEnt (main source files: OmegaMaxEnt_main.cpp, OmegaMaxEnt_data.h, OmegaMaxEnt_data.cpp, graph_2D.h, graph_2D.cpp, generique.h, generique.cpp)

 Copyright (C) 2015 Dominic Bergeron (dominic.bergeron@usherbrooke.ca)

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

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

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

#ifndef INCLUDEDEF_H
#define INCLUDEDEF_H

#include "fftw3.h"
#include <cmath>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdio>
#include <complex>
#include <string>
#include <cstring>
#include <ctime>
#include <limits>
#include <stdio.h>
#include <sys/stat.h>
#include <random>
#include <gsl/gsl_linalg.h>
#include "armadillo"

#ifndef PI
#define PI acos((double)-1.0)
#endif

#ifndef EPSILON
#define EPSILON numeric_limits<double>::epsilon()
#endif

#ifndef DBL_MIN
#define DBL_MIN numeric_limits<double>::min()
#endif

#ifndef DBL_MAX
#define DBL_MAX numeric_limits<double>::max()
#endif

#ifndef INF
#define INF numeric_limits<double>::infinity()
#endif

using namespace std;

typedef complex<double> dcomplex;
typedef unsigned int  uint;

#endif

max_stars_repo_name: amstremblay/OmegaMaxEnt

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: amstremblay/OmegaMaxEnt

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: amstremblay/OmegaMaxEnt

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
/*
  NAME:
     normalize_row
  PURPOSE:
     normalize a row (or column) of a matrix given as logs
  CALLING SEQUENCE:
     normalize_row(gsl_matrix * q, int row,bool isrow,bool noweight,
     double weight)
  INPUT:
     q            - matrix
     row          - row to be normalized
     isrow        - is it a row or a column
     noweight    - add a weight to all of the values?
     weight       - weight to be added to all of the values
  OUTPUT:
     normalization factor (i.e. logsum)
  REVISION HISTORY:
     2008-09-21 - Written Bovy
     2010-04-01 - Added noweight and weight inputs to allow the qij to have
                  weights - Bovy
*/
#include <math.h>
#include <gsl/gsl_matrix.h>
#include "proj_gauss_mixtures.h"

double normalize_row(gsl_matrix * q, int row, bool isrow,
		     bool noweight, double weight){
  double loglike;
  if (isrow)
    loglike = logsum(q,row,true);
  else
    loglike = logsum(q,row,false);

  int dd;
  if (isrow)
    for (dd = 0; dd != q->size2; ++dd) {
      if ( noweight )
	gsl_matrix_set(q,row,dd,gsl_matrix_get(q,row,dd)-loglike);
      else
	gsl_matrix_set(q,row,dd,gsl_matrix_get(q,row,dd)-loglike+weight);
    }
  else
    for (dd = 0; dd != q->size1; ++dd) {
      if ( noweight )
	gsl_matrix_set(q,dd,row,gsl_matrix_get(q,dd,row)-loglike);
      else
	gsl_matrix_set(q,dd,row,gsl_matrix_get(q,dd,row)-loglike+weight);
    }
  if ( ! noweight ) loglike*= exp(weight);

  return loglike;
}

max_stars_repo_name: surbut/mashr

max_issues_repo_name: surbut/mashr

max_forks_repo_name: surbut/mashr


text:
#pragma once

#include <list>

#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

#include "../../../Utils/Random.h"
#include "../../Expression.h"
#include "../RandomDistributionDirectory.h"

double rand_Gauss(const std::list<expression>& args, const Variables& vars){
    if (args.size() == 1){
        double sigma = args.front()->value(vars);
        return gsl_ran_gaussian(rng.get(), sigma);
    }
    return GSL_NAN;
}

double pdf_Gauss(double x, const std::list<expression>& args, const Variables& vars){
    if (args.size() == 1){
        double sigma = args.front()->value(vars);
        return gsl_ran_gaussian_pdf(x, sigma);
    }
    return GSL_NAN;
}

double cdf_P_Gauss(double x, const std::list<expression>& args, const Variables& vars){
    if (args.size() == 1){
        double sigma = args.front()->value(vars);
        return gsl_cdf_gaussian_P(x, sigma);
    }
    return GSL_NAN;
}

double cdf_Q_Gauss(double x, const std::list<expression>& args, const Variables& vars){
    if (args.size() == 1){
        double sigma = args.front()->value(vars);
        return gsl_cdf_gaussian_Q(x, sigma);
    }
    return GSL_NAN;
}

double cdf_Pinv_Gauss(double x, const std::list<expression>& args, const Variables& vars){
    if (args.size() == 1){
        double sigma = args.front()->value(vars);
        return gsl_cdf_gaussian_Pinv(x, sigma);
    }
    return GSL_NAN;
}

double cdf_Qinv_Gauss(double x, const std::list<expression>& args, const Variables& vars){
    if (args.size() == 1){
        double sigma = args.front()->value(vars);
        return gsl_cdf_gaussian_Qinv(x, sigma);
    }
    return GSL_NAN;
}

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once
#include <iosfwd>
#include <vector>
#include <algorithm>
#include <string>
#include <cstring>
#include <gsl/gsl>
#include "onnxruntime_config.h"

#ifndef DISABLE_ABSEIL
// Need to include abseil inlined_vector.h header directly here
// as hash tables cause CUDA 10.2 compilers to fail. inlined_vector.h is fine.
#ifdef _MSC_VER
#pragma warning(push)
// C4127: conditional expression is constant
#pragma warning(disable : 4127)
// C4324: structure was padded due to alignment specifier
// Usage of alignas causes some internal padding in places.
#pragma warning(disable : 4324)
#endif

#include <absl/container/inlined_vector.h>

#ifdef _MSC_VER
#pragma warning(pop)
#endif
#endif  // DISABLE_ABSEIL

namespace onnxruntime {
#ifdef __GNUC__
#pragma GCC diagnostic push
#ifdef HAS_NULL_DEREFERENCE
#pragma GCC diagnostic ignored "-Wnull-dereference"
#endif
#endif

constexpr size_t kTensorShapeSmallBufferElementsSize = 5;

#ifndef DISABLE_ABSEIL
// Use this type to build a shape and then create TensorShape.
using TensorShapeVector = absl::InlinedVector<int64_t, kTensorShapeSmallBufferElementsSize>;
#else
class TensorShapeVector : public std::vector<int64_t> {
  using Base = std::vector<int64_t>;

 public:
   using Base::Base;
};

#endif  // DISABLE_ABSEIL

inline TensorShapeVector ToShapeVector(const gsl::span<const int64_t>& span) {
  TensorShapeVector out;
  out.reserve(span.size());
  out.assign(span.begin(), span.end());
  return out;
}

inline gsl::span<const int64_t> ToConstSpan(const TensorShapeVector& vec) {
  return gsl::make_span(vec);
}

class TensorShape {
  // We use negative numbers for unknown symbolic dimension. Each negative
  // number represents a unique symbolic dimension.
 public:
  TensorShape() = default;

  TensorShape(const TensorShape& other) : TensorShape(other.GetDims()) {}
  TensorShape& operator=(const TensorShape& other);
  TensorShape& operator=(const gsl::span<const int64_t>& dims) {
    *this = TensorShape(dims);
    return *this;
  }

  TensorShape(TensorShape&& other) noexcept { operator=(std::move(other)); }
  TensorShape& operator=(TensorShape&& other) noexcept;

  TensorShape(gsl::span<const int64_t> dims);
  TensorShape(const TensorShapeVector& dims) : TensorShape(gsl::make_span(dims)) {}
  TensorShape(std::initializer_list<int64_t> dims) : TensorShape(gsl::make_span(dims.begin(), dims.end())) {}
  TensorShape(const int64_t* dimension_sizes, size_t dimension_count) : TensorShape(gsl::span<const int64_t>(dimension_sizes, dimension_count)) {}
  TensorShape(const std::vector<int64_t>& dims, size_t start, size_t end) : TensorShape(gsl::span<const int64_t>(&dims[start], end - start)) {}

  // Create a TensorShape that points to an existing buffer internally. As no copy is made, 'data' must remain valid for the life of the TensorShape
  static const TensorShape FromExistingBuffer(const std::vector<int64_t>& data) {
    return TensorShape(External{}, gsl::span<int64_t>(const_cast<int64_t*>(data.data()), data.size()));
  }

  /**
     Return the dimension specified by <idx>.
  */
  int64_t operator[](size_t idx) const { return values_[idx]; }
  int64_t& operator[](size_t idx) { return values_[idx]; }

  bool operator==(const TensorShape& other) const noexcept { return GetDims() == other.GetDims(); }
  bool operator!=(const TensorShape& other) const noexcept { return !(*this == other); }

  size_t NumDimensions() const noexcept {
    return values_.size();
  }

  /**
     Copy dims into an array with given size
  */
  void CopyDims(int64_t* dims, size_t num_dims) const {
    memcpy(dims, values_.data(), sizeof(int64_t) * std::min(num_dims, NumDimensions()));
  }

  /**
     Copy dims from a specific start dim into an array with given size
     `start_dim` is expected to be in the inclusive range [0, NumDimensions() - 1]
     and this function does no checks to ensure that
  */
  void CopyDims(int64_t* dims, size_t start_dim, size_t num_dims) const {
    memcpy(dims, values_.data() + start_dim, sizeof(int64_t) * std::min(num_dims, NumDimensions() - start_dim));
  }

  /**
     Return underlying vector representation.
  */
  gsl::span<const int64_t> GetDims() const { return values_; }

  TensorShapeVector AsShapeVector() const {
    return ToShapeVector(values_);
  }

  /**
   * Return the total number of elements. Returns 1 for an empty (rank 0) TensorShape.
   *
   * May return -1
   */
  int64_t Size() const;

  /**
     Return the total number of elements up to the specified dimension.
     If the dimension interval is empty (dimension == 0), return 1.
     @param dimension Return size up to this dimension. Value must be between 0 and this->NumDimensions(), inclusive.
  */
  int64_t SizeToDimension(size_t dimension) const;

  /**
     Return the total number of elements from the specified dimension to the end of the tensor shape.
     If the dimension interval is empty (dimension == this->NumDimensions()), return 1.
     @param dimension Return size from this dimension to the end. Value must be between 0 and this->NumDimensions(),
                      inclusive.
  */
  int64_t SizeFromDimension(size_t dimension) const;

  /**
     Return a new TensorShape of the dimensions from dimstart to dimend.
  */
  TensorShape Slice(size_t dimstart, size_t dimend) const;

  /**
     Return a new TensorShape of the dimensions from dimstart to end.
  */
  TensorShape Slice(size_t dimstart) const { return Slice(dimstart, values_.size()); }

  /**
     output dimensions nicely formatted
  */
  std::string ToString() const;

  /**
     Calculate size between start and end.
     Assumes start and end are between 0 and this->NumDimensions(), inclusive, and that
     start < end.
  */
  int64_t SizeHelper(size_t start, size_t end) const;

  /**
     empty shape or 1D shape (1) is regarded as scalar tensor
  */
  bool IsScalar() const {
    size_t len = values_.size();
    return len == 0 || (len == 1 && values_[0] == 1);
  }

 private:
  struct External {};
  TensorShape(External, gsl::span<int64_t> buffer) : values_{buffer} {}

  void Allocate(size_t size);

  gsl::span<int64_t> values_;
  int64_t small_buffer_[kTensorShapeSmallBufferElementsSize];
  std::unique_ptr<int64_t[]> allocated_buffer_;

  friend struct ProviderHostImpl;  // So that the shared provider interface can access Allocate
};
#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif
// operator<< to nicely output to a stream
std::ostream& operator<<(std::ostream& out, const TensorShape& shape);

}  // namespace onnxruntime

max_stars_repo_name: mszhanyi/onnxruntime

max_issues_repo_name: mszhanyi/onnxruntime

max_forks_repo_name: mszhanyi/onnxruntime


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_permutation.h>

gsl_permutation *
gsl_permutation_alloc (const size_t n)
{
  gsl_permutation * p;

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

  p = (gsl_permutation *) malloc (sizeof (gsl_permutation));

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

  p->data = (size_t *) malloc (n * sizeof (size_t));

  if (p->data == 0)
    {
      free (p);         /* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for permutation data",
                        GSL_ENOMEM, 0);
    }

  p->size = n;

  return p;
}

gsl_permutation *
gsl_permutation_calloc (const size_t n)
{
  size_t i;

  gsl_permutation * p =  gsl_permutation_alloc (n);

  if (p == 0)
    return 0;

  /* initialize permutation to identity */

  for (i = 0; i < n; i++)
    {
      p->data[i] = i;
    }

  return p;
}

void
gsl_permutation_init (gsl_permutation * p)
{
  const size_t n = p->size ;
  size_t i;

  /* initialize permutation to identity */

  for (i = 0; i < n; i++)
    {
      p->data[i] = i;
    }
}

void
gsl_permutation_free (gsl_permutation * p)
{
  free (p->data);
  free (p);
}


text:
#pragma once

#include "halley/net/connection/iconnection.h"
#include "halley/net/connection/network_packet.h"

#ifdef _MSC_VER
#pragma warning(disable: 4834)
#endif
#define BOOST_SYSTEM_NO_DEPRECATED
#define BOOST_ERROR_CODE_HEADER_ONLY
#include <boost/asio.hpp>

#include <deque>
#include <array>
#include <string>
#include <gsl/gsl>

namespace Halley
{
	class NetworkService;
	using UDPEndpoint = boost::asio::ip::udp::endpoint;
	using UDPSocket = boost::asio::ip::udp::socket;

	class AsioUDPConnection : public IConnection
	{
	public:
		AsioUDPConnection(UDPSocket& socket, UDPEndpoint remote);

		void close() override;
		ConnectionStatus getStatus() const override { return status; }
		void send(OutboundNetworkPacket&& packet) override;
		bool receive(InboundNetworkPacket& packet) override;

		bool matchesEndpoint(const UDPEndpoint& remoteEndpoint) const;
		void onReceive(gsl::span<const gsl::byte> data);
		void setError(const std::string& cs);

		void open(short connectionId);
		void onOpen(short connectionId);
		void terminateConnection();
		short getConnectionId() const { return connectionId; }

	private:
		UDPSocket& socket;
		UDPEndpoint remote;
		ConnectionStatus status;
		short connectionId;

		std::deque<OutboundNetworkPacket> pendingSend;
		std::deque<InboundNetworkPacket> pendingReceive;
		std::array<gsl::byte, 2048> sendBuffer;
		std::string error;

		void sendNext();
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
//! \file heat_driver.h
//! Driver for Magnolia's heat transfer solver
#ifndef ENRICO_SURROGATE_HEAT_DRIVER_H
#define ENRICO_SURROGATE_HEAT_DRIVER_H

#include "driver.h"

#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
#include <gsl/gsl>
#include <mpi.h>

#include <cstddef>

namespace enrico {

class SurrogateHeatDriver : public Driver {
public:
  //! Initializes heat-fluids surrogate with the given MPI communicator.
  //!
  //! \param comm  The MPI communicator used to initialze the surrogate
  //! \param node  XML node containing settings for surrogate
  explicit SurrogateHeatDriver(MPI_Comm comm, pugi::xml_node node);

  //! Solves the heat-fluids surrogate solver
  void solve_step() final;

  //! Number of rings in fuel and clad
  //! \return Number of rings
  std::size_t n_rings() { return n_fuel_rings_ + n_clad_rings_; }

  //! Write data to VTK
  void write_step(int timestep, int iteration) final;

  xt::xtensor<double, 1> temperature() const;

  double temperature(int pin, int axial, int ring) const;

  // Data on fuel pins
  xt::xtensor<double, 2> pin_centers_; //!< (x,y) values for center of fuel pins
  xt::xtensor<double, 1> z_;           //!< Bounding z-values for axial segments
  std::size_t n_pins_;                 //!< number of fuel pins
  std::size_t n_axial_;                //!< number of axial segments

  // Dimensions for a single fuel pin axial segment
  double clad_outer_radius_;     //!< clad outer radius in [cm]
  double clad_inner_radius_;     //!< clad inner radius in [cm]
  double pellet_radius_;         //!< fuel pellet radius in [cm]
  std::size_t n_fuel_rings_{20}; //!< number of fuel rings
  std::size_t n_clad_rings_{2};  //!< number of clad rings

  // solver variables and settings
  double tol_;                    //!< tolerance on convergence
  xt::xtensor<double, 3> source_; //!< heat source for each (axial segment, ring)
  xt::xtensor<double, 1> r_grid_clad_; //!< radii of each clad ring in [cm]
  xt::xtensor<double, 1> r_grid_fuel_; //!< radii of each fuel ring in [cm]

  // visualization
  std::string viz_basename_{
    "heat_surrogate"}; //!< base filename for visualization files (default: magnolia)
  std::string viz_iterations_{
    "none"};                    //!< visualization iterations to write (none, all, final)
  std::string viz_data_{"all"}; //!< visualization data to write
  std::string viz_regions_{"all"}; //!< visualization regions to write
  size_t vtk_radial_res_{20};      //!< radial resolution of resulting vtk files

private:
  //! Create internal arrays used for heat equation solver
  void generate_arrays();

  //!< temperature in [K] for each (axial segment, ring)
  xt::xtensor<double, 3> temperature_;

}; // end SurrogateHeatDriver

} // namespace enrico

#endif //ENRICO_SURROGATE_HEAT_DRIVER_H

max_stars_repo_name: aprilnovak/enrico

max_issues_repo_name: aprilnovak/enrico

max_forks_repo_name: aprilnovak/enrico


text:
/**
 *
 * @file qwrapper_cpotrf.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_clauum(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int n, int nb,
                       PLASMA_Complex32_t *A, int lda)
{
    DAG_CORE_LAUUM;
    QUARK_Insert_Task(quark, CORE_clauum_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &n,     VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    A,             INOUT,
        sizeof(int),                        &lda,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_clauum_quark = PCORE_clauum_quark
#define CORE_clauum_quark PCORE_clauum_quark
#endif
void CORE_clauum_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int N;
    PLASMA_Complex32_t *A;
    int LDA;

    quark_unpack_args_4(quark, uplo, N, A, LDA);
    LAPACKE_clauum_work(LAPACK_COL_MAJOR, lapack_const(uplo), N, A, LDA);
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_core.h>
#include <gbpCosmo_NFW_etc.h>
#include <gsl/gsl_sf_expint.h>

void set_NFW_params(double M, double z, int mode, cosmo_info **cosmo, double *c_vir, double *R_vir) {
    if(mode != NFW_MODE_DEFAULT)
        SID_exit_error("Unknown mode (%d) in set_NFW_params()", SID_ERROR_LOGIC, mode);

    switch(ADaPS_exist(*cosmo, "M_WDM")) {
        case GBP_FALSE: {
            double Omega_M  = ((double *)ADaPS_fetch(*cosmo, "Omega_M"))[0];
            double h_Hubble = ((double *)ADaPS_fetch(*cosmo, "h_Hubble"))[0];

            // Mass-concentration from Munoz-Cuartas et al 2010
            double w     = 0.029;
            double m     = 0.097;
            double alpha = -110.001;
            double beta  = 2469.720;
            double gamma = 16.885;
            double a_z   = w * z - m;
            double b_z   = alpha / (z + gamma) + beta / pow(z + gamma, 2.);
            double Delta = Delta_vir(z, *cosmo);
            Delta        = 200.;

            (*c_vir) = take_alog10(a_z * take_log10(M / (M_SOL / h_Hubble)) + b_z);
            (*R_vir) = R_Delta_z(M, Delta, z, *cosmo); // Bullock et al '01
        } break;
        case GBP_TRUE:
            SID_exit_error("ENS not working.", SID_ERROR_LOGIC);
            //(*c_vir)=c_ENS(M,z,*cosmo);          // Eke, Navarro and Steinmetz
            //(*R_vir)=R_Delta_z(M,200.,z,*cosmo); // R_200
            break;
    }
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode