text:
/* segSiteHMM v. 0.1 !!!
/
/
/ Andrew Kern 3/14/07
*/

#include "../hmm/hmm.h"
#include "../hmm/adkGSL.h"
#include "../hmm/popGenTools.h"
#include "segSites.h"
#include "segSiteHMM.h"
#include <math.h>
#include <assert.h>
#include <string.h>
#include <gsl/gsl_vector.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>

#define MAXSNPS 100000000

char modeFlag, *outfileName, smoothFlag;
int snpNumber, maxSampleSize, states, verboseFlag = 0, fixedEmisFlag = 0, sampleSize;
gsl_vector *sampleSizes, *piStart, *piEnd;
struct snp data[MAXSNPS];
gsl_matrix *m, *sfsCounts, *emisMatrix;
double gamma1, theta1;
gsl_vector *gammas;
double priorCons = 0;
double omega = 0;
int outputLevel = 0;

int main(int argc, char *argv[]){
  gsl_matrix *alphas, *betas, *posts, *logPosts;
	gsl_vector  *vits, *locs, *include, *include2;
	gsl_matrix  *emisLogs;
	double score1, sum, pFor;
	HMM *h;
	FILE *outfile;
	int i, j, **elements, elementCount, transitionPowerFlag;

	printf("segSiteHMM v0.3 -- thank you for playing thermonuclear war\n");
	printf("getting parameters/data\n");
	getParameters(argc,argv);

	//find max sampleSize
	maxSampleSize = 0;
	for(j = 0; j < snpNumber; j++){
	  if (data[j].n > maxSampleSize){
	    maxSampleSize = data[j].n;
	  }
	}
	//transition gaps?
	transitionPowerFlag = 0;

	//make sampleSize vector and summarize sfs
	sampleSizes = sampleSizeVector(data,snpNumber,maxSampleSize);
	sfsCounts = summarizeSFSBool(maxSampleSize, data, snpNumber);

	//was smoothing used?
	if(smoothFlag == 'f'){
	  //init with equal starting probs
	  sum = 1.0 / states;
	  for (i = 0; i < states; i++){
	    gsl_vector_set(piStart, i, sum);
	    gsl_vector_set(piEnd, i, sum);
	  }
	}
	else if(smoothFlag == 's'){
	  //set start and end probs based on smoothing
	  gsl_vector_set(piStart, 1, priorCons);
	  gsl_vector_set(piStart, 0, 1.0 - gsl_vector_get(piStart, 1));
	  gsl_vector_set(piEnd, 0, gsl_vector_get(piStart,0));
	  gsl_vector_set(piEnd, 1, gsl_vector_get(piStart,1));

	  //define transition probs based on omega and gamma
	  gsl_matrix_set(m, 1, 0, 1.0 / omega);
	  gsl_matrix_set(m, 1, 1, 1.0 - gsl_matrix_get(m, 1, 0));
	  gsl_matrix_set(m, 0, 1, (gsl_matrix_get(m, 1, 0) / gsl_vector_get(piStart, 0)) \
			 -  gsl_matrix_get(m, 1,0));
	  gsl_matrix_set(m, 0, 0, 1.0 - gsl_matrix_get(m, 0, 1));
	}

	printf("path is %d SNPs long\n",snpNumber);
	locs = gsl_vector_alloc(snpNumber);
	//fill up locs vector for use in calculating transition powers
	for(i = 0; i < snpNumber; i++){
	  gsl_vector_set(locs, i, data[i].pos);
	}

	//initialize HMM
	printf("initializing model\n");
	h = newHMM(m, piStart,piEnd,locs);
	initHMM(h);
	printHMM(h);

	switch(modeFlag){
	case 't' :
	  printf("training mode\nEM iterations follow:\n");
	  fflush(stdout);
	  emAlg(h, gammas,snpNumber, data, transitionPowerFlag);
	  //output selection coefficient
	  strcat(outfileName,".alphas");
	  outfile = fopen(outfileName, "w");
	  gsl_vector_fprintf(outfile,gammas,"%f");
	  fclose(outfile);
	  //output transition matrix
	  strcat(outfileName,".trainedTransMat");
	  printTransitions(h, outfileName);
	  break;
	case 'b' :
	  printf("train and decode mode\nEM iterations follow:\n");
	  fflush(stdout);
	  emAlg(h, gammas,snpNumber, data, transitionPowerFlag);
	  	if(outputLevel>0){
		//output selection coefficient
		  strcat(outfileName,".alphas");
		  outfile = fopen(outfileName, "w");
		  fprintf(outfile,"%f",gsl_vector_get(gammas,1));
		  fclose(outfile);
		  //output transition matrix
		  strcat(outfileName,".trainedTransMat");
		  printTransitions(h, outfileName);
		}
	  strcat(outfileName,".decoding");
	  printf("Viterbi decoding mode\nelement predicition will be stored in a file called \"%s\"\nPosterior Probs and Viterbi calls follow:\n",outfileName);
	  vits = gsl_vector_alloc(snpNumber);
          if(!fixedEmisFlag)
          {
	    emisLogs = gsl_matrix_alloc(h->nstates, snpNumber);
	    calculateEmissions(h, emisLogs, gammas, theta1, snpNumber, data);
          }
	  alphas = gsl_matrix_alloc(h->nstates,snpNumber);
	  betas = gsl_matrix_alloc(h->nstates,snpNumber);
	  posts = gsl_matrix_alloc(h->nstates, snpNumber);
	  logPosts =  gsl_matrix_alloc(h->nstates, snpNumber);
	  pFor = forwardAlg(h, alphas, emisLogs, snpNumber, transitionPowerFlag);
	  backwardAlg(h, betas, emisLogs, snpNumber, transitionPowerFlag);
	  posteriorProbsReduced(h, posts, pFor, alphas, betas, emisLogs, snpNumber, transitionPowerFlag);
	  viterbi(h, vits, emisLogs, snpNumber, transitionPowerFlag);

	  //print out likelihood
	  printf("Likelihood = %lf\n",pFor);

	if (outputLevel > 0){
		//output posteriors and viterbi calls
		outfile = fopen(outfileName, "w");
		for(i = 0; i < emisLogs->size2; i++){
			fprintf(outfile,"%d",(int)gsl_vector_get(h->otherData,i));
			for(j = 0; j < h->nstates; j++){
				fprintf(outfile,"\t%lf",gsl_matrix_get(posts, j, i));
			}
			fprintf(outfile,"\t%d\n",(int)gsl_vector_get(vits, i));
		}
		fclose(outfile);
  	}
	  //score elements
	  strcat(outfileName,".elements");
	  scoreElements(h, emisLogs, vits, outfileName, transitionPowerFlag);
	  break;
	case 'v' :
	  printf("Viterbi decoding mode\nelement predicition will be stored in a file called \"%s\"\nPosterior Probs and Viterbi calls follow:\n",outfileName);
	  vits = gsl_vector_alloc(snpNumber);
	  emisLogs = gsl_matrix_alloc(h->nstates, snpNumber);
          if (fixedEmisFlag)
          {
              initializeEmissions(emisLogs, emisMatrix, snpNumber, data);
          }
          else
          {
	      calculateEmissions(h, emisLogs, gammas, theta1, snpNumber, data);
          }
	  alphas = gsl_matrix_alloc(h->nstates,snpNumber);
	  betas = gsl_matrix_alloc(h->nstates,snpNumber);
	  posts = gsl_matrix_alloc(h->nstates, snpNumber);
	  logPosts =  gsl_matrix_alloc(h->nstates, snpNumber);

	  pFor = forwardAlg(h, alphas, emisLogs, snpNumber, transitionPowerFlag);
	  backwardAlg(h, betas, emisLogs, snpNumber, transitionPowerFlag);
	  posteriorProbsReduced(h, posts, pFor, alphas, betas, emisLogs, snpNumber, transitionPowerFlag);
	  viterbi(h, vits, emisLogs, snpNumber, transitionPowerFlag);

	  //print out likelihood
	  printf("Likelihood = %lf\n",pFor);

	if(outputLevel > 0){
	//output posteriors and viterbi calls
		outfile = fopen("snpPosts", "w");
		for(i = 0; i < emisLogs->size2; i++){
			fprintf(outfile,"%d",(int)gsl_vector_get(h->otherData,i));
			for(j = 0; j < h->nstates; j++){
				fprintf(outfile,"\t%lf",gsl_matrix_get(posts, j, i));
			}
			fprintf(outfile,"\t%d\n",(int)gsl_vector_get(vits, i));
		}
		fclose(outfile);
		scoreElements(h, emisLogs, vits, outfileName, transitionPowerFlag);
		break;
		}
	}

	freeHMM(h);
	return 0;
}

void usage(){
    fprintf(stderr,"usage: pgHMM - population genetic HMMs\nrun modes:\n\t-t train mode -- EM training for model\n\t-v viterbi decoding mode\n\t-b both train and decode\n\ninput control\n\t-S stateNumber\n\t-d dataFile\n\t-p transitionMatrix\n\t-s startProbs\n\t-e endProbs\n\t-g prior probability of being in selected state\n\t-E emissionMatrixFile sampleSize -- specifies a file containing an SFS for each state (including frequencies of monomorphic sites); entires must sum to 1, and data must have constant sampleSize\n\t-w omega smoothing parameter (note either transitionMatrix or gamma and omega need to be defined)\n\t-a beta1 beta2 -- value for selected classes\n\t-u 4Nu initial value\n\t-o outfile name (trained matrix or elements)\n\t-V verbose\n\t-O output level (default= 0)\n");
    exit(1);
}

/* sets the parameters/run mode and then opens data */
void getParameters(int argc, char *argv[]){
  FILE *infile;
  int  i, j, n, args;
  long int pos;

  if (argc < 2){
    usage();
  }
  else{
    args = 1;
    while(args < argc){
      switch(argv[args][1]){
      case 't' :
	/* train mode? */
	modeFlag = 't';
	smoothFlag = 'f';
	break;
      case 'v' :
	/* viterbi decode */
	modeFlag = 'v';
	smoothFlag = 'f';
	break;
      case 'b' :
	/*both train and decode */
	modeFlag = 'b';
	smoothFlag = 'f';
	break;
      case 'a' :
	/*set alphas for selected class */
	for(j = 1; j < states; j++){
	  gsl_vector_set(gammas,j,atof(argv[++args]));
	}
	break;
      case 'u' :
	//set theta
	theta1 = atof(argv[++args]);
	break;
	case 'O' :
	//set outputLevel
	outputLevel = atoi(argv[++args]);
	break;
      case 'S' :
	//set state number
	states = atoi(argv[++args]);
	m = gsl_matrix_alloc(states,states);
	gsl_matrix_set_zero(m);
	//initialize selection params and allocate stuff
	gammas = gsl_vector_alloc(states);
	piStart = gsl_vector_alloc(states);
	piEnd = gsl_vector_alloc(states);
	break;
      case 'd' :
      	/* open data file, errors? */
	infile = fopen(argv[++args], "r");
	printf("%s is datafile\n",argv[args]);
	if (infile == NULL){
		fprintf(stderr,"Error opening infile! ARRRRR!!!!\n");
		exit(1);
	}
	/* go through infile and get snp info; expect pos, i, n */
	snpNumber = 0;
	maxSampleSize = 0;
	while (fscanf(infile, "%ld\t%d\t%d", &pos, &i, &n) != EOF){
		data[snpNumber].pos = pos;
		data[snpNumber].i = i;
		data[snpNumber].n = n;
		snpNumber++;
		if(n > maxSampleSize){
		  maxSampleSize = n;
		}
	}
	fclose(infile);
	break;
	//either whole transition matrix needs to be defined or gamma and omega
      case 'p':
	//get transition matrix
	infile = fopen(argv[++args], "r");
	if (infile == NULL){
	  fprintf(stderr,"Error opening matrixfile! ARRRRR!!!!\n");
	  exit(1);
	}
	printf("%s is matrixfile\n",argv[args]);
	gsl_matrix_fscanf(infile, m);
	//gsl_matrix_fprintf(stdout,m,"%f");
	fclose(infile);
	break;
      case's':
	//get start vector
		infile = fopen(argv[++args], "r");
	if (infile == NULL){
	  fprintf(stderr,"Error opening start prob. file! ARRRRR!!!!\n");
	  exit(1);
	}
	printf("%s is start prob. file\n",argv[args]);
	gsl_vector_fscanf(infile, piStart);
	fclose(infile);
	smoothFlag = 'piSet';
	break;
      case'e':
	//get end vector
		infile = fopen(argv[++args], "r");
	if (infile == NULL){
	  fprintf(stderr,"Error opening end prob. file! ARRRRR!!!!\n");
	  exit(1);
	}
	printf("%s is end prob. file\n",argv[args]);
	gsl_vector_fscanf(infile, piEnd);
	smoothFlag = 'piSet';
	fclose(infile);
	break;
      case 'E':
        //get emission matrix
        infile = fopen(argv[++args],"r");
        sampleSize = atoi(argv[++args]);
        if (infile == NULL){
          fprintf(stderr,"Error opening end emission prob. file! ARRRRR!!!!\n");
          exit(1);
        }
        printf("%s is emission prob. file; sampleSize=%d\n",argv[args-1],sampleSize);
	emisMatrix = gsl_matrix_alloc(states,sampleSize);
        gsl_matrix_fscanf(infile, emisMatrix);
        printf("finished reading emission prob. matrix\n");
        for (i=0;i<states;i++)
        {
            for (j=0;j<sampleSize;j++)
            {
                printf("%f ",gsl_matrix_get(emisMatrix,i,j));
            }
            printf("\n");
        }
        fixedEmisFlag = 1;
        fclose(infile);
        break;

      case 'g':
	//get gamma for transitions
	priorCons = atof(argv[++args]);
	printf("prior probability of conservation defined\n");
	break;
      case 'w':
	//get omega for transitions
	omega =  atof(argv[++args]);
	smoothFlag = 's';
	break;
      case 'o':
	//get outfile name
	outfileName = argv[++args];
	break;
      case 'V':
	//verbose training
	verboseFlag = 1;
	break;
      }
      args++;
    }
  }
  if (omega == 0 && gsl_matrix_get(m, 0, 0) == 0){
    usage();
    exit(1);
  }
}

/* calculateEmissions- this fills up a matrix of Nstates * L observations where the elements are
   the computed allele frequency probabilities drawn from the stationary distribution of Wright's genic
   selection model. The population is assumed to be at equilibrium and all snps are considered to be
    independent draws from the same distributions */
void calculateEmissions(HMM *hmm, gsl_matrix *emisLogs, gsl_vector *gammas, double theta,\
			int snpNumber, struct snp theData[MAXSNPS]){
  int i, j;
  gsl_matrix *probs;
  gsl_vector *pVector;

  //error check
  assert(emisLogs->size1 == hmm->nstates);
  assert(emisLogs->size2 == snpNumber);
  //alloc pVector
  pVector = gsl_vector_alloc(hmm->nstates);
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  //calculate emissions and fill up matrix
  for(i = 0; i < emisLogs->size1; i++){
    gsl_vector_set(pVector,0,gsl_vector_get(gammas,i));
    gsl_vector_set(pVector,1,theta);
    //make prob matrix (note this is in logs)
    siteProbMatrix(probs,pVector,maxSampleSize, sampleSizes);
    //go through snps, calculate probs
    for(j=0; j < emisLogs->size2; j++){
      gsl_matrix_set(emisLogs, i, j, gsl_matrix_get(probs, data[j].n, data[j].i));
    }
  }
  gsl_vector_free(pVector);
  gsl_matrix_free(probs);
}

void initializeEmissions(gsl_matrix *emisLogs, gsl_matrix *emisMatrix, int snpNumber, struct snp theData[MAXSNPS]){
  int i, j;
  int freq;
  double logProb;

  //calculate emissions and fill up matrix
  for(i = 0; i < emisLogs->size1; i++){
    for(j=0; j < emisLogs->size2; j++){
      //need to get the allele frequency here, and then extract the probs from the corresponding emisMatrix elements
      freq = data[j].i;
      logProb = log(gsl_matrix_get(emisMatrix, i, freq));
      //if (emisLogs->size2 - j <= 10)
      //{
      //    fprintf(stderr,"emission probability for state %d, position %d: %f; log=%f\n",i,j,gsl_matrix_get(emisMatrix, i, freq),logProb);
      //}
      gsl_matrix_set(emisLogs, i, j, logProb);
    }
  }
}

/* emAlg- this is mean to train the emissions (the parameter gamma) using the EM algorithm.
 it returns the trained up gamma value  */
double emAlg(HMM *hmm, gsl_vector *startGammas, int snpNumber, struct snp theData[MAXSNPS], int transitionPowerFlag){
  gsl_matrix *alphas, *betas, *posts, *logPosts, *emisLogs, *tempTrans, *powTrans, **jpMatArray;
  gsl_vector *rowSums, *estimateVector;
  double  lastLik, newLik, delta, dummy, val, pFor, pBack;
  gsl_function likFunc;
   struct my_jointLik_params paramSet;
  int i, j,  k, count;

  //initialize and allocate
  emisLogs = gsl_matrix_alloc(hmm->nstates,snpNumber);
  if(fixedEmisFlag)
  {
      initializeEmissions(emisLogs, emisMatrix, snpNumber, theData);
  }
  else
  {
      calculateEmissions(hmm, emisLogs, startGammas, theta1, snpNumber, theData);
  }
  posts = gsl_matrix_alloc(hmm->nstates,snpNumber);
  logPosts =  gsl_matrix_alloc(hmm->nstates,snpNumber);
  alphas = gsl_matrix_alloc(hmm->nstates,snpNumber);
  betas = gsl_matrix_alloc(hmm->nstates,snpNumber);
  rowSums = gsl_vector_alloc(hmm->nstates);
  delta = 1e-4;
  newLik = 1 ;
  lastLik = 0;
  tempTrans = gsl_matrix_alloc(hmm->nstates,hmm->nstates);

  //make jpMatArray
  jpMatArray = malloc(sizeof(gsl_matrix) * snpNumber);
  for(i = 0; i < snpNumber; i++){
    //each is nstate x nstate transition matrix
    jpMatArray[i] = gsl_matrix_alloc(hmm->nstates,hmm->nstates);
    gsl_matrix_set_all(jpMatArray[i],0.0);
  }

				     //setup parameter estimate vector
  estimateVector = gsl_vector_alloc(hmm->nstates); // 1 mutation parameter, n-1 selection parameters
  gsl_vector_set(estimateVector,0,theta1);
  for(i = 1; i < hmm->nstates; i++){
    gsl_vector_set(estimateVector,i,gsl_vector_get(startGammas,i));
  }

  //initializing joint lik params
  paramSet.data = theData;
  paramSet.snpNumber = snpNumber;
  paramSet.nstates = hmm->nstates;
  paramSet.maxSampleSize = maxSampleSize;
  paramSet.sampleSizeVector = sampleSizes;
  paramSet.sfsBools = sfsCounts;
  paramSet.startingPoint = estimateVector;

  //main loop
  count = 0;
  while ( fabs(newLik - lastLik) >  delta){
    transitionPowerFlag = 0;

    //E step
    //set lastLik to newLik
    dummy = newLik;
    lastLik = dummy;

    //recalculate emissions
    if (!fixedEmisFlag)
    {
        calculateEmissions(hmm, emisLogs, startGammas, theta1, snpNumber, theData);
    }

    //run forwards/backwards and calculate posterior probs
    pFor = forwardAlg(hmm, alphas, emisLogs, snpNumber, transitionPowerFlag);
    pBack = backwardAlg(hmm, betas, emisLogs, snpNumber,transitionPowerFlag);
    //print current status
    printf("%d\tlik= %f\t\ttheta_hat= %f", count,pFor,theta1);
    for(i = 1; i < hmm->nstates; i++){
      printf("\tgamma_hat= %f",gsl_vector_get(startGammas, i));
    }
    printf("\n");
    fflush(stdout);
    //check probs; calculate posterior probs; set "weights"
    assert(fabs(pFor - pBack) < 0.1);
    posteriorProbsReduced(hmm, posts, pFor, alphas, betas, emisLogs, snpNumber,transitionPowerFlag);
    posteriorProbsReducedLog(hmm, logPosts, pFor, alphas, betas, emisLogs, snpNumber,transitionPowerFlag);

    //reestimate the new shiny way
    reestimateTransitions(hmm, posts, jpMatArray,  pFor, alphas, betas, emisLogs, snpNumber,transitionPowerFlag);

    //update start and end vectors based on posteriors
    for(j = 0; j < hmm->nstates; j++){
      gsl_vector_set(hmm->piStart,j,gsl_matrix_get(posts,j,0));
      gsl_vector_set(hmm->piEnd,j,gsl_matrix_get(posts,j,snpNumber - 1));
      gsl_vector_set(hmm->piStart_scores, j, log(gsl_vector_get(hmm->piStart,j)));
      gsl_vector_set(hmm->piEnd_scores, j, log(gsl_vector_get(hmm->piEnd,j)));
    }

    // M step for other params
    if (!fixedEmisFlag)
    {
        paramSet.posts = posts;
        paramSet.startingPoint = estimateVector;
        likFunc.params = &paramSet;
        jointMLEstNState(estimateVector, &newLik, &paramSet);

        //maximization
        theta1 = gsl_vector_get(estimateVector,0);
        for(j = 1; j < hmm->nstates; j++){
          gsl_vector_set(startGammas,j, gsl_vector_get(estimateVector, j));
        }
    }
    newLik = pFor;
    count += 1;
    if (verboseFlag) printHMM(hmm);
  }

 //free jpMatArray
  for(i = 0; i < snpNumber; i++){
    //each is nstate x nstate transition matrix
    gsl_matrix_free(jpMatArray[i]);
  }
  //free(jpMatArray);

  gsl_matrix_free(emisLogs);
  gsl_matrix_free(posts);
  gsl_matrix_free(alphas);
  gsl_matrix_free(betas);
  gsl_vector_free(rowSums);
  gsl_vector_free(estimateVector);
  printHMM(hmm);
  return(gsl_vector_get(startGammas,1));
}

/*This returns an array of arrays of ints which indicate the Viterbi decoded elements from the model.
  The arrays returned are in the format startIndex, stopIndex */
int **elementArray(gsl_vector *vits, double state, int *pCount){
  int i, **array, elementCount, size, last;

  //allocate storage for element arrays
  array = malloc(vits->size * sizeof(int *));
  for(i = 0; i < vits->size; i++){
    array[i] = malloc(3 * sizeof(int));
  }
  //set last to initial obs
  elementCount = 0;
  size = state; //this is the first one we're about to find!
  //find first occurence of state
  for(last = 0; last < vits->size; last++){
    if (gsl_vector_get(vits,last) == state){
      break;
    }
  }

  //iterate through Viterbi's, find runs
  for(i = last+1; i < vits->size; i++){
    if (gsl_vector_get(vits,last) == gsl_vector_get(vits,i)){
      if (size == 0){ //new one?
	last = i;
      }
      size++;
      if ( i == vits->size - 1 && size > 1){
	array[elementCount][0] = last;
	array[elementCount][1] = i;
	elementCount++;
      }
    }
    else{
      if (size > 1){
	array[elementCount][0] = last;
	array[elementCount][1] = i - 1;
	array[elementCount][2] = size;
	elementCount++;
	size = 0;
	last = i - 1 ;
      }
      else{
	size = 0;
      }
    }
  }
  *pCount = elementCount;
  return(array);
}

void scoreElements(HMM *h, gsl_matrix *emisLogs, gsl_matrix *vits, char *outfileName, int transitionPowerFlag){
  gsl_vector *include, *include2;
  FILE *outfile;
  int elementCount, **elements, i, j;
  double score1;

  //open file to print to
  outfile = fopen(outfileName, "w");
  if (outfile == NULL){
    fprintf(stderr,"Error opening outfile! ARRRRR!!!!\n");
    exit(1);
  }
  //print header
  fprintf(outfile,"state\tstart\tstop\tlog_odds_score\tlength\n");

  //alloc state inclusion vectors
  include = gsl_vector_alloc(h->nstates);
  include2 = gsl_vector_alloc(h->nstates);

  //iterate across states
  for(j = 0; j < h->nstates; j++){
    //get state elements from Viterbi
    elementCount = 0;
    elements = elementArray(vits, j, &elementCount);
    //set up state sets for element scoring
    gsl_vector_set_all(include,0);
    gsl_vector_set(include, j, 1);
    gsl_vector_set_all(include2,1);
    gsl_vector_set(include2, j, 0);
    //score elements
    for(i=0; i < elementCount; i++){
      score1 = hmm_subpath_log_odds(h, emisLogs, include,include2, elements[i][0], \
				    elements[i][1] - elements[i][0] + 1 , transitionPowerFlag);
      fprintf(outfile,"%d\t%d\t%d\t%f\t%d\n",j,(int) gsl_vector_get(h->otherData,elements[i][0]), \
	      (int) gsl_vector_get(h->otherData,elements[i][1]),score1, (int) elements[i][2]);
      free(elements[i]);
    }
    free(elements);
  }
  fclose(outfile);
}

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
/**
 *
 * @file pcunmqr_blgtrd.c
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @date 2011-05-15
 * @generated c Tue Jan  7 11:45:13 2014
 *
 **/
#include <lapacke.h>
#include <cblas.h>
#include "common.h"
#undef REAL
#define COMPLEX

#define E(_m, _n) (E   + (_m) + LDE * (_n) )
#define V(_m)     (V   + (_m))
#define T(_m)     (T   + (_m))
#define TAU(_m)   (TAU + (_m))
/***************************************************************************
 *  Parallel apply Q2 from bulgechasing matrices - static scheduling
 *  Lower case is treated
 **/
    /*
     * side == PlasmaLeft:
     *     meaning apply E = Q*E = (q_1*q_2*.....*q_n) * E ==> so
     *     traverse Vs in reverse order (forward) from q_n to q_1 Also
     *     E is splitten by block of col over cores because each apply
     *     consist in a block of row (horizontal block)
     */
    /*
     * side == PlasmaRight:
     *     meaning apply E = E*Q = E * (q_1*q_2*.....*q_n) ==> so
     *     traverse Vs in normal order (forward) from q_1 to q_n Also
     *     E is splitten by block of row over core because each apply
     *     consist in a block of col (vertical block)
     */
/***************************************************************************/
void plasma_pcunmqr_blgtrd(plasma_context_t *plasma)
{
    int my_core_id = PLASMA_RANK;
    int cores_num  = plasma->world_size;

    /*===========================*/
    int N, NB, NE, LDE, Vblksiz, WANTZ;
    PLASMA_enum side;
    PLASMA_enum trans;
    PLASMA_Complex32_t *V;
    PLASMA_Complex32_t *T;
    PLASMA_Complex32_t *TAU;
    PLASMA_Complex32_t *E;
    PLASMA_sequence *sequence;
    PLASMA_request *request;

    /*===========================
     *   local variables
     *===========================*/
    PLASMA_Complex32_t *WORK;
    int LDT, LDV;
    int Vm, Vn, mt, nt;
    int myrow, mycol, blkj, blki;
    int firstrow, nbcolinvolvd;
    int blkid, vpos, taupos, tpos;
    int chunkid, nbchunk, colpercore, corest, corelen, len, col;
    int coreid, allcoresnb, maxrequiredcores;
    int lchunkid, rchunkid, halfchunk, nbportion, sw;
    int LWORK;
    int standalonework = 0 ;
    int versionL, versionR;

    plasma_unpack_args_14(side, trans, N, NB, NE, Vblksiz, WANTZ, V, T, TAU, E, LDE, sequence, request);
    if (sequence->status != PLASMA_SUCCESS)
        return;

    /* Quick return */
    if ( N == 0 ) {
        return;
    }
    if ( NB == 0 ) {
        return;
    }
    if ( NE == 0 ) {
        return;
    }
    /* ==========================================
     * some infor for developer
     * Initialisation and checking nb of cores
     * ==========================================*/
    /* we have to 2 algo for left (113 114) and 2 algo for right (91 92)
     * which correspond to versionL versionR.
     * They ae very similar (detail explained in tech report and matlab code)
     * however version 114 and 92 improve locality.
     * while version 113 is used in case WNATZ=1 (construct Q2) which allow
     * the construction to be done in an optimized way taking into
     * consideration that the matrix is Identity so making less flops.
     *
     */
    versionL = 113;
    versionR = 92;
    LDT      = Vblksiz;
    LDV      = NB+Vblksiz-1;
    /* use colpercore =  N/cores_num; :if i want to split E into
     * cores_num chunk so I will have large chunk for each core.
     * However I prefer to split E into chunk of small size where
     * I guarantee that blas3 occur and the size of chunk remain into
     * cache, I think it is better. than I will have different chunk of
     * small size per core and i will do each chunk till the end and
     * then move to the second one for data locality
     *
     * version v1: for each chunck it apply all the V's then move to
     *                    the other chunck. the locality here inside each
     *                    chunck meaning that thread "t" apply V_k then move
     *                    to V_k+1 which overlap with V_k meaning that the
     *                    E_k+1 overlap with E_k. so here is the
     *                    locality however thread t had to read V_k+1 and
     *                    T_k+1 at each apply. note that all thread if they
     *                    run at same speed they might reading the same V_k
     *                    and T_k at the same time.
     *
     * version v2: for each V_k and T_k thread t apply those V_k and
     *                    T_k to E_k for all its chunck, then move to V_k+1
     *                    T_k+1 and apply them to E_k+1 on all the chunck,
     *                    and so on. the difference is that, thread keep V_k
     *                    and T_K while move over the E_k.
     *                    both version look like similar in perf.
     *
     *  THIS IS v1 CODE
     * */

    if(WANTZ==1)
        colpercore =  plasma_ceildiv(NE,2*cores_num);
    else
        colpercore =  plasma_ceildiv(NE,cores_num);

    if(colpercore>1000)
        colpercore =  plasma_ceildiv(colpercore,10);
    else if(colpercore>800)
        colpercore =  plasma_ceildiv(colpercore,8);
    else if(colpercore>600)
        colpercore =  plasma_ceildiv(colpercore,6);
    else if(colpercore>400)
        colpercore =  plasma_ceildiv(colpercore,4);
    else if(colpercore>200)
        colpercore =  plasma_ceildiv(colpercore,2);
    if(colpercore>200)
        colpercore=120;
    if(colpercore<30)
        colpercore=32;
    /*colpercore = N make the code sequential running on thread=0;*/
    nbchunk          =  plasma_ceildiv(NE, colpercore);
    allcoresnb       = cores_num;
    maxrequiredcores = nbchunk;
    if( maxrequiredcores < 1 )
        maxrequiredcores = 1;
    if(allcoresnb > maxrequiredcores)
        allcoresnb = maxrequiredcores;

#if defined (ENABLE_DEBUG)
    if(my_core_id==0)
        printf("  APPLY Q_v1   parallel with threads %d   nbchunk %d  colpercore %d  N %d   NB %d   Vblksiz %d SIDE %c TRANS %c versionL %d versionR %d WANTZ %d \n",
               cores_num,nbchunk,colpercore,N,NB,Vblksiz,lapack_const(side),lapack_const(trans),versionL,versionR,WANTZ);
#endif
    /* =========================================
     * case NB = 1  special case.
     * just make the element real for complex
     * =========================================
     */
    if ( NB == 1 ) {
        /* NOTE in CASE USED STANDALONE FUNCTION
         * In COMPLEX matrix Z need to be scaled by the TAU
         * generated during the make off-diagonal elements real
         */
        /*
         * In case where the first stage has been done we are sure
         * that all element of E are real so no need to go through
         * NB=1.  However in case where this function need to be used
         * for only band matrices meaning only stage2 has been called
         * then it require to make all off-diagonal elements real and
         * so remove the return from here and from the bulgechasing
         * function
         */
        if(WANTZ==1){
            PLASMA_Complex32_t zone = 1.0;
            memset(E, 0, LDE*N*sizeof(PLASMA_Complex32_t));
            for(sw=0; sw<NE; sw++)
                E[sw+sw*LDE] = zone;
        }
        if(standalonework==0){
            return;
        }
        else{
#ifdef COMPLEX
            for (chunkid = 0; chunkid<nbchunk; chunkid++) {
                coreid  = chunkid%allcoresnb;
                corest  = chunkid*colpercore;
                corelen = min(colpercore, (NE-(chunkid*colpercore)));
                if( my_core_id==coreid ){
                    if( side==PlasmaLeft ){
                        for (mycol =1; mycol<NE; mycol++){
                            cblas_cscal(corelen, TAU(mycol), E(mycol, corest), LDE);
                        }
                    }else{
                        PLASMA_Complex32_t ztmp;
                        /*Upper case need to be checked*/
                        for (mycol = corest; mycol<corest+corelen; mycol++){
                            ztmp = conjf(*TAU(mycol));
                            cblas_cscal(N, &ztmp, E(0,mycol), 1);
                        }
                    }
                }
            }
#endif
            return;
        }
    }
    /* =========================================
     *       case NB >1  main code
     * =========================================
     */
    LWORK = 2*colpercore*max(Vblksiz,64);
    WORK  = (PLASMA_Complex32_t*) plasma_private_alloc(plasma, LWORK, PlasmaComplexFloat);
    /* WANTZ = 1 meaning E is IDENTITY so form Q using optimized update.
     *         So we use the reverse order from small q to large one,
     *         so from q_n to q_1 so Left update to Identity.
     *         Use version 113 because in 114 we need to update the whole matrix and not in icreasing order.
     * WANTZ = 2 meaning E is a full matrix and need to be updated from Left or Right so use normal update
     * */
    if ( WANTZ==1 ) {
        versionL = 113;
        halfchunk = plasma_ceildiv(nbchunk,2);
        for (lchunkid = 0; lchunkid<halfchunk; lchunkid++) {
            rchunkid  = (nbchunk-1) - lchunkid;
            nbportion = lchunkid == rchunkid ? 1 : 2;
            /* sw is the switch between left and right side chunk
             * it is used to have same load distribution of work */

            for (sw = 0; sw<nbportion; sw++) {
                chunkid = sw == 0 ? lchunkid : rchunkid;
                coreid  = lchunkid%allcoresnb;
                corest  = chunkid*colpercore;
                corelen = min(colpercore, (NE-(chunkid*colpercore)));
                if( my_core_id == coreid ) {
                    /*
                     * Version 113:
                     * loop over the block_col (nt) and for each find the
                     * number of tiles (mt) in this block_col. then loop over mt, find
                     * the size of the V's(Vm,Vn) and apply it to the corresponding
                     * portion of E.
                     */
                    nt  = plasma_ceildiv((N-1),Vblksiz);
                    for (blkj=nt-1; blkj>=0; blkj--) {
                        /* the index of the first row on the top of block (blkj) */
                        firstrow = blkj * Vblksiz + 1;
                        /*find the number of tile for this block */
                        if( blkj == nt-1 )
                            mt = plasma_ceildiv( N -  firstrow,    NB);
                        else
                            mt = plasma_ceildiv( N - (firstrow+1), NB);
                        /*loop over the tiles find the size of the Vs and apply it */
                        for (blki=mt; blki>0; blki--) {
                            /*calculate the size of each losange of Vs= (Vm,Vn)*/
                            myrow     = firstrow + (mt-blki)*NB;
                            mycol     = blkj*Vblksiz;
                            Vm = min( NB+Vblksiz-1, N-myrow);
                            if( ( blkj == nt-1 ) && ( blki == mt ) ){
                                Vn = min (Vblksiz, Vm);
                            } else {
                                Vn = min (Vblksiz, Vm-1);
                            }
                            /*calculate the pointer to the Vs and the Ts.
                             * Note that Vs and Ts have special storage done
                             * by the bulgechasing function*/
                            findVTpos(N,NB,Vblksiz,mycol,myrow, &vpos, &taupos, &tpos, &blkid);
                            if((Vm>0)&&(Vn>0)){
                                col = max(mycol,corest);
                                len = corelen - (col-corest);
                                if(side==PlasmaLeft){
                                    if( len > 0 )
                                        CORE_clarfb_gemm(
                                            side, trans,
                                            PlasmaForward, PlasmaColumnwise,
                                            Vm, len, Vn, V(vpos), LDV, T(tpos), LDT, E(myrow,col), LDE,  WORK, len);
                                }
                                else{
                                    if( len > 0 )
                                        CORE_clarfb_gemm(
                                            side, trans,
                                            PlasmaForward, PlasmaColumnwise,
                                            len, Vm, Vn, V(vpos), LDV, T(tpos), LDT, E(col, myrow), LDE,  WORK, len);
                                }
                            }
                        }
                    }
                } /* END my_core_id=coreid */
            } /* END of sw  */
        } /* END loop over the chunk */
    } /* END if WANTZ=1 */
    else{
        /*
         * WANTZ != 1
         */
        for (chunkid = 0; chunkid<nbchunk; chunkid++) {
            coreid  = chunkid%allcoresnb;
            corest  = chunkid*colpercore;
            corelen = min(colpercore, (NE-(chunkid*colpercore)));
            if(corelen < 0)
                corelen = 0;
            if( my_core_id == coreid ) {
                /*
                 * PlasmaLeft
                 */
                //if( side == PlasmaLeft ) {
                if( versionL == 113 ) {
                    /*
                     * Version 113:
                     * loop over the block_col (nt) and for each find the
                     * number of tiles (mt) in this block_col. then loop over mt, find
                     * the size of the V's(Vm,Vn) and apply it to the corresponding
                     * portion of E.
                     */
                    if( versionL == 113 ) {
                        nt  = plasma_ceildiv((N-1),Vblksiz);
                        for (blkj=nt-1; blkj>=0; blkj--) {
                            /* the index of the first row on the top of block (blkj) */
                            firstrow = blkj * Vblksiz + 1;
                            /*find the number of tile for this block */
                            if( blkj == nt-1 )
                                mt = plasma_ceildiv( N -  firstrow,    NB);
                            else
                                mt = plasma_ceildiv( N - (firstrow+1), NB);
                            /*loop over the tiles find the size of the Vs and apply it */
                            for (blki=mt; blki>0; blki--) {
                                /*calculate the size of each losange of Vs= (Vm,Vn)*/
                                myrow     = firstrow + (mt-blki)*NB;
                                mycol     = blkj*Vblksiz;
                                Vm = min( NB+Vblksiz-1, N-myrow);
                                if( ( blkj == nt-1 ) && ( blki == mt ) ){
                                    Vn = min (Vblksiz, Vm);
                                } else {
                                    Vn = min (Vblksiz, Vm-1);
                                }
                                /*calculate the pointer to the Vs and the Ts.
                                 * Note that Vs and Ts have special storage done
                                 * by the bulgechasing function*/
                                findVTpos(N,NB,Vblksiz,mycol,myrow, &vpos, &taupos, &tpos, &blkid);
                                if((Vm>0)&&(Vn>0)){
                                    if( side == PlasmaLeft ){
                                        CORE_clarfb_gemm(
                                            PlasmaLeft, trans,
                                            PlasmaForward, PlasmaColumnwise,
                                            Vm, corelen, Vn, V(vpos), LDV, T(tpos), LDT, E(myrow,corest), LDE,  WORK, corelen);

                                    }else{
                                        CORE_clarfb_gemm(
                                            PlasmaRight, trans,
                                            PlasmaForward, PlasmaColumnwise,
                                            corelen, Vm, Vn, V(vpos), LDV, T(tpos), LDT, E(corest, myrow), LDE,  WORK, corelen);
                                    }
                                }
                            }
                        }
                    }
                    /*
                     * Version 114:
                     * loop over the block_row (mt) and for each find diagonally the
                     * number of tiles (nt) in this block_row. then loop over nt, find
                     * the size of the V's(Vm,Vn) and apply it to the corresponding
                     * portion of E.
                     */
                    else {
                        mt    = plasma_ceildiv((N-1),NB);
                        for (blki = mt; blki>0; blki--) {
                            /* nbcolinvolvd = number of column corresponding to this block_row (blki) */
                            nbcolinvolvd = min(N-1, blki*NB);
                            /*find the number of tile for this block (diagonal row of tiles) */
                            nt = plasma_ceildiv(nbcolinvolvd,Vblksiz);
                            /*loop over the tiles find the size of the Vs and apply it */
                            for (blkj = nt-1; blkj>=0; blkj--) {
                                /* the index of the first row of the first col meaning
                                 * the block on the top left (blki) */
                                firstrow = (mt-blki)*NB+1;
                                /*calculate the size of each losange of Vs= (Vm,Vn)*/
                                myrow    = firstrow + blkj*Vblksiz;
                                mycol    = blkj*Vblksiz;
                                Vm = min( NB+Vblksiz-1, N-myrow);
                                if( ( blkj == nt-1 ) && ( blki == mt ) ){
                                    Vn = min (Vblksiz, Vm);
                                }else{
                                    Vn = min (Vblksiz, Vm-1);
                                }
                                /*calculate the pointer to the Vs and the Ts.
                                 * Note that Vs and Ts have special storage done
                                 * by the bulgechasing function*/
                                findVTpos(N,NB,Vblksiz,mycol,myrow, &vpos, &taupos, &tpos, &blkid);
                                if((Vm>0)&&(Vn>0)) {
                                    CORE_clarfb_gemm(
                                        PlasmaLeft, trans,
                                        PlasmaForward, PlasmaColumnwise,
                                        Vm, corelen, Vn, V(vpos), LDV, T(tpos), LDT, E(myrow,corest), LDE,  WORK, corelen);
                                }
                            }
                        }
                    }
                }
                /*
                 * PlasmaRight
                 */
                else {
                    /*
                     * Version 91:
                     */
                    if( versionR == 91 ) {
                        nt  = plasma_ceildiv((N-1),Vblksiz);
                        for (blkj=0; blkj<nt; blkj++) {
                            /* the index of the first myrow on the top of block (blkj) */
                            firstrow = blkj * Vblksiz + 1;
                            /*find the number of tile for this block */
                            if( blkj == nt-1 )
                                mt = plasma_ceildiv( N -  firstrow,    NB);
                            else
                                mt = plasma_ceildiv( N - (firstrow+1), NB);
                            /*loop over the tiles find the size of the Vs and apply it */
                            for (blki=1; blki<=mt; blki++) {
                                /*calculate the size of each losange of Vs= (Vm,Vn)*/
                                myrow  = firstrow + (mt-blki)*NB;
                                Vm = min( NB+Vblksiz-1, N-myrow);
                                if( ( blkj == nt-1 ) && ( blki == mt ) )
                                {
                                    Vn = min (Vblksiz, Vm);
                                }else{
                                    Vn = min (Vblksiz, Vm-1);
                                }
                                /*calculate the pointer to the Vs and the Ts.
                                 * Note that Vs and Ts have special storage done
                                 * by the bulgechasing function*/
                                mycol     = blkj*Vblksiz;
                                findVTpos(N,NB,Vblksiz,mycol,myrow, &vpos, &taupos, &tpos, &blkid);
                                if((Vm>0)&&(Vn>0)){
                                    CORE_clarfb_gemm(
                                        PlasmaRight,
                                        trans,
                                        PlasmaForward,
                                        PlasmaColumnwise,
                                        corelen, Vm, Vn, V(vpos), LDV, T(tpos), LDT, E(corest,myrow), LDE,  WORK, corelen);
                                }
                            }
                        }
                    }
                    /*trans
                     * Version 92:
                     */
                    else {
                        mt    = plasma_ceildiv((N-1),NB);
                        for (blki = 1; blki<=mt; blki++) {
                            /* nbcolinvolvd = number of column corresponding to this block_row (blki) */
                            nbcolinvolvd = min(N-1, blki*NB);
                            /*find the number of tile for this block (diagonal row of tiles) */
                            nt = plasma_ceildiv(nbcolinvolvd,Vblksiz);
                            /*loop over the tiles find the size of the Vs and apply it */
                            for (blkj = 0; blkj<nt; blkj++) {
                                /* the index of the first row of the first col meaning
                                 * the block on the top left (blki) */
                                firstrow = (mt-blki)*NB+1;
                                /*calculate the size of each losange of Vs= (Vm,Vn)*/
                                myrow    = firstrow + blkj*Vblksiz;
                                mycol    = blkj*Vblksiz;
                                Vm = min( NB+Vblksiz-1, N-myrow);
                                if( ( blkj == nt-1 ) && ( blki == mt ) ){
                                    Vn = min (Vblksiz, Vm);
                                }else{
                                    Vn = min (Vblksiz, Vm-1);
                                }
                                /*calculate the pointer to the Vs and the Ts.
                                 * Note that Vs and Ts have special storage done
                                 * by the bulgechasing function*/
                                findVTpos(N,NB,Vblksiz,mycol,myrow, &vpos, &taupos, &tpos, &blkid);
                                if((Vm>0)&&(Vn>0)) {
                                    CORE_clarfb_gemm(
                                        PlasmaRight, trans,
                                        PlasmaForward, PlasmaColumnwise,
                                        corelen, Vm, Vn, V(vpos), LDV, T(tpos), LDT, E(corest,myrow), LDE,  WORK, corelen);
                                }
                            }
                        }
                    }
                }
            } /* END my_core_id=coreid */
        } /* END loop over the chunk */
    } /* END ELSE of WANTZ == 1 */
    plasma_private_free(plasma, WORK);
}
#undef E
#undef V
#undef T
#undef TAU
/***************************************************************************/
#undef COMPLEX

ext: c


text:
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <math.h>
#include <options/options.h>

#include "../csm/csm_all.h"

struct ld_exp_tro1_params {
	int seed;

	double max_xy_error;
	double max_theta_error_deg;

	const char* file_input;
	const char* file_output1;
	const char* file_output2;

	int num_per_scan;

	int debug;
};

const char * banner =
	"This program prepares the data for one of the experiments. \n\n"
		"The input is any sensor log (Carmen or JSON format) \n"
		"The output are two files that contain laser_ref and laser_sens\n"
		"(you have to match the i-th scan in the first file with the i-th\n"
		" in the second).\n\n"
		"The two files contain exactly the same date but for the 'odometry' field\n"
		"The odometry error is uniform in the intervals given.\n";

int main(int argc, const char ** argv) {
	sm_set_program_name(argv[0]);

	struct ld_exp_tro1_params p;

	options_banner(banner);

	struct option* ops = options_allocate(10);
	options_double(ops, "max_xy_error", &p.max_xy_error, 10.0, "Maximum error for x,y (m)");
	options_double(ops, "max_theta_error_deg", &p.max_theta_error_deg, 10.0, "Maximum error for orientation (deg)");
	options_int   (ops, "seed", &p.seed, 0, "Seed for random number generator (if 0, use GSL_RNG_SEED env. variable).");

	options_int(ops, "num_per_scan", &p.num_per_scan, 10, "Number of trials for each scan.");

	options_string(ops, "in", &p.file_input, "stdin", "Input file ");
	options_string(ops, "out1", &p.file_output1, "stdout", "Output file for first scan");
	options_string(ops, "out2", &p.file_output2, "stdout", "Output file for second scan");

	options_int(ops, "debug", &p.debug, 0, "Shows debug information");

	if(!options_parse_args(ops, argc, argv)) {
		options_print_help(ops, stderr);
		return -1;
	}

	sm_debug_write(p.debug);

	gsl_rng_env_setup();
	gsl_rng * rng = gsl_rng_alloc (gsl_rng_ranlxs0);
	if(p.seed != 0)
	gsl_rng_set(rng, (unsigned int) p.seed);

	/* Open the two output files (possibly the same one) */

	FILE * in = open_file_for_reading(p.file_input);
	if(!in) return -3;

	FILE * out1 = open_file_for_writing(p.file_output1);
	if(!out1) return -2;

	FILE * out2;
	if(!strcmp(p.file_output1, p.file_output2)) {
		out1 = out2;
	} else {
		out2 = open_file_for_writing(p.file_output2);
		if(!out2) return -2;
	}

	/* Read laser data from input file */
	LDP ld; int count=0;
	while( (ld = ld_read_smart(in))) {
		count++;
		if(!ld_valid_fields(ld))  {
			sm_error("Invalid laser data (#%d in file)\n", count);
			continue;
		}

		for(int n=0; n < p.num_per_scan; n++) {
			ld->true_pose[0] = 0;
			ld->true_pose[1] = 0;
			ld->true_pose[2] = 0;

			ld->odometry[0] = 0;
			ld->odometry[1] = 0;
			ld->odometry[2] = 0;

			ld_write_as_json(ld, out1);

			ld->odometry[0] = 2*(gsl_rng_uniform(rng)-0.5) * p.max_xy_error;
			ld->odometry[1] = 2*(gsl_rng_uniform(rng)-0.5) * p.max_xy_error;
			ld->odometry[2] = 2*(gsl_rng_uniform(rng)-0.5) * deg2rad(p.max_theta_error_deg);

			ld_write_as_json(ld, out2);
		}

		ld_free(ld);
	}

	return 0;
}

max_stars_repo_name: ozaslan/basics

max_issues_repo_name: ozaslan/basics

max_forks_repo_name: ozaslan/basics


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

#include <gsl/gsl_ieee_utils.h>

#include "test_longley.c"
#include "test_filip.c"
#include "test_pontius.c"
#include "test_brown.c"
#include "test_enso.c"
#include "test_kirby2.c"
#include "test_hahn1.c"
#include "test_nelson.c"
#include "test_fn.c"
#include "test_estimator.c"

void
test_lmder (gsl_multifit_function_fdf * f, double x0[],
            double * X, double F[], double * cov);

void
test_fdf (const char * name, gsl_multifit_function_fdf * f,
          double x0[], double x[], double sumsq,
          double sigma[]);

int
main (void)
{
  gsl_ieee_env_setup();

  test_longley();
  test_filip();
  test_pontius();
  test_estimator();

  {
    gsl_multifit_function_fdf f = make_fdf (&brown_f, &brown_df, &brown_fdf,
                                            brown_N, brown_P, 0);

    test_lmder(&f, brown_x0, &brown_X[0][0], brown_F, &brown_cov[0][0]);
  }

  {
    gsl_multifit_function_fdf f = make_fdf (&enso_f, &enso_df, &enso_fdf,
                                            enso_N, enso_P, 0);

    test_fdf("nist-ENSO", &f, enso_x0, enso_x, enso_sumsq, enso_sigma);
  }

  {
    gsl_multifit_function_fdf f = make_fdf (&kirby2_f, &kirby2_df, &kirby2_fdf,
                                            kirby2_N, kirby2_P, 0);

    test_fdf("nist-kirby2", &f, kirby2_x0, kirby2_x, kirby2_sumsq, kirby2_sigma);
  }

  {
    gsl_multifit_function_fdf f = make_fdf (&hahn1_f, &hahn1_df, &hahn1_fdf,
                                            hahn1_N, hahn1_P, 0);

    test_fdf("nist-hahn1", &f, hahn1_x0, hahn1_x, hahn1_sumsq, hahn1_sigma);
  }

#ifdef JUNK
  {
    gsl_multifit_function_fdf f = make_fdf (&nelson_f, &nelson_df, &nelson_fdf,
                                            nelson_N, nelson_P, 0);

    test_fdf("nist-nelson", &f, nelson_x0, nelson_x, nelson_sumsq, nelson_sigma);
  }
#endif

  /* now summarize the results */

  exit (gsl_test_summary ());
}

void
test_lmder (gsl_multifit_function_fdf * f, double x0[],
            double * X, double F[], double * cov)
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  const size_t n = f->n;
  const size_t p = f->p;

  int status;
  size_t iter = 0, i;

  gsl_vector_view x = gsl_vector_view_array (x0, p);

  T = gsl_multifit_fdfsolver_lmsder;
  s = gsl_multifit_fdfsolver_alloc (T, n, p);
  gsl_multifit_fdfsolver_set (s, f, &x.vector);

  do
    {
      status = gsl_multifit_fdfsolver_iterate (s);

      for (i = 0 ; i < p; i++)
        {
          gsl_test_rel (gsl_vector_get (s->x, i), X[p*iter+i], 1e-5,
                        "lmsder, iter=%u, x%u", iter, i);
        }

      gsl_test_rel (gsl_blas_dnrm2 (s->f), F[iter], 1e-5,
                    "lmsder, iter=%u, f", iter);

      iter++;
    }
  while (iter < 20);

  {
    size_t i, j;
    gsl_matrix * covar = gsl_matrix_alloc (4, 4);
    gsl_multifit_covar (s->J, 0.0, covar);

    for (i = 0; i < 4; i++)
      {
        for (j = 0; j < 4; j++)
          {
            gsl_test_rel (gsl_matrix_get(covar,i,j), cov[i*p + j], 1e-7,
                          "gsl_multifit_covar cov(%d,%d)", i, j) ;
          }
      }

    gsl_matrix_free (covar);
  }

  gsl_multifit_fdfsolver_free (s);

}

void
test_fdf (const char * name, gsl_multifit_function_fdf * f,
          double x0[], double x_final[],
          double f_sumsq, double sigma[])
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  const size_t n = f->n;
  const size_t p = f->p;

  int status;
  size_t iter = 0;

  gsl_vector_view x = gsl_vector_view_array (x0, p);

  T = gsl_multifit_fdfsolver_lmsder;
  s = gsl_multifit_fdfsolver_alloc (T, n, p);
  gsl_multifit_fdfsolver_set (s, f, &x.vector);

  do
    {
      status = gsl_multifit_fdfsolver_iterate (s);

#ifdef DEBUG
       printf("iter = %d  status = %d  |f| = %.18e x = \n",
         iter, status, gsl_blas_dnrm2 (s->f));

         gsl_vector_fprintf(stdout, s->x, "%.8e");
#endif
      status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-7);

      iter++;
    }
  while (status == GSL_CONTINUE && iter < 1000);

  {
    size_t i;
    gsl_matrix * covar = gsl_matrix_alloc (p, p);
    gsl_multifit_covar (s->J, 0.0, covar);

    for (i = 0 ; i < p; i++)
      {
        gsl_test_rel (gsl_vector_get (s->x, i), x_final[i], 1e-5,
                      "%s, lmsder, x%u", name, i);
      }

    {
      double s2 = pow(gsl_blas_dnrm2 (s->f), 2.0);

      gsl_test_rel (s2, f_sumsq, 1e-5, "%s, lmsder, |f|^2", name);

      for (i = 0; i < p; i++)
        {
          double ei = sqrt(s2/(n-p))*sqrt(gsl_matrix_get(covar,i,i));
          gsl_test_rel (ei, sigma[i], 1e-4,
                        "%s, sigma(%d)", name, i) ;
        }
    }

    gsl_matrix_free (covar);
  }

  gsl_multifit_fdfsolver_free (s);
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
  NAME:
     minmax
  PURPOSE:
     returns the (finite) minimum and the maximum of a vector which is the row/column of a matrix
  CALLING SEQUENCE:
     minmax(gsl_matrix * q, int row, bool isrow, double * min, double * max, bool partial, int partial_indx[3])
  INPUT:
     q            - matrix which holds the vector
     row          - row/column which holds the vector
     isrow        - is it the row or is it the column
  OUTPUT:
     min   - finite minimum
     max   - finite maximum
  REVISION HISTORY:
     2008-09-21 - Written Bovy
*/
#include <gsl/gsl_matrix.h>
#include <float.h>
#include <stdbool.h>
#include <proj_gauss_mixtures.h>

void minmax(gsl_matrix * q, int row, bool isrow, double * min,
	    double * max){
  *max = -DBL_MAX;
  *min = DBL_MAX;
  int dd;
  double temp;
  if (isrow) {
    for (dd = 0; dd != q->size2; ++dd){
	temp  = gsl_matrix_get(q,row,dd);
	if (temp > *max && bovy_isfin(temp))
	  *max = temp;
	if (temp < *min && bovy_isfin(temp))
	  *min = temp;
    }
  }
  else {
    for (dd = 0; dd != q->size1; ++dd){
	temp  = gsl_matrix_get(q,dd,row);
	if (temp > *max && bovy_isfin(temp))
	  *max = temp;
	if (temp < *min && bovy_isfin(temp))
	  *min = temp;
    }
  }

  return ;
}


text:
//Roots of polynomial using compan matrix, as in Matlab/Octave.
//This does roots for each vector X.
//The roots are from an eig using LAPACKE,
//and are tested to match Octave, except sometimes for sorting order.

//For the complex case with non-contiguous vecs in X and Y,
//the answers are correct except for a mysterious sign change.

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

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

int poly2roots_s (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim);
int poly2roots_d (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim);
int poly2roots_c (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim);
int poly2roots_z (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim);

int poly2roots_s (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in poly2roots_s: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    if (N==0u || Lx<2u) {}
    else if (Lx==2u)
    {
        if (Lx==N) { *Y = -*(X+1); *++Y = 0.0f; }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v, ++X, ++Y) { ++X; *Y = -*X; *++Y = 0.0f; }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B)
                {
                    for (size_t b=B; b>0u; --b, ++X, ++Y) { *Y = -*(X+K); *++Y = 0.0f; }
                }
            }
        }
    }
    else
    {
        const size_t Ly = Lx - 1u;
        float x0;
        const int job = 'E', compz = 'N';  //eigenvalues only
        const lapack_int ldh = (int)Ly, n = (int)Ly, ldz = 1;
        const lapack_int ilo = 1, ihi = n;  //avoids balancing
        lapack_int info;
        float *compan, *wr, *wi, zz[1];
        if (!(wr=(float *)malloc((size_t)(n)*sizeof(float)))) { fprintf(stderr,"error in poly2roots_s: problem with malloc. "); perror("malloc"); return 1; }
        if (!(wi=(float *)malloc((size_t)(n)*sizeof(float)))) { fprintf(stderr,"error in poly2roots_s: problem with malloc. "); perror("malloc"); return 1; }
        if (!(compan=(float *)calloc((size_t)(n*n),sizeof(float)))) { fprintf(stderr,"error in poly2roots_s: problem with calloc. "); perror("calloc"); return 1; }

        if (Lx==N)
        {
            ++compan;
            for (size_t l=Ly; l>0u; --l, compan+=Lx) { *compan = 1.0f; }
            compan -= Lx*Ly + 1u;
            x0 = -*X++;
            for (size_t l=Ly; l>0u; --l, ++X, compan+=Ly) { *compan = *X / x0; }
            compan -= Ly*Ly;
            info = LAPACKE_shseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2roots_s: lapacke decomposition failed\n"); return 1; }
            for (size_t l=Ly; l>0u; --l) { *Y++ = *wr++; *Y++ = *wi++; }
            wr -= Ly; wi -= Ly;
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v)
                {
                    for (size_t l=Ly*Ly; l>0u; --l, ++compan) { *compan = 0.0f; }
                    compan -= Lx;
                    for (size_t l=Ly; l>2u; --l, compan-=Lx) { *compan = 1.0f; }
                    *compan-- = 1.0f;
                    x0 = -*X++;
                    for (size_t l=Ly; l>0u; --l, ++X, compan+=Ly) { *compan = *X / x0; }
                    compan -= Ly*Ly;
                    info = LAPACKE_shseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
                    if (info) { fprintf(stderr,"error in poly2roots_s: lapacke decomposition failed\n"); return 1; }
                    for (size_t l=Ly; l>0u; --l) { *Y++ = *wr++; *Y++ = *wi++; }
                    wr -= Ly; wi -= Ly;
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=2u*B*(Ly-1u))
                {
                    for (size_t b=B; b>0u; --b, X-=K*Lx-1u, Y-=2u*K*Ly-2u)
                    {
                        for (size_t l=Ly*Ly; l>0u; --l, ++compan) { *compan = 0.0f; }
                        compan -= Lx;
                        for (size_t l=Ly; l>2u; --l, compan-=Lx) { *compan = 1.0f; }
                        *compan-- = 1.0f;
                        x0 = -*X; X += K;
                        for (size_t l=Ly; l>0u; --l, X+=K, compan+=Ly) { *compan = *X / x0; }
                        compan -= Ly*Ly;
                        info = LAPACKE_shseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
                        if (info) { fprintf(stderr,"error in poly2roots_s: lapacke decomposition failed\n"); return 1; }
                        for (size_t l=Ly; l>0u; --l, Y+=2u*K-1u) { *Y = *wr++; *++Y = *wi++; }
                        wr -= Ly; wi -= Ly;
                    }
                }
            }
        }

        free(compan); free(wr); free(wi);
    }

    return 0;
}

int poly2roots_d (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in poly2roots_d: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    if (N==0u || Lx<2u) {}
    else if (Lx==2u)
    {
        if (Lx==N) { *Y = -*(X+1); *++Y = 0.0; }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v, ++X, ++Y) { ++X; *Y = -*X; *++Y = 0.0; }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B)
                {
                    for (size_t b=B; b>0u; --b, ++X, ++Y) { *Y = -*(X+K); *++Y = 0.0; }
                }
            }
        }
    }
    else
    {
        const size_t Ly = Lx - 1u;
        double x0;

        const int job = 'E', compz = 'N';  //eigenvalues only
        const lapack_int ldh = (int)Ly, n = (int)Ly, ldz = 1;
        const lapack_int ilo = 1, ihi = n;  //avoids balancing
        lapack_int info;
        double *compan, *wr, *wi, zz[1];
        if (!(wr=(double *)malloc((size_t)(n)*sizeof(double)))) { fprintf(stderr,"error in poly2roots_d: problem with malloc. "); perror("malloc"); return 1; }
        if (!(wi=(double *)malloc((size_t)(n)*sizeof(double)))) { fprintf(stderr,"error in poly2roots_d: problem with malloc. "); perror("malloc"); return 1; }
        if (!(compan=(double *)calloc((size_t)(n*n),sizeof(double)))) { fprintf(stderr,"error in poly2roots_d: problem with calloc. "); perror("calloc"); return 1; }

        if (Lx==N)
        {
            ++compan;
            for (size_t l=Ly; l>0u; --l, compan+=Lx) { *compan = 1.0; }
            compan -= Lx*Ly + 1u;
            x0 = -*X++;
            for (size_t l=Ly; l>0u; --l, ++X, compan+=Ly) { *compan = *X / x0; }
            compan -= Ly*Ly;
            info = LAPACKE_dhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2roots_d: lapacke decomposition failed\n"); return 1; }
            for (size_t l=Ly; l>0u; --l) { *Y++ = *wr++; *Y++ = *wi++; }
            wr -= Ly; wi -= Ly;
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v)
                {
                    for (size_t l=Ly*Ly; l>0u; --l, ++compan) { *compan = 0.0; }
                    compan -= Lx;
                    for (size_t l=Ly; l>2u; --l, compan-=Lx) { *compan = 1.0; }
                    *compan-- = 1.0;
                    x0 = -*X++;
                    for (size_t l=Ly; l>0u; --l, ++X, compan+=Ly) { *compan = *X / x0; }
                    compan -= Ly*Ly;
                    info = LAPACKE_dhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
                    if (info) { fprintf(stderr,"error in poly2roots_d: lapacke decomposition failed\n"); return 1; }
                    for (size_t l=Ly; l>0u; --l) { *Y++ = *wr++; *Y++ = *wi++; }
                    wr -= Ly; wi -= Ly;
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=2u*B*(Ly-1u))
                {
                    for (size_t b=B; b>0u; --b, X-=K*Lx-1u, Y-=2u*K*Ly-2u)
                    {
                        for (size_t l=Ly*Ly; l>0u; --l, ++compan) { *compan = 0.0; }
                        compan -= Lx;
                        for (size_t l=Ly; l>2u; --l, compan-=Lx) { *compan = 1.0; }
                        *compan-- = 1.0;
                        x0 = -*X; X += K;
                        for (size_t l=Ly; l>0u; --l, X+=K, compan+=Ly) { *compan = *X / x0; }
                        compan -= Ly*Ly;
                        info = LAPACKE_dhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
                        if (info) { fprintf(stderr,"error in poly2roots_d: lapacke decomposition failed\n"); return 1; }
                        for (size_t l=Ly; l>0u; --l, Y+=2u*K-1u) { *Y = *wr++; *++Y = *wi++; }
                        wr -= Ly; wi -= Ly;
                    }
                }
            }
        }

        free(compan); free(wr); free(wi);
    }

    return 0;
}

int poly2roots_c (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in poly2roots_c: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    if (N==0u || Lx<2u) {}
    else
    {
        const size_t Ly = Lx - 1u;
        const int job = 'E', compz = 'N';  //eigenvalues only
        const lapack_int ldh = (int)Ly, n = (int)Ly, ldz = 1;
        const lapack_int ilo = 1, ihi = n;  //avoids balancing
        lapack_int info;
        float *compan, zz[2], scr, sci;
        if (!(compan=(float *)calloc((size_t)(4*n*n),sizeof(float)))) { fprintf(stderr,"error in poly2roots_c: problem with calloc. "); perror("calloc"); return 1; }

        if (Lx==N)
        {
            compan += 2;
            for (size_t l=Ly; l>0u; --l, compan+=2u*Lx) { *compan = 1.0f; }
            compan -= 2u*Lx*Ly + 2u;
            scr = 1.0f/(*X**X+*(X+1)**(X+1)); sci = *(X+1)*scr; scr *= -*X;
            X += 2;
            for (size_t l=Ly; l>0u; --l, X+=2, compan+=2u*Ly-1u)
            {
                *compan = *X*scr - *(X+1)*sci;
                *++compan = *X*sci + *(X+1)*scr;
            }
            compan -= 2u*Ly*Ly;
            info = LAPACKE_chseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_float *)compan,ldh,(lapack_complex_float *)Y,(lapack_complex_float *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2roots_c: lapacke decomposition failed\n"); return 1; }
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v, Y+=2u*Ly)
                {
                    if (Ly>1u)
                    {
                        for (size_t l=0u; l<2u*Ly*Ly; ++l, ++compan) { *compan = 0.0f; }
                        compan -= 2u*Lx;
                        for (size_t l=Ly; l>2u; --l, compan-=2u*Lx) { *compan = 1.0f; }
                        *compan = 1.0f; compan -= 2;
                    }
                    scr = 1.0f/(*X**X+*(X+1)**(X+1)); sci = *(X+1)*scr; scr *= -*X;
                    X += 2;
                    for (size_t l=Ly; l>0u; --l, X+=2, compan+=2u*Ly-1u)
                    {
                        *compan = *X*scr - *(X+1)*sci;
                        *++compan = *X*sci + *(X+1)*scr;
                    }
                    compan -= 2u*Ly*Ly;
                    info = LAPACKE_chseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_float *)compan,ldh,(lapack_complex_float *)Y,(lapack_complex_float *)zz,ldz);
                    if (info) { fprintf(stderr,"error in poly2roots_c: lapacke decomposition failed\n"); return 1; }
                }
            }
            else
            {
                float *roots;
                if (!(roots=(float *)malloc((size_t)(2*n)*sizeof(float)))) { fprintf(stderr,"error in poly2roots_c: problem with malloc. "); perror("malloc"); return 1; }

                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(Ly-1u))
                {
                    for (size_t b=B; b>0u; --b, X-=2u*K*Lx-2u, Y-=2u*K*Ly-2u)
                    {
                        if (Ly>1u)
                        {
                            for (size_t l=0u; l<2u*Ly*Ly; ++l, ++compan) { *compan = 0.0f; }
                            compan -= 2u*Lx;
                            for (size_t l=Ly; l>2u; --l, compan-=2u*Lx) { *compan = 1.0f; }
                            *compan = 1.0f; compan -= 2;
                        }
                        scr = 1.0f/(*X**X+*(X+1)**(X+1)); sci = *(X+1)*scr; scr *= -*X;
                        X += 2u*K;
                        for (size_t l=Ly; l>0u; --l, X+=2u*K, compan+=2u*Ly-1u)
                        {
                            *compan = *X*scr - *(X+1)*sci;
                            *++compan = *X*sci + *(X+1)*scr;
                        }
                        compan -= 2u*Ly*Ly;
                        info = LAPACKE_chseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_float *)compan,ldh,(lapack_complex_float *)roots,(lapack_complex_float *)zz,ldz);
                        if (info) { fprintf(stderr,"error in poly2roots_c: lapacke decomposition failed\n"); return 1; }
                        for (size_t l=Ly; l>0u; --l, ++roots, Y+=2u*K-1u) { *Y++ = *roots++; *Y = *roots; }
                        roots -= 2u*Ly;
                    }
                }
                free(roots);
            }
        }
        free(compan);
    }

    return 0;
}

int poly2roots_z (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in poly2roots_z: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    if (N==0u || Lx<2u) {}
    else
    {
        const size_t Ly = Lx - 1u;
        const int job = 'E', compz = 'N';  //eigenvalues only
        const lapack_int ldh = (int)Ly, n = (int)Ly, ldz = 1;
        const lapack_int ilo = 1, ihi = n;  //avoids balancing
        lapack_int info;
        double *compan, zz[2], scr, sci;
        if (!(compan=(double *)calloc((size_t)(4*n*n),sizeof(double)))) { fprintf(stderr,"error in poly2roots_z: problem with calloc. "); perror("calloc"); return 1; }

        if (Lx==N)
        {
            compan += 2;
            for (size_t l=Ly; l>0u; --l, compan+=2u*Lx) { *compan = 1.0; }
            compan -= 2u*Lx*Ly + 2u;
            scr = 1.0/(*X**X+*(X+1)**(X+1)); sci = *(X+1)*scr; scr *= -*X;
            X += 2;
            for (size_t l=Ly; l>0u; --l, X+=2, compan+=2u*Ly-1u)
            {
                *compan = *X*scr - *(X+1)*sci;
                *++compan = *X*sci + *(X+1)*scr;
            }
            compan -= 2u*Ly*Ly;
            info = LAPACKE_zhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_double *)compan,ldh,(lapack_complex_double *)Y,(lapack_complex_double *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2roots_z: lapacke decomposition failed\n"); return 1; }
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v, Y+=2u*Ly)
                {
                    if (Ly>1u)
                    {
                        for (size_t l=0u; l<2u*Ly*Ly; ++l, ++compan) { *compan = 0.0; }
                        compan -= 2u*Lx;
                        for (size_t l=Ly; l>2u; --l, compan-=2u*Lx) { *compan = 1.0; }
                        *compan = 1.0; compan -= 2;
                    }
                    scr = 1.0/(*X**X+*(X+1)**(X+1)); sci = *(X+1)*scr; scr *= -*X;
                    X += 2;
                    for (size_t l=Ly; l>0u; --l, X+=2, compan+=2u*Ly-1u)
                    {
                        *compan = *X*scr - *(X+1)*sci;
                        *++compan = *X*sci + *(X+1)*scr;
                    }
                    compan -= 2u*Ly*Ly;
                    info = LAPACKE_zhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_double *)compan,ldh,(lapack_complex_double *)Y,(lapack_complex_double *)zz,ldz);
                    if (info) { fprintf(stderr,"error in poly2roots_z: lapacke decomposition failed\n"); return 1; }
                }
            }
            else
            {
                double *roots;
                if (!(roots=(double *)malloc((size_t)(2*n)*sizeof(double)))) { fprintf(stderr,"error in poly2roots_z: problem with malloc. "); perror("malloc"); return 1; }

                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(Ly-1u))
                {
                    for (size_t b=B; b>0u; --b, X-=2u*K*Lx-2u, Y-=2u*K*Ly-2u)
                    {
                        if (Ly>1u)
                        {
                            for (size_t l=0u; l<2u*Ly*Ly; ++l, ++compan) { *compan = 0.0; }
                            compan -= 2u*Lx;
                            for (size_t l=Ly; l>2u; --l, compan-=2u*Lx) { *compan = 1.0; }
                            *compan = 1.0; compan -= 2;
                        }
                        scr = 1.0/(*X**X+*(X+1)**(X+1)); sci = *(X+1)*scr; scr *= -*X;
                        X += 2u*K;
                        for (size_t l=Ly; l>0u; --l, X+=2u*K, compan+=2u*Ly-1u)
                        {
                            *compan = *X*scr - *(X+1)*sci;
                            *++compan = *X*sci + *(X+1)*scr;
                        }
                        compan -= 2u*Ly*Ly;
                        info = LAPACKE_zhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_double *)compan,ldh,(lapack_complex_double *)roots,(lapack_complex_double *)zz,ldz);
                        if (info) { fprintf(stderr,"error in poly2roots_z: lapacke decomposition failed\n"); return 1; }
                        for (size_t l=Ly; l>0u; --l, ++roots, Y+=2u*K-1u) { *Y++ = *roots++; *Y = *roots; }
                        roots -= 2u*Ly;
                    }
                }
                free(roots);
            }
        }
        free(compan);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
/*System includes*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

/*GSL includes*/
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_multimin.h>
#include <pthread.h>

/*User includes*/

#include "Subsample.h"

static char *usage[] = {"Subsample - Subsamples an emprical rarefraction curve by resampling\n",
                        "Required parameters:\n",
			"   -in   filename      parameter file  \n",
			"   -i    integer       sample to subsample from\n",
			"   -r    integer       number of samples to generate\n",
			"   -n    integer       new sample size\n",
                        "Optional:\n",
			"   -seed long          seed random number generator\n",
			"   -v                  verbose\n"};

static int  nLines   = 9;

static int  verbose  = FALSE;

void setCP(double* adP, double *adC, int *anN, int nS, int nL)
{
  int i = 0;

  for(i = 0; i < nS; i++){
    adP[i] = ((double) anN[i])/((double) nL);
  }

  adC[i] = 0;
  for(i = 1; i <= nS; i++){
    adC[i] = adP[i - 1] + adC[i - 1];
  }
}

int sampleSpecies(gsl_rng *ptGSLRNG, double *adC, int nS)
{
  double dRand = gsl_rng_uniform(ptGSLRNG);
  int    nRet  = ((double) nS)*dRand;

  while(dRand <= adC[nRet]) nRet--;

  while(dRand > adC[nRet + 1]) nRet++;

  return nRet;
}

void writeSample(t_Params *ptParams, int nSample, int nR, int nS, int* anF)
{
  FILE *ofp = NULL;
  int i   = 0;
  char szFileName[MAX_LINE_LENGTH];

  sprintf(szFileName, "%s_%d_%d.dat",ptParams->szOutFileStub,nR, nSample);

  ofp = fopen(szFileName, "w");

  for(i = 0; i < nS; i++){
    fprintf(ofp,"%d %d\n",i, anF[i]);
  }

  fclose(ofp);
}

int main(int argc, char* argv[])
{
  int  i = 0, j = 0, r = 0;
  int  nCount = 0, nRS = 0, nD = 0;
  t_Params tParams;
  int      nS = 0, nL = 0;
  double   *adP = NULL;
  double   *adC = NULL;
  int      *anN = NULL;
  int      **aanF = NULL;
  gsl_rng            *ptGSLRNG = NULL;
  const gsl_rng_type *ptGSLRNGType = NULL;

  gsl_rng_env_setup();

  gsl_set_error_handler_off();

  /*get command line params*/
  getCommandLineParams(&tParams, argc, argv);

  ptGSLRNGType = gsl_rng_default;
  ptGSLRNG     = gsl_rng_alloc(ptGSLRNGType);

  gsl_rng_set(ptGSLRNG, tParams.lSeed);

  /*read in abundance distribution*/
  readAbundanceData(tParams.nI, tParams.szInputFile,&anN,&nS);

  for(i = 0; i < nS; i++){
    nL += anN[i];
  }

  adP = (double *) malloc(sizeof(double)*nS);
  adC = (double *) malloc(sizeof(double)*(nS + 1));
  aanF = (int    **) malloc(sizeof(int*)*nS);
  for(i = 0; i < nS; i++){
    aanF[i] = (int *) malloc(tParams.nR*sizeof(int));
  }

  setCP(adP, adC, anN, nS, nL);

  while(r < tParams.nR){

    for(i = 0; i < nS; i++){
      aanF[i][r] = 0;
    }

    for(i = 0; i < tParams.nN; i++){
      nD = sampleSpecies(ptGSLRNG, adC, nS);

      aanF[nD][r]++;
    }/*sample*/

    r++;
  }
  gsl_rng_free(ptGSLRNG);

  printf("OTU,");
  for(r = 0; r < tParams.nR - 1; r++){
    printf("r%d,",r);
  }
  printf("r%d\n",r);

  for(i = 0; i < nS; i++){
    printf("C%d,",i);
    for(r = 0; r < tParams.nR - 1; r++){
      printf("%d,",aanF[i][r]);
    }
    printf("%d\n",aanF[i][r]);
  }

  for(i = 0; i < nS; i++){
    free(aanF[i]);
  }
  free(aanF);
  free(anN);
  free(adC);
  free(adP);

  exit(EXIT_SUCCESS);
}

void writeUsage(FILE* ofp)
{
  int i = 0;
  char *line;

  for(i = 0; i < nLines; i++){
    line = usage[i];
    fputs(line,ofp);
  }
}

char *extractParameter(int argc, char **argv, char *param,int when)
{
  int i = 0;

  while((i < argc) && (strcmp(param,argv[i]))){
    i++;
  }

  if(i < argc - 1){
    return(argv[i + 1]);
  }

  if((i == argc - 1) && (when == OPTION)){
    return "";
  }

  if(when == ALWAYS){
    fprintf(stdout,"Can't find asked option %s\n",param);
  }

  return (char *) NULL;
}

void getCommandLineParams(t_Params *ptParams,int argc,char *argv[])
{
  char *szTemp = NULL;
  char *cError = NULL;

  /*get parameter file name*/
  ptParams->szInputFile  = extractParameter(argc,argv, INPUT_FILE,ALWAYS);
  if(ptParams->szInputFile == NULL)
    goto error;

  /*get out file stub*/
  //ptParams->szOutFileStub  = extractParameter(argc,argv,OUT_FILE_STUB,ALWAYS);
  //if(ptParams->szOutFileStub == NULL)
    //goto error;

  szTemp  = extractParameter(argc,argv,SAMPLE,OPTION);
  if(szTemp != NULL){
    ptParams->nI = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nI = 1;
  }

  szTemp  = extractParameter(argc,argv,N_SAMPLES,OPTION);
  if(szTemp != NULL){
    ptParams->nR = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nR = 10;
  }

  szTemp  = extractParameter(argc,argv,N_SIZE,OPTION);
  if(szTemp != NULL){
    ptParams->nN = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nN = 1000;
  }

  /*get out file stub*/
  szTemp  = extractParameter(argc,argv,SEED,OPTION);
  if(szTemp != NULL){
    ptParams->lSeed = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->lSeed = 0;
  }

  /*verbosity*/
  szTemp = extractParameter(argc, argv, VERBOSE, OPTION);
  if(szTemp != NULL){
    verbose = TRUE;
  }

  return;

 error:
  writeUsage(stdout);
  exit(EXIT_FAILURE);
}

void readAbundanceData(int nI, const char *szFile, int **panN, int *pnS)
{
  int  *anN = NULL;
  int  i = 0, j = 0, nS = 0;
  char szLine[MAX_LINE_LENGTH];
  FILE* ifp = NULL;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    fgets(szLine, MAX_LINE_LENGTH, ifp);

    while(fgets(szLine, MAX_LINE_LENGTH, ifp) != NULL){
    	nS++;
    }
    fclose(ifp);

    ifp = fopen(szFile, "r");
    fgets(szLine, MAX_LINE_LENGTH, ifp);
    anN = (int *) malloc(nS*sizeof(int));

    for(i = 0; i < nS; i++){

      fgets(szLine, MAX_LINE_LENGTH, ifp);

      szTok = strtok(szLine, DELIM);

      for(j = 0; j < nI; j++){
 	szTok = strtok(NULL, DELIM);
      }
      anN[i] = strtol(szTok,&pcError,10);

      if(*pcError != '\0'){
	goto formatError;
      }
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  (*pnS) = nS;
  (*panN) = anN;
  return;

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

int compare_doubles(const void* a, const void* b)
{
  double* arg1 = (double *) a;
  double* arg2 = (double *) b;
  if( *arg1 < *arg2 ) return -1;
  else if( *arg1 == *arg2 ) return 0;
  else return 1;
}

max_stars_repo_name: chrisquince/DiversityEstimates

max_issues_repo_name: chrisquince/DiversityEstimates

max_forks_repo_name: chrisquince/DiversityEstimates


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

#include <config.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

gsl_spline2d *
gsl_spline2d_alloc(const gsl_interp2d_type * T, size_t xsize, size_t ysize)
{
  double * array_mem;
  gsl_spline2d * interp;

  if (xsize < T->min_size || ysize < T->min_size)
    {
      GSL_ERROR_NULL("insufficient number of points for interpolation type", GSL_EINVAL);
    }

  interp = calloc(1, sizeof(gsl_spline2d));
  if (interp == NULL)
    {
      GSL_ERROR_NULL("failed to allocate space for gsl_spline2d struct", GSL_ENOMEM);
    }

  interp->interp_object.type = T;
  interp->interp_object.xsize = xsize;
  interp->interp_object.ysize = ysize;
  if (interp->interp_object.type->alloc == NULL)
    {
      interp->interp_object.state = NULL;
    }
  else
    {
      interp->interp_object.state = interp->interp_object.type->alloc(xsize, ysize);
      if (interp->interp_object.state == NULL)
        {
          gsl_spline2d_free(interp);
          GSL_ERROR_NULL("failed to allocate space for gsl_spline2d state", GSL_ENOMEM);
        }
    }

  /*
   * Use one contiguous block of memory for all three data arrays.
   * That way the code fails immediately if there isn't sufficient space for everything,
   * rather than allocating one or two and then having to free them.
   */
  array_mem = (double *)calloc(xsize + ysize + xsize * ysize, sizeof(double));
  if (array_mem == NULL)
    {
      gsl_spline2d_free(interp);
      GSL_ERROR_NULL("failed to allocate space for data arrays", GSL_ENOMEM);
    }

  interp->xarr = array_mem;
  interp->yarr = array_mem + xsize;
  interp->zarr = array_mem + xsize + ysize;

  return interp;
} /* gsl_spline2d_alloc() */

int
gsl_spline2d_init(gsl_spline2d * interp, const double xarr[],
                  const double yarr[], const double zarr[],
                  size_t xsize, size_t ysize)
{
  int status = gsl_interp2d_init(&(interp->interp_object), xarr, yarr, zarr, xsize, ysize);

  memcpy(interp->xarr, xarr, xsize * sizeof(double));
  memcpy(interp->yarr, yarr, ysize * sizeof(double));
  memcpy(interp->zarr, zarr, xsize * ysize * sizeof(double));

  return status;
} /* gsl_spline2d_init() */

void
gsl_spline2d_free(gsl_spline2d * interp)
{
  RETURN_IF_NULL(interp);

  if (interp->interp_object.type->free)
    interp->interp_object.type->free(interp->interp_object.state);

  /*
   * interp->xarr points to the beginning of one contiguous block of memory
   * that holds interp->xarr, interp->yarr, and interp->zarr. So it all gets
   * freed with one call. cf. gsl_spline2d_alloc() implementation
   */
  if (interp->xarr)
    free(interp->xarr);

  free(interp);
} /* gsl_spline2d_free() */

double
gsl_spline2d_eval(const gsl_spline2d * interp, const double x,
                  const double y, gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  return gsl_interp2d_eval(&(interp->interp_object), interp->xarr, interp->yarr,
                           interp->zarr, x, y, xa, ya);
}

int
gsl_spline2d_eval_e(const gsl_spline2d * interp, const double x,
                    const double y, gsl_interp_accel * xa, gsl_interp_accel * ya,
                    double * z)
{
  return gsl_interp2d_eval_e(&(interp->interp_object), interp->xarr, interp->yarr,
                             interp->zarr, x, y, xa, ya, z);
}

double
gsl_spline2d_eval_extrap(const gsl_spline2d * interp, const double x,
                         const double y, gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  return gsl_interp2d_eval_extrap(&(interp->interp_object), interp->xarr, interp->yarr,
                                  interp->zarr, x, y, xa, ya);
}

int
gsl_spline2d_eval_extrap_e(const gsl_spline2d * interp, const double x,
                           const double y, gsl_interp_accel * xa, gsl_interp_accel * ya,
                           double * z)
{
  return gsl_interp2d_eval_extrap_e(&(interp->interp_object), interp->xarr, interp->yarr,
                                    interp->zarr, x, y, xa, ya, z);
}

double
gsl_spline2d_eval_deriv_x(const gsl_spline2d * interp, const double x,
                          const double y, gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  return gsl_interp2d_eval_deriv_x(&(interp->interp_object), interp->xarr, interp->yarr,
                                   interp->zarr, x, y, xa, ya);
}

int
gsl_spline2d_eval_deriv_x_e(const gsl_spline2d * interp, const double x,
                            const double y, gsl_interp_accel * xa, gsl_interp_accel * ya,
                            double * z)
{
  return gsl_interp2d_eval_deriv_x_e(&(interp->interp_object), interp->xarr, interp->yarr,
                                     interp->zarr, x, y, xa, ya, z);
}

double
gsl_spline2d_eval_deriv_y(const gsl_spline2d * interp, const double x,
                          const double y, gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  return gsl_interp2d_eval_deriv_y(&(interp->interp_object), interp->xarr, interp->yarr,
                                   interp->zarr, x, y, xa, ya);
}

int
gsl_spline2d_eval_deriv_y_e(const gsl_spline2d * interp, const double x,
                            const double y, gsl_interp_accel * xa, gsl_interp_accel * ya,
                            double * z)
{
  return gsl_interp2d_eval_deriv_y_e(&(interp->interp_object), interp->xarr, interp->yarr,
                                     interp->zarr, x, y, xa, ya, z);
}

double
gsl_spline2d_eval_deriv_xx(const gsl_spline2d * interp, const double x,
                           const double y, gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  return gsl_interp2d_eval_deriv_xx(&(interp->interp_object), interp->xarr, interp->yarr,
                                    interp->zarr, x, y, xa, ya);
}

int
gsl_spline2d_eval_deriv_xx_e(const gsl_spline2d * interp, const double x,
                             const double y, gsl_interp_accel * xa, gsl_interp_accel * ya,
                             double * z)
{
  return gsl_interp2d_eval_deriv_xx_e(&(interp->interp_object), interp->xarr, interp->yarr,
                                      interp->zarr, x, y, xa, ya, z);
}

double
gsl_spline2d_eval_deriv_yy(const gsl_spline2d * interp, const double x,
                           const double y, gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  return gsl_interp2d_eval_deriv_yy(&(interp->interp_object), interp->xarr, interp->yarr,
                                    interp->zarr, x, y, xa, ya);
}

int
gsl_spline2d_eval_deriv_yy_e(const gsl_spline2d * interp, const double x,
                             const double y, gsl_interp_accel * xa, gsl_interp_accel * ya,
                             double * z)
{
  return gsl_interp2d_eval_deriv_yy_e(&(interp->interp_object), interp->xarr, interp->yarr,
                                      interp->zarr, x, y, xa, ya, z);
}

double
gsl_spline2d_eval_deriv_xy(const gsl_spline2d * interp, const double x,
                           const double y, gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  return gsl_interp2d_eval_deriv_xy(&(interp->interp_object), interp->xarr, interp->yarr,
                                    interp->zarr, x, y, xa, ya);
}

int
gsl_spline2d_eval_deriv_xy_e(const gsl_spline2d * interp, const double x,
                             const double y, gsl_interp_accel * xa, gsl_interp_accel * ya,
                             double * z)
{
  return gsl_interp2d_eval_deriv_xy_e(&(interp->interp_object), interp->xarr, interp->yarr,
                                      interp->zarr, x, y, xa, ya, z);
}

size_t
gsl_spline2d_min_size(const gsl_spline2d * interp)
{
  return gsl_interp2d_min_size(&(interp->interp_object));
}

const char *
gsl_spline2d_name(const gsl_spline2d * interp)
{
  return gsl_interp2d_name(&(interp->interp_object));
}

int
gsl_spline2d_set(const gsl_spline2d * interp, double zarr[],
                 const size_t i, const size_t j, const double z)
{
  return gsl_interp2d_set(&(interp->interp_object), zarr, i, j, z);
} /* gsl_spline2d_set() */

double
gsl_spline2d_get(const gsl_spline2d * interp, const double zarr[],
                 const size_t i, const size_t j)
{
  return gsl_interp2d_get(&(interp->interp_object), zarr, i, j);
} /* gsl_spline2d_get() */

max_stars_repo_name: igormcoelho/optstats

max_issues_repo_name: igormcoelho/optstats

max_forks_repo_name: igormcoelho/optstats