text:
/* bspline/gsl_bspline.h
 *
 * Copyright (C) 2006 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_BSPLINE_H__
#define __GSL_BSPLINE_H__

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

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

__BEGIN_DECLS

typedef struct
{
  size_t k;      /* spline order */
  size_t km1;    /* k - 1 (polynomial order) */
  size_t l;      /* number of polynomial pieces on interval */
  size_t nbreak; /* number of breakpoints (l + 1) */
  size_t n;      /* number of bspline basis functions (l + k - 1) */

  gsl_vector *knots;  /* knots vector */
  gsl_vector *deltal; /* left delta */
  gsl_vector *deltar; /* right delta */
  gsl_vector *B;      /* temporary spline results */
} gsl_bspline_workspace;

gsl_bspline_workspace *
gsl_bspline_alloc(const size_t k, const size_t nbreak);

void gsl_bspline_free(gsl_bspline_workspace *w);

size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * w);
size_t gsl_bspline_order (gsl_bspline_workspace * w);
size_t gsl_bspline_nbreak (gsl_bspline_workspace * w);
double gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w);

int
gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w);

int gsl_bspline_knots_uniform(const double a, const double b,
                              gsl_bspline_workspace *w);

int
gsl_bspline_eval(const double x, gsl_vector *B,
                 gsl_bspline_workspace *w);

__END_DECLS

#endif /* __GSL_BSPLINE_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef GSL_INTERFACE_H

#ifdef HAVE_GSL
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

namespace GSL{
	static std::string errMsg;
	static int gsl_errno;
	void error_handler(const char* reason, const char* file, int line, int gsl_errno){
		GSL::gsl_errno=gsl_errno;
		errMsg=file+std::string(":")+std::to_string(line)+": "+reason;
	}

	void reset_gsl_error(){ gsl_errno=0; }

	///\class
	///\brief Container for GSL workspace to be use with the integrators.
	class IntegrateWorkspace {
	private:

	public:
		gsl_integration_workspace* ws;
		IntegrateWorkspace(size_t limit) {
			ws=gsl_integration_workspace_alloc(limit);
		}
		~IntegrateWorkspace() {
			gsl_integration_workspace_free(ws);
		}
	};

	///\brief One dimensional integral using GSL.
	/// @param ws GSL integration workspace.
	/// @param f Function to integrate.
	/// @param a Lower integration limit.
	/// @param b Upper integration limit.
	/// @param acc Accuracy parameter.
	/// @param max_iter Maximum number of iterations to perform the integral.
	template<typename FunctionType>
	double integrate(IntegrateWorkspace& ws, FunctionType&& f, double a, double b, double acc=1e-7, unsigned int max_iter=10000){
		//precondition copied from gsl-1.16/integration/qag.c:122
		//assuming epsabs=0 and epsrel=acc
		if ((acc < 50 * GSL_DBL_EPSILON || acc < 0.5e-28))
			throw std::runtime_error("Specified tolerance too small");

		using FPtr=decltype(&f);
		double (*wrapper)(double,void*)=[](double x, void* params)->double{
			auto& f=*static_cast<FPtr>(params);
			return(f(x));
		};

		double result, error;
		gsl_function F;
		F.function = wrapper;
		F.params = &f;

		gsl_integration_qag(&F, a, b, 0, acc, max_iter, GSL_INTEG_GAUSS15, ws.ws, &result, &error);

		return(result);
	}

	///\brief One dimensional integral using GSL.
	/// @param f Function to integrate.
	/// @param a Lower integration limit.
	/// @param b Upper integration limit.
	/// @param acc Accuracy parameter.
	/// @param max_iter Maximum number of iterations to perform the integral.
	template<typename FunctionType>
	double integrate(FunctionType&& f, double a, double b, double acc=1e-7, unsigned int max_iter=10000, size_t memory_alloc=10000){
		IntegrateWorkspace ws(memory_alloc);
		return integrate(ws, f, a, b, acc, max_iter);
	}
} //namespace GSL
#endif //HAVE_GSL

#endif //GSL_INTERFACE_H

ext: h

max_stars_repo_name: cnweaver/AdaptiveQuad

max_issues_repo_name: cnweaver/AdaptiveQuad

max_forks_repo_name: cnweaver/AdaptiveQuad


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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_histogram2d.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>

#define N 107
#define M 239

static void my_error_handler (const char *reason, const char *file,
                              int line, int err);
static int status = 0;

void
test2d_trap (void)
{
  gsl_histogram2d *h;
  double result, lower, upper;
  size_t i, j;

  gsl_set_error_handler (&my_error_handler);

  gsl_ieee_env_setup ();

  status = 0;
  h = gsl_histogram2d_calloc (0, 10);
  gsl_test (!status, "gsl_histogram_calloc traps zero-width histogram");
  gsl_test (h != 0,
            "gsl_histogram2d_calloc returns NULL for zero-width histogram");

  status = 0;
  h = gsl_histogram2d_calloc (10, 0);
  gsl_test (!status, "gsl_histogram_calloc traps zero-length histogram");
  gsl_test (h != 0,
            "gsl_histogram2d_calloc returns NULL for zero-length histogram");

  status = 0;
  h = gsl_histogram2d_calloc_uniform (0, 10, 0.0, 1.0, 0.0, 1.0);
  gsl_test (!status,
            "gsl_histogram2d_calloc_uniform traps zero-width histogram");
  gsl_test (h != 0,
    "gsl_histogram2d_calloc_uniform returns NULL for zero-width histogram");

  status = 0;
  h = gsl_histogram2d_calloc_uniform (10, 0, 0.0, 1.0, 0.0, 1.0);
  gsl_test (!status,
            "gsl_histogram2d_calloc_uniform traps zero-length histogram");
  gsl_test (h != 0,
   "gsl_histogram2d_calloc_uniform returns NULL for zero-length histogram");

  status = 0;
  h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 1.0, 1.0);
  gsl_test (!status,
            "gsl_histogram2d_calloc_uniform traps equal endpoints");
  gsl_test (h != 0,
         "gsl_histogram2d_calloc_uniform returns NULL for equal endpoints");

  status = 0;
  h = gsl_histogram2d_calloc_uniform (10, 10, 1.0, 1.0, 0.0, 1.0);
  gsl_test (!status,
            "gsl_histogram2d_calloc_uniform traps equal endpoints");
  gsl_test (h != 0,
         "gsl_histogram2d_calloc_uniform returns NULL for equal endpoints");

  status = 0;
  h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 2.0, 1.0);
  gsl_test (!status,
            "gsl_histogram2d_calloc_uniform traps invalid range");
  gsl_test (h != 0,
            "gsl_histogram2d_calloc_uniform returns NULL for invalid range");

  status = 0;
  h = gsl_histogram2d_calloc_uniform (10, 10, 2.0, 1.0, 0.0, 1.0);
  gsl_test (!status,
            "gsl_histogram2d_calloc_uniform traps invalid range");
  gsl_test (h != 0,
            "gsl_histogram2d_calloc_uniform returns NULL for invalid range");

  h = gsl_histogram2d_calloc_uniform (N, M, 0.0, 1.0, 0.0, 1.0);

  status = gsl_histogram2d_accumulate (h, 1.0, 0.0, 10.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps x at xmax");

  status = gsl_histogram2d_accumulate (h, 2.0, 0.0, 100.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps x above xmax");

  status = gsl_histogram2d_accumulate (h, -1.0, 0.0, 1000.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps x below xmin");

  status = gsl_histogram2d_accumulate (h, 0.0, 1.0, 10.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps y at ymax");

  status = gsl_histogram2d_accumulate (h, 0.0, 2.0, 100.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps y above ymax");

  status = gsl_histogram2d_accumulate (h, 0.0, -1.0, 1000.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps y below ymin");

  status = gsl_histogram2d_increment (h, 1.0, 0.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps x at xmax");

  status = gsl_histogram2d_increment (h, 2.0, 0.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps x above xmax");

  status = gsl_histogram2d_increment (h, -1.0, 0.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps x below xmin");

  status = gsl_histogram2d_increment (h, 0.0, 1.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps y at ymax");

  status = gsl_histogram2d_increment (h, 0.0, 2.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps y above ymax");

  status = gsl_histogram2d_increment (h, 0.0, -1.0);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps y below ymin");

  result = gsl_histogram2d_get (h, N, 0);
  gsl_test (result != 0, "gsl_histogram2d_get traps x index at nx");

  result = gsl_histogram2d_get (h, N + 1, 0);
  gsl_test (result != 0, "gsl_histogram2d_get traps x index above nx");

  result = gsl_histogram2d_get (h, 0, M);
  gsl_test (result != 0, "gsl_histogram2d_get traps y index at ny");

  result = gsl_histogram2d_get (h, 0, M + 1);
  gsl_test (result != 0, "gsl_histogram2d_get traps y index above ny");

  status = gsl_histogram2d_get_xrange (h, N, &lower, &upper);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_xrange traps index at nx");

  status = gsl_histogram2d_get_xrange (h, N + 1, &lower, &upper);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_xrange traps index above nx");

  status = gsl_histogram2d_get_yrange (h, M, &lower, &upper);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_yrange traps index at ny");

  status = gsl_histogram2d_get_yrange (h, M + 1, &lower, &upper);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_yrange traps index above ny");

  status = 0;
  gsl_histogram2d_find (h, -0.01, 0.0, &i, &j);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps x below xmin");

  status = 0;
  gsl_histogram2d_find (h, 1.0, 0.0, &i, &j);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps x at xmax");

  status = 0;
  gsl_histogram2d_find (h, 1.1, 0.0, &i, &j);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps x above xmax");

  status = 0;
  gsl_histogram2d_find (h, 0.0, -0.01, &i, &j);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps y below ymin");

  status = 0;
  gsl_histogram2d_find (h, 0.0, 1.0, &i, &j);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps y at ymax");

  status = 0;
  gsl_histogram2d_find (h, 0.0, 1.1, &i, &j);
  gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps y above ymax");

  gsl_histogram2d_free (h);
}

static void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  if (0)
    printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
  status = 1;
}

ext: c

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
/* hmm.c v0.0

 	Andrew D. Kern	8/24/05
*/

#include "hmm.h"
#include "adkGSL.h"
#include <math.h>
#include <assert.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_randist.h>

/* returns a pointer to a new instance of an HMM.
Currently an HMM is defined by a transition matrix (pMatrix),
and a vector representing the starting probabilities.
use of initHMM() fills the log score matrices. 	*/

HMM* newHMM(gsl_matrix *pMatrix, gsl_vector *piStart, gsl_vector *piEnd, gsl_vector *otherData ) {
  HMM *hmm = (HMM*)malloc(sizeof(HMM));
  int i;

  hmm->transition_matrix = pMatrix;
  hmm->piStart = piStart;
  hmm->piEnd = piEnd;
  hmm->nstates = pMatrix->size1;
  hmm->transition_score_matrix = NULL;
  hmm->piStart_scores = NULL;
  hmm->piEnd_scores = NULL;
  hmm->otherData = otherData;

  /* if piStart are NULL, make them uniform */
  if (piStart == NULL) {
    hmm->piStart = gsl_vector_alloc(pMatrix->size1);
    for (i = 0; i < pMatrix->size1; i++)
      gsl_vector_set(hmm->piStart, i, 1.0/pMatrix->size1);
  }
  if (piEnd == NULL) {
    hmm->piEnd = gsl_vector_alloc(pMatrix->size1);
    for (i = 0; i < pMatrix->size1; i++)
      gsl_vector_set(hmm->piEnd, i, 1.0/pMatrix->size1);
  }

  return hmm;
}

/* initialize an HMM, this makes sure that things are kosher,
it checks that the row sums of the pMatrix == 1, that the sum
of the piStart vector == 1. then it initializes scores */

void initHMM(HMM *hmm){
	gsl_matrix *pMatrixLogs;
	gsl_vector *piStartLogs, *piEndLogs;
	int i, j;
	double piEndSum, piStartSum, prob;

	//check piStartSum == 1
	piStartSum = gsl_vector_sum(hmm->piStart, hmm->nstates);
	assert(piStartSum == 1.0);

	//make log transformed start probs
	piStartLogs = gsl_vector_alloc(hmm->nstates);
	for (i = 0; i < hmm->nstates; i++){
		prob = gsl_vector_get(hmm->piStart, i);
		gsl_vector_set(piStartLogs, i, log(prob));
		}

	//set hmm piEnd_scores
	hmm->piStart_scores = piStartLogs;

	//check piEndSum == 1
	piEndSum = gsl_vector_sum(hmm->piEnd, hmm->nstates);
	assert(piEndSum == 1.0);

	//make log transformed start probs
	piEndLogs = gsl_vector_alloc(hmm->nstates);
	for (i = 0; i < hmm->nstates; i++){
		prob = gsl_vector_get(hmm->piEnd, i);
		gsl_vector_set(piEndLogs, i, log(prob));
		}

	//set hmm piEnd_scores
	hmm->piEnd_scores = piEndLogs;

	//check row sums for markov matrix
	for(i = 0; i < hmm->nstates; i++){
		assert(gsl_matrix_row_sum(hmm->transition_matrix, i, hmm->nstates) == 1.0);
		}
	//make log transformed transition matrix
	pMatrixLogs = gsl_matrix_alloc(hmm->nstates, hmm->nstates);
	for (i = 0; i < hmm->nstates; i++){
		for (j = 0; j < hmm->nstates; j++){
			prob = gsl_matrix_get(hmm->transition_matrix, i, j);
			gsl_matrix_set(pMatrixLogs, i, j, log(prob));
			}
		}

	//set hmm transition_score_matrix
	hmm->transition_score_matrix = pMatrixLogs;
}

/* frees an entire hmm */
void freeHMM(HMM *hmm){

	if (hmm->transition_matrix)
		gsl_matrix_free(hmm->transition_matrix);
	if (hmm->transition_score_matrix)
		gsl_matrix_free(hmm->transition_score_matrix);
	if (hmm->piStart)
		gsl_vector_free(hmm->piStart);
	if (hmm->piStart_scores)
		gsl_vector_free(hmm->piStart_scores);
	if (hmm->piEnd)
		gsl_vector_free(hmm->piEnd);
	if (hmm->piEnd_scores)
		gsl_vector_free(hmm->piEnd_scores);
	if (hmm->otherData)
		gsl_vector_free(hmm->otherData);
	free(hmm);
}

/*output hmm values: nstates, transition matrix, piStarts */
void printHMM(HMM *hmm){
  printf("hmm output\nstates: %d \ntransition matrix:\n", hmm->nstates);
  gsl_matrix_fprintf(stdout, hmm->transition_matrix, "%f");
  printf("piStart vector:\n");
  gsl_vector_fprintf(stdout, hmm->piStart, "%f");
  printf("piEnd vector:\n");
  gsl_vector_fprintf(stdout, hmm->piEnd, "%f");
  if (hmm->transition_score_matrix != NULL){
    printf("transition matrix (logs):\n");
    gsl_matrix_fprintf(stdout, hmm->transition_score_matrix, "%f");
  }
  if (hmm->piStart_scores != NULL){
    printf("piStart vector (logs):\n");
    gsl_vector_fprintf(stdout, hmm->piStart_scores, "%f");
  }
  if (hmm->piEnd_scores != NULL){
    printf("piEnd vector (logs):\n");
    gsl_vector_fprintf(stdout, hmm->piEnd_scores, "%f");
  }
  if (hmm->otherData != NULL){
    //    printf("other data:\n");
    //gsl_vector_fprintf(stdout, hmm->otherData, "%f");
  }

}

/*output transition matrix of hmm values to file */
void printTransitions(HMM *h, char *outfileName){
  FILE *outfile;
  int i,j;

   outfile = fopen(outfileName, "w");
   if (outfile == NULL){
     fprintf(stderr,"Error opening outfile! ARRRRR!!!!\n");
     exit(1);
   }
   for(i = 0; i < h->nstates; i++){
     for(j = 0; j < h->nstates; j++){
       fprintf(outfile,"%f\t",gsl_matrix_get(h->transition_matrix,i,j));
     }
     fprintf(outfile,"\n");
   }
   fclose(outfile);
}

/* this normalizes the transition matrix so that rows sum to one */
void hmm_normalize_transitions(HMM *hmm){
  int i, j;
  double rowSum, tempVal;

  for(i=0; i < hmm->nstates; i++){
    rowSum = 0;
    for(j=0; j < hmm->nstates; j++){
      rowSum += gsl_matrix_get(hmm->transition_matrix,i,j);
    }
    if(rowSum == 0){
      gsl_matrix_set(hmm->transition_matrix, i, i, 1.0);
    }
    else{
      for(j= 0; j < hmm->nstates;j++){
	tempVal =  gsl_matrix_get(hmm->transition_matrix, i, j) / rowSum;
	gsl_matrix_set(hmm->transition_matrix, i, j, tempVal);
      }
    }
  }
}

void hmm_logify_transitions(HMM *hmm){
  int i, j;
  double prob;
  //set log transition scores
  for (i = 0; i < hmm->nstates; i++){
    for (j = 0; j < hmm->nstates; j++){
      prob = gsl_matrix_get(hmm->transition_matrix, i, j);
      gsl_matrix_set(hmm->transition_score_matrix, i, j, log(prob));
    }
  }
  //set log start and end scores
  for (i = 0; i < hmm->nstates; i++){
    prob = gsl_vector_get(hmm->piStart, i);
    gsl_vector_set(hmm->piStart_scores, i, log(prob));
    prob = gsl_vector_get(hmm->piEnd, i);
    gsl_vector_set(hmm->piEnd_scores, i, log(prob));
  }
}

/* forward algorithm- this fills a matrix (alphas) of Nstates * L observations
with the forward log probabilities. a matrix of emission probs (log!!!) with the same
dimensions as the forward matrix must be allocated and passed to this function.
returns probFor. */

double forwardAlg(HMM *hmm, gsl_matrix *alphas, gsl_matrix *emisLogs, int length, int transitionPowerFlag){
	int i, j, power, last;
	double prob, emit, prevSum, pFor;
	gsl_vector_view prevCol;
	gsl_vector *tempProbs;
	gsl_matrix *logTransitions;

	assert(length > 0 && hmm != NULL);
	assert(alphas->size1 == hmm->nstates && alphas->size2 == length);
	assert(emisLogs->size1 == hmm->nstates && emisLogs->size2 == length);

//initialization
	for (i = 0; i < hmm->nstates; i++){
		prob = gsl_vector_get(hmm->piStart_scores, i);
		emit = gsl_matrix_get(emisLogs, i, 0);
		gsl_matrix_set(alphas, i, 0, (prob + emit));
	}

//are we using powers of the transition matrix? if so check to see we have location data
	if (transitionPowerFlag){
		assert(hmm->otherData != NULL);
	//set "last" equal to the first location (stored in otherData) here's a totally non general assumption... FIX ME!!!
		last = gsl_vector_get(hmm->otherData,0);
	}
	else {
		last = 0; /* not used */
	}

//fill em up!
	for (j = 1; j < length; j++){			//iterate over obs
		for (i = 0; i < hmm->nstates; i++){		//iterate over states
			emit = gsl_matrix_get(emisLogs, i, j);
			/* are we taking powers here? */
			if (transitionPowerFlag){
			//get power from data
				power = gsl_vector_get(hmm->otherData,j) - last;

				//get power of matrix in log scale
				logTransitions = gsl_matrix_power_logs(hmm->transition_matrix, power);
				prevSum = forwardSumPrevious(hmm, alphas, logTransitions, i, j);
				gsl_matrix_free(logTransitions);
			}
			/* no powers */
			else{
				prevSum = forwardSumPrevious(hmm, alphas, hmm->transition_score_matrix, i, j);
			//fprintf (stderr,"forwardAlg state %2d, index %2d, log.prevSum %g, emit %g\n",i,j,prevSum,emit);
			}
			gsl_matrix_set(alphas, i, j, (emit + prevSum));
		}
		if (transitionPowerFlag){
			last = gsl_vector_get(hmm->otherData,j);
		}
	}
//return P(x)
	tempProbs = gsl_vector_alloc(hmm->nstates);
	prevCol = gsl_matrix_column(alphas,  length - 1);
	for(i = 0; i < hmm->nstates; i++){
		gsl_vector_set(tempProbs, i, gsl_vector_get(&prevCol.vector, i));
	}
	pFor = log_sum(tempProbs);
	gsl_vector_free(tempProbs);
	return(pFor);
}

/* forward algorithm offset- identical to the normal forward algorithm except starts at an arbitrary
offset with respect to the observations. this is useful for getting probabilities from path subsets */

double forwardAlgOffset(HMM *hmm, gsl_matrix *alphas, gsl_matrix *emisLogs, int length, int offset, int transitionPowerFlag){
  int i, j, power, last;
  double prob, emit, prevSum, pFor;
  gsl_vector_view prevCol;
  gsl_vector *tempProbs;
  gsl_matrix *logTransitions;

  assert(length > 0 && hmm != NULL);
  assert(alphas->size1 == hmm->nstates && alphas->size2 == length);

  //initialization
  for (i = 0; i < hmm->nstates; i++){
    prob = gsl_vector_get(hmm->piStart_scores, i);
    emit = gsl_matrix_get(emisLogs, i, offset);
    gsl_matrix_set(alphas, i, 0, (prob + emit));
  }

  //are we using powers of the transition matrix? if so check to see we have location data
  if (transitionPowerFlag){
    assert(hmm->otherData != NULL);
    //set "last" equal to the first location (stored in otherData) here's a totally non general assumption... FIX ME!!!
    last = gsl_vector_get(hmm->otherData,offset);
  }
  else {
      last = 0; /* not used */
  }

  //fill em up!
  for (j = 1; j < length; j++){			//iterate over obs
    for (i = 0; i < hmm->nstates; i++){		//iterate over states
      emit = gsl_matrix_get(emisLogs, i, j + offset);
      /* are we taking powers here? */
      if (transitionPowerFlag){
	//get power from data
	power = gsl_vector_get(hmm->otherData,j + offset) - last;

	//get power of matrix in log scale
	logTransitions = gsl_matrix_power_logs(hmm->transition_matrix, power);
	prevSum = forwardSumPrevious(hmm, alphas, logTransitions, i, j);
	gsl_matrix_free(logTransitions);
      }
      /* no powers */
      else{
	prevSum = forwardSumPrevious(hmm, alphas, hmm->transition_score_matrix, i, j);
      }
      gsl_matrix_set(alphas, i, j, (emit + prevSum));
    }
    if (transitionPowerFlag){
	last = gsl_vector_get(hmm->otherData,j+offset);
    }
  }
  //return P(x)
  tempProbs = gsl_vector_alloc(hmm->nstates);
  prevCol = gsl_matrix_column(alphas,  length - 1);
  for(i = 0; i < hmm->nstates; i++){
    gsl_vector_set(tempProbs, i, gsl_vector_get(&prevCol.vector, i));
  }
  pFor = log_sum(tempProbs);
  gsl_vector_free(tempProbs);
  return(pFor);
}

/* backward algorithm- this fills a matrix (betas) of Nstates * L observations
with the backward log probabilities. a matrix of emission probs (log!!!) with the same
dimensions as the backward matrix must be allocated and passed to this function.
returns probBack. */

double backwardAlg(HMM *hmm, gsl_matrix *betas, gsl_matrix *emisLogs, int length, int transitionPowerFlag){
  int i, j, last, power;
  gsl_vector_view nextCol;
  gsl_vector *tempProbs;
  gsl_matrix *logTransitions;
  double pBack;

  assert(length > 0 && hmm != NULL);
  assert(betas->size1 == hmm->nstates && betas->size2 == length);
  assert(emisLogs->size1 == hmm->nstates && emisLogs->size2 == length);

  //initialization
  for (i = 0; i < hmm->nstates; i++){
    //prob = gsl_vector_get(hmm->piEnd_scores, i);
    //emit = gsl_matrix_get(emisLogs, i, length - 1);
    gsl_matrix_set(betas, i, length - 1, 0);		//initialize to prob = 1.0, log = 0
  }

  //are we using powers of the transition matrix?
  if (transitionPowerFlag){
    assert(hmm->otherData != NULL);
    //set "last" equal to the last location (stored in otherData) here's a totally non general assumption... FIX ME!!!
    last = gsl_vector_get(hmm->otherData,length - 1);
  }
  else {
      last = 0; /* not used */
  }

  //fill em up!
  for (j =  length - 2; j >=0; j--){			//iterate over obs
    for (i = 0; i < hmm->nstates; i++){		//iterate over states
      /* are we taking powers here? */
      if (transitionPowerFlag){
	//get power from data
	power = last - gsl_vector_get(hmm->otherData,j);

	//get power of matrix in log scale
	logTransitions = gsl_matrix_power_logs(hmm->transition_matrix, power);
	gsl_matrix_set(betas, i, j, backwardSumNext(hmm, betas,logTransitions,emisLogs, i, j));
	gsl_matrix_free(logTransitions);
      }
      else{
	gsl_matrix_set(betas, i, j, backwardSumNext(hmm, betas,hmm->transition_score_matrix,emisLogs, i, j));
      }
    }
    if (transitionPowerFlag){
	last = gsl_vector_get(hmm->otherData,j);
    }
  }
  //return P(x)
  tempProbs = gsl_vector_alloc(hmm->nstates);
  nextCol = gsl_matrix_column(betas, 0);
  for(i = 0; i < hmm->nstates; i++){
    gsl_vector_set(tempProbs, i, (gsl_vector_get(&nextCol.vector, i) + gsl_matrix_get(emisLogs, i, 0) + gsl_vector_get(hmm->piStart_scores, i)));
  }
  pBack = log_sum(tempProbs);
  gsl_vector_free(tempProbs);
  return(pBack);
}

/* posterior probs- this fills a matrix (postProbMat) of nstates * L observations with the
posterior prob (not log !!) of being in that state at time i given the observations and the model.
calls forwardAlg and backwardAlg, checks that prob of observations determined by
forward and backward are (nearly) equal, and also (nearly) equal the sum of each column
before normalization to postProbMat.
a matrix of emission probs (log!!!) of dimensions nstates * L observations must be supplied. */

void posteriorProbs(HMM *hmm, gsl_matrix *postProbMat, gsl_matrix *emisLogs, int length, int transitionPowerFlag){
	gsl_matrix *alphas, *betas;
	gsl_vector *tempVec;
	int i, j;
	double pFor, pBack, sumProb;

	//allocate matrices for forwards and backwards algs.
	alphas = gsl_matrix_alloc(hmm->nstates, length);
	betas = gsl_matrix_alloc(hmm->nstates, length);

	//run forwards/backwards
	pFor = forwardAlg(hmm, alphas, emisLogs, length, transitionPowerFlag);
	pBack = backwardAlg(hmm, betas, emisLogs, length, transitionPowerFlag);

	//check that Forward and Backward look okay
	//assert(fabs(pFor - pBack) < 0.01);

	//do the posterior probability thing- okay sir pull the switch
	tempVec = gsl_vector_alloc(hmm->nstates);
	for(j=0; j < length; j++){
		//first get sum for denominator
		for(i = 0; i < hmm->nstates; i++){
			gsl_vector_set(tempVec, i, (gsl_matrix_get(alphas,i,j) +  gsl_matrix_get(betas,i,j)));
		}
		sumProb = log_sum(tempVec);

		//check that sum looks like pFor
		assert(fabs(pFor - sumProb) < 0.01);

		//then calculate marginal and set value (not log)
		for(i = 0; i < hmm->nstates; i++){
			gsl_matrix_set(postProbMat, i, j, exp(gsl_matrix_get(alphas,i,j) +  gsl_matrix_get(betas,i,j) - sumProb));
		}
	}
	//cleanup
	gsl_matrix_free(alphas);
	gsl_matrix_free(betas);
	gsl_vector_free(tempVec);
}

/* posterior probs reduced- this fills a matrix of nstates * L observations with the
posterior probs as above except this time takes alphas and betas as args */

void posteriorProbsReduced(HMM *hmm, gsl_matrix *postProbMat, double pFor, gsl_matrix *alphas, gsl_matrix *betas, gsl_matrix *emisLogs, int length, int transitionPowerFlag){
  gsl_vector *tempVec;
  int i, j;
  double  sumProb;

  //do the posterior probability thing- okay sir pull the switch
  tempVec = gsl_vector_alloc(hmm->nstates);
  for(j=0; j < length; j++){
    //first get sum for denominator
    for(i = 0; i < hmm->nstates; i++){
      gsl_vector_set(tempVec, i, (gsl_matrix_get(alphas,i,j) +  gsl_matrix_get(betas,i,j)));
    }
    sumProb = log_sum(tempVec);

    //check that sum looks like pFor
      //assert(fabs(pFor - sumProb) < 0.01);

    //then calculate marginal and set value (not log)
    for(i = 0; i < hmm->nstates; i++){
      gsl_matrix_set(postProbMat, i, j, exp(gsl_matrix_get(alphas,i,j) +  gsl_matrix_get(betas,i,j) - sumProb));
    }
  }
  //cleanup
   gsl_vector_free(tempVec);
}

/* posterior probs reduced log- this fills a matrix of nstates * L observations with the
posterior probs as above except this time takes alphas and betas as args */

void posteriorProbsReducedLog(HMM *hmm, gsl_matrix *postProbMat, double pFor, gsl_matrix *alphas, gsl_matrix *betas, gsl_matrix *emisLogs, int length, int transitionPowerFlag){
  gsl_vector *tempVec;
  int i, j;
  double  sumProb;

  //do the posterior probability thing- okay sir pull the switch
  tempVec = gsl_vector_alloc(hmm->nstates);
  for(j=0; j < length; j++){
    //first get sum for denominator
    for(i = 0; i < hmm->nstates; i++){
      gsl_vector_set(tempVec, i, (gsl_matrix_get(alphas,i,j) +  gsl_matrix_get(betas,i,j)));
    }
    sumProb = log_sum(tempVec);

    //check that sum looks like pFor
    //assert(fabs(pFor - sumProb) < 0.01);

    //then calculate marginal and set value
    for(i = 0; i < hmm->nstates; i++){
      gsl_matrix_set(postProbMat, i, j, gsl_matrix_get(alphas,i,j) +  gsl_matrix_get(betas,i,j) - sumProb);
    }
  }
  //cleanup
   gsl_vector_free(tempVec);
}

/* meant to computer prob{state_i_t, state_j_t+1 | data, model} -- for restimating transitions, takes as an arg *jpMatArray an
array of gsl_matrix structs representing joint posterior probabilities -- log scale so need log posts */

void jointPosteriorProbsReduced(HMM *hmm, gsl_matrix **jpMatArray, double pFor, gsl_matrix *alphas, gsl_matrix *betas, gsl_matrix *emisLogs, int length, int transitionPowerFlag){

  int i, j, t;
  double  sumProb, val = 0;
  gsl_matrix *powTrans;

  //for each obs
  for(t=0; t < length - 1; t++){
    gsl_matrix_set_all(jpMatArray[t],0);
    //iterate i,j states
    for(i = 0; i < hmm->nstates; i++){
      for(j = 0; j < hmm->nstates; j++){
	//transition stuff
	if (transitionPowerFlag){

	  if (gsl_vector_get(hmm->otherData,t+1) -  gsl_vector_get(hmm->otherData, t) > 1){
	    powTrans = gsl_matrix_power_logs(hmm->transition_matrix, gsl_vector_get(hmm->otherData,t+1) \
					     -  gsl_vector_get(hmm->otherData, t));
	    val = gsl_matrix_get(alphas, i, t) + gsl_matrix_get(powTrans, i, j) \
	      + gsl_matrix_get(emisLogs, j, t+1) + gsl_matrix_get(betas, j, t+1) - pFor;
	    gsl_matrix_free(powTrans);
	  }
	  else{
	    val = gsl_matrix_get(alphas, i, t) + gsl_matrix_get(hmm->transition_score_matrix, i, j) \
	      + gsl_matrix_get(emisLogs, j, t+1) + gsl_matrix_get(betas, j, t+1) - pFor;
	  }
	}
	else{
	  val = gsl_matrix_get(alphas, i, t) + gsl_matrix_get(hmm->transition_score_matrix, i, j) \
		    + gsl_matrix_get(emisLogs, j, t+1) + gsl_matrix_get(betas, j, t+1) - pFor;
	}
	//set i,j value in gsl_matrix at t
	gsl_matrix_set(jpMatArray[t], i, j, val);
      }
    }
  }
}

/*restimateTransitions-- uses posteriors and joint posteriors to reestimate the transition matrix based on EM */
void reestimateTransitions(HMM *hmm, gsl_matrix *postProbMat, gsl_matrix **jpMatArray, double pFor, gsl_matrix *alphas, \
			   gsl_matrix *betas, gsl_matrix *emisLogs, int length, int transitionPowerFlag){
  double sum_gamma, sum_cosi, val;
  int i, j, t,k;
  gsl_vector *tmpVec1, *tmpVec2;
  gsl_matrix *tempTrans, *powTrans;
  gsl_vector *rowSums;

  /* this is the joint posterior method- currently not operational */
  /*
  //calculate joints, no reason to do it twice...
  //  jointPosteriorProbsReduced(hmm, jpMatArray, pFor, alphas, betas, emisLogs, length ,transitionPowerFlag);
  //allocate tmpVecs
  tmpVec1 = gsl_vector_alloc(length - 2);
  tmpVec2 = gsl_vector_alloc(length - 2);
  //iterate over states
  for(i = 0; i < hmm->nstates; i++){
    for(j = 0; j < hmm->nstates; j++){
      gsl_vector_set_all(tmpVec1,0);
      gsl_vector_set_all(tmpVec2,0);
      //sum up numerator and denom
      for(t = 1; t < length - 1; t++){
	//	gsl_vector_set(tmpVec1, t - 1, gsl_matrix_get(postProbMat,i,t));
	//gsl_vector_set(tmpVec2, t - 1, gsl_matrix_get(jpMatArray[t],i,j));
	sum_gamma += gsl_matrix_get(postProbMat,i,t);
	sum_cosi += gsl_matrix_get(jpMatArray[t],i,j);
      }

      //gsl_matrix_set(hmm->transition_matrix, i, j, exp(log_sum(tmpVec1) - log_sum(tmpVec2)));
      // hmm_logify_transitions(hmm);
      gsl_matrix_set(hmm->transition_matrix, i, j, sum_cosi / sum_gamma);
       gsl_matrix_set(hmm->transition_score_matrix, i, j, log(sum_cosi / sum_gamma));
    }
  }
  */

  /* older method */
  tempTrans = gsl_matrix_alloc(hmm->nstates,hmm->nstates);
  rowSums = gsl_vector_alloc(hmm->nstates);

  gsl_matrix_set_all(tempTrans,0.0);
  gsl_vector_set_all(rowSums, 0.0);
  for(i = 1; i < length - 1; i++){
    for(j = 0; j < hmm->nstates; j++){
      for(k = 0; k < hmm->nstates; k++){
	if (transitionPowerFlag){
	  if (gsl_vector_get(hmm->otherData,i + 1) -  gsl_vector_get(hmm->otherData, i ) > 1){
	    powTrans = gsl_matrix_power_logs(hmm->transition_matrix, \
					     gsl_vector_get(hmm->otherData,i + 1) -  gsl_vector_get(hmm->otherData, i));
	    val = exp(gsl_matrix_get(alphas, j, i) + gsl_matrix_get(powTrans, j, k) \
		      + gsl_matrix_get(emisLogs, k, i+1) + gsl_matrix_get(betas, k, i+1) - pFor);
	    gsl_matrix_free(powTrans);
	  }
	  else{
	    val = exp(gsl_matrix_get(alphas, j, i) + gsl_matrix_get(hmm->transition_score_matrix, j, k) \
		      + gsl_matrix_get(emisLogs, k, i+1) + gsl_matrix_get(betas, k, i+1) - pFor);
	  }
	}
	else{
	    val = exp(gsl_matrix_get(alphas, j, i) + gsl_matrix_get(hmm->transition_score_matrix, j, k) \
		      + gsl_matrix_get(emisLogs, k, i+1) + gsl_matrix_get(betas, k, i+1) - pFor);
	}
      	gsl_matrix_set(tempTrans, j, k, gsl_matrix_get(tempTrans,j,k) + val);
	gsl_vector_set(rowSums, j, gsl_vector_get(rowSums, j) + val);
      }
    }
  }

  //update transition scores
  for(j = 0; j < hmm->nstates; j++){
    for(k = 0; k < hmm->nstates; k++){
      val = gsl_matrix_get(tempTrans, j, k) / gsl_vector_get(rowSums, j);
      gsl_matrix_set(hmm->transition_matrix, j, k, val);
      gsl_matrix_set(hmm->transition_score_matrix, j, k, log(val));
    }
  }
  gsl_matrix_free(tempTrans);
  gsl_vector_free(rowSums);
}

/*restimateTransitions_vanilla-- uses posteriors to reestimate the transition matrix based on EM */
void reestimateTransitions_vanilla(HMM *hmm, double pFor, gsl_matrix *alphas, \
gsl_matrix *betas, gsl_matrix *emisLogs, int length){
	double val;
	int i, j,k;
	gsl_vector *rowSums;
	gsl_matrix *tempTrans;

/* older method */
  	tempTrans = gsl_matrix_alloc(hmm->nstates,hmm->nstates);
	gsl_matrix_set_all(tempTrans,0.0);
	rowSums = gsl_vector_alloc(hmm->nstates);
	gsl_vector_set_all(rowSums, 0.0);
	for(i = 1; i < length - 1; i++){
		for(j = 0; j < hmm->nstates; j++){
			for(k = 0; k < hmm->nstates; k++){
				val = exp(gsl_matrix_get(alphas, j, i) + gsl_matrix_get(hmm->transition_score_matrix, j, k) \
					+ gsl_matrix_get(emisLogs, k, i+1) + gsl_matrix_get(betas, k, i+1) - pFor);

				gsl_matrix_set(tempTrans, j, k, gsl_matrix_get(tempTrans,j,k) + val);
				gsl_vector_set(rowSums, j, gsl_vector_get(rowSums, j) + val);
			}
		}
	}

//update transition scores
	for(j = 0; j < hmm->nstates; j++){
		for(k = 0; k < hmm->nstates; k++){
			val = gsl_matrix_get(tempTrans, j, k) / gsl_vector_get(rowSums, j);
			gsl_matrix_set(hmm->transition_matrix, j, k, val);
			gsl_matrix_set(hmm->transition_score_matrix, j, k, log(val));
		}
	}
	gsl_matrix_free(tempTrans);
	gsl_vector_free(rowSums);
}

/* posteriorDecode-- fills a vector (posteriorPath) of length L observations
with the most probably state sequence based on posterior probabilities */
void posteriorDecode(HMM *hmm, gsl_vector *posteriorPath, gsl_matrix *posts, int length){
	int i, j, maxState;
	double maxProb;

	assert(length > 0 && hmm != NULL);
	assert(posteriorPath->size == length);
	assert(posts->size1 == hmm->nstates && posts->size2 == length);

	for(i=0; i < length; i++){
		maxProb = 0;
		maxState = 0;
		for(j=0;j < hmm->nstates;j++){
			if (gsl_matrix_get(posts, j, i) > maxProb){
				maxProb = gsl_matrix_get(posts, j, i);
				maxState = j;
			}
		}
		gsl_vector_set(posteriorPath,i,maxState);
	}

}

/* viterbi algorithm- this fills a vector (viterbiPath) of length L observations with
the most probable state sequence. it returns the prob of the most probable path.
a matrix of emission probs (log!!!) of dimensions nstates * L observations must be supplied. */

  double viterbi(HMM *hmm, gsl_vector *viterbiPath, gsl_matrix *emisLogs, int length, int transitionPowerFlag){
  gsl_matrix *deltas, *psi, *logTransitions;
	gsl_vector_view col;
	int i, j, lastState, tmp, last, power;
	double tempProb,prob, emit;
	int point;

	assert(length > 0 && hmm != NULL);
	assert(viterbiPath->size == length);
	assert(emisLogs->size1 == hmm->nstates && emisLogs->size2 == length);

	//if using transition powers, then check that we have other data
	if (transitionPowerFlag){
	  assert(hmm->otherData != NULL);
	  //initialize last
	  last = gsl_vector_get(hmm->otherData,0);
	}
	else {
	    last = 0; /* not used */
	}

	//allocate matrices for deltas (log prob of most prob path to stateN,obsL) and psis (the backpointers)
	deltas = gsl_matrix_alloc(hmm->nstates, length);
	psi    = gsl_matrix_alloc(hmm->nstates, length);

	//initialization
	for (i = 0; i < hmm->nstates; i++){
		prob = gsl_vector_get(hmm->piStart_scores, i);
		emit = gsl_matrix_get(emisLogs, i, 0);
		gsl_matrix_set(deltas, i, 0, (prob + emit));
		gsl_matrix_set(psi, i, 0, 0);
		}

	point = 666;
	//recursion
	for (j = 1; j < length; j++){
	  for(i = 0; i < hmm->nstates; i++){
	    /* transition powers? */
	    if (transitionPowerFlag){
	      power = gsl_vector_get(hmm->otherData,j) - last;
	      logTransitions = gsl_matrix_power_logs(hmm->transition_matrix, power);
	      tempProb = viterbiMaxPrevious(hmm, deltas, logTransitions, i, j, &point);
	      gsl_matrix_free(logTransitions);
	    }
	    else{
	      tempProb =  viterbiMaxPrevious(hmm, deltas, hmm->transition_matrix, i, j, &point);
	    }
	    gsl_matrix_set(deltas, i, j, tempProb + gsl_matrix_get(emisLogs, i, j));
	    gsl_matrix_set(psi, i ,j , point);
	  }
	  if (transitionPowerFlag){
	      last = gsl_vector_get(hmm->otherData,j);
	  }
	}

	//termination
	col = gsl_matrix_column(deltas, length - 1);
	lastState = (int) gsl_vector_max_index(&col.vector);

	//backtrack
	gsl_vector_set(viterbiPath, length - 1, lastState);
	for (j = length - 2; j >= 0; j--){
		tmp = gsl_matrix_get(psi, lastState, j + 1);
		lastState =  tmp;
		gsl_vector_set(viterbiPath, j , (int) lastState);
		}

	return gsl_vector_max(&col.vector);
	gsl_matrix_free(deltas);
	gsl_matrix_free(psi);
}

/* this returns the sum of the previous states for the forward probabilities */

double forwardSumPrevious(HMM *hmm, gsl_matrix *alphas, gsl_matrix *transition_score_matrix_i, int state, int obsIndex){
	gsl_vector *tempProbs;
	gsl_vector_view prevCol;
	int i;
	double tempProb;

	tempProbs = gsl_vector_alloc(hmm->nstates);
	prevCol = gsl_matrix_column(alphas, obsIndex - 1);
	for(i = 0; i < hmm->nstates; i++){
		tempProb = gsl_vector_get(&prevCol.vector, i) + gsl_matrix_get(transition_score_matrix_i, i, state);
		gsl_vector_set(tempProbs, i, tempProb);
		}
	tempProb = log_sum(tempProbs);
	gsl_vector_free(tempProbs);
	return(tempProb);
}

/* this returns the sum of the "next" states for the backward probabilities */

double backwardSumNext(HMM *hmm, gsl_matrix *betas, gsl_matrix *transition_score_matrix_j,gsl_matrix *emisLogs, int state, int obsIndex){
	gsl_vector *tempProbs;
	gsl_vector_view nextCol;
	int j;
	double tempProb;

	tempProbs = gsl_vector_alloc(hmm->nstates);
	nextCol = gsl_matrix_column(betas, obsIndex + 1);
	for(j = 0; j < hmm->nstates; j++){		//sum over begin states...
		tempProb = gsl_vector_get(&nextCol.vector, j) + gsl_matrix_get(transition_score_matrix_j, state, j) + gsl_matrix_get(emisLogs, j,obsIndex + 1) ;
		gsl_vector_set(tempProbs, j, tempProb);
		}
	tempProb = log_sum(tempProbs);
	gsl_vector_free(tempProbs);
	return(tempProb);
}

/* this returns the max of the previous states for the viterbi algorithm.
it also sets the back pointer (point) to that state. */

double viterbiMaxPrevious(HMM *hmm, gsl_matrix *deltas, gsl_matrix *transition_score_matrix_i, int state, int obsIndex, int *point){
	gsl_vector *tempProbs;
	gsl_vector_view prevCol;
	int i;
	double tempProb;

	tempProbs = gsl_vector_alloc(hmm->nstates);
	prevCol = gsl_matrix_column(deltas, obsIndex - 1);
	for(i = 0; i < hmm->nstates; i++){
		tempProb = gsl_vector_get(&prevCol.vector, i) + gsl_matrix_get(transition_score_matrix_i, i, state);
		gsl_vector_set(tempProbs, i, tempProb);
		}
	*point = gsl_vector_max_index(tempProbs);
	tempProb = gsl_vector_max(tempProbs);
	gsl_vector_free(tempProbs);
	return(tempProb);
}

/* This computes the log likelihood of an subpath of the data, for a specified set of the states defined
    in the vector stateInclusion (1 for include, 0 for exclude). */
double hmm_subpath_score(HMM *hmm, gsl_matrix *emisLogs, gsl_vector *stateInclusion, int begin, int length, int transitionPowerFlag){
  gsl_matrix *alphas, *real_trans ;
  gsl_vector *real_begin;
  int i, j, stateCount;
  double score;

  //initialize some matrices
  alphas = gsl_matrix_alloc(hmm->nstates, length);
  real_trans = gsl_matrix_alloc(hmm->nstates, hmm->nstates);
  real_begin = gsl_vector_alloc(hmm->nstates);

  //count number of states under consideration
  stateCount = 0;
  for(i = 0; i < hmm->nstates; i++){
    stateCount += gsl_vector_get(stateInclusion, i);
  }

  //keep track of original transition matrix and "begin" vector
  gsl_matrix_memcpy(real_trans,hmm->transition_matrix);
  gsl_vector_memcpy(real_begin, hmm->piStart);

  //set begin states to uniform distribution over states considered in subset
  for(i = 0; i < hmm->nstates; i++){
    if (gsl_vector_get(stateInclusion,i)){
      gsl_vector_set(hmm->piStart, i, 1.0 / stateCount);
    }
    else{
      gsl_vector_set(hmm->piStart, i, 0.0);
    }
  }

  //renomalize transition matrix to reflect states included.
  //go through transition matrix, shutting off those which aren't in subset
  for(i = 0; i < hmm->nstates; i++){
    for(j = 0; j < hmm->nstates; j++){
      if (gsl_vector_get(stateInclusion,i) == 0 || gsl_vector_get(stateInclusion,j) == 0){
	gsl_matrix_set( hmm->transition_matrix, i, j, 0.0);
      }
    }
  }

  //one is the magic number...
  hmm_normalize_transitions(hmm);
  hmm_logify_transitions(hmm);

  //run the forward algorithm to get our answer!
  score = forwardAlgOffset(hmm, alphas, emisLogs, length, begin, transitionPowerFlag);

  //cleanup
  gsl_matrix_memcpy(hmm->transition_matrix, real_trans);
  gsl_vector_memcpy(hmm->piStart, real_begin);
  hmm_logify_transitions(hmm);
  gsl_matrix_free(alphas);
  gsl_matrix_free(real_trans);
  gsl_vector_free(real_begin);
  return(score);
}

/* Calculates the log odds score for two competing sets of states, over a subset of the data
   by comparing the ratio of their likelihoods */
double hmm_subpath_log_odds(HMM *hmm, gsl_matrix *emisLogs, gsl_vector *stateInclusion1,
			     gsl_vector *stateInclusion2,int begin, int length,
			    int transitionPowerFlag){
  double score1, score2;
  score1 = hmm_subpath_score(hmm, emisLogs, stateInclusion1, begin,length,transitionPowerFlag);
  score2 =  hmm_subpath_score(hmm, emisLogs, stateInclusion2, begin,length,transitionPowerFlag);
  //printf("scores: %f %f\n",score1,score2);
  return(score1 - score2);
}

/* simulatePath-- returns a vector snpNumber long that represents an instance
of the simulated markov state path through the data. takes as arguments,an hmm,
 and a pointer to void which actually is an gsl random numb. gen, as well as the observation
 number */
gsl_vector *simulatePath(HMM *h, void *r, int obsNumber){
  int i, j, pastState;
  double sum, rand;
  gsl_rng * rn = (gsl_rng *) r;
  gsl_vector *path;

  path = gsl_vector_alloc(obsNumber);

  /*choose initial state, set paths[0] */
  rand = gsl_rng_uniform(rn);
  sum = 0;
  for(i = 0; i < h->nstates; i++){
    sum += gsl_vector_get(h->piStart, i);
    printf("%f %f\n",rand, sum);
    if (rand <= sum){
      printf("true\n");
      gsl_vector_set(path,0,i);
      rand = 2.0;
    }
  }

  /* now go through obs, using transition matrix to choose states */

  for(i = 1; i < obsNumber; i++){
    pastState = gsl_vector_get(path, i-1);
    rand = gsl_rng_uniform(rn);
    sum = 0;
    for(j = 0; j < h->nstates; j++){
      sum += gsl_matrix_get(h->transition_matrix, pastState, j);
      if (rand <= sum){
	gsl_vector_set(path,i,j);
	rand = 2.0;
      }
    }
  }
  return(path);
}

/* simulatePathSpaced-- just like simulatePath except takes two
added parameters, a pointer to a gsl_vector which will contain the
simulated snp locations, and a double for the exponential mean */
gsl_vector *simulatePathSpaced(HMM *h, void *r, int obsNumber, gsl_vector *locs, double mu){
  int i, j, pastState;
  double sum, rand;
  gsl_rng * rn = (gsl_rng *) r;
  gsl_vector *path;
  gsl_matrix *matPower;

  path = gsl_vector_alloc(obsNumber);

  /*choose initial state, set paths[0] */
  rand = gsl_rng_uniform(rn);
  sum = 0;
  for(i = 0; i < h->nstates; i++){
    sum += gsl_vector_get(h->piStart, i);
    printf("%f %f\n",rand, sum);
    if (rand <= sum){
      printf("true\n");
      gsl_vector_set(path,0,i);
      rand = 2.0;
    }
  }

  /*set initial location to 1 by convention */
  gsl_vector_set(locs,0,1);
  //now go through obs and simulate spacing -- curently drawn from exp
  for(i = 1; i < obsNumber; i++){
    rand = gsl_ran_exponential(rn, mu);
    gsl_vector_set(locs, i, gsl_vector_get(locs, i - 1) + ceil(rand));
  }

  /* now go through obs, using transition matrix and powers to choose states */
  for(i = 1; i < obsNumber; i++){
    pastState = gsl_vector_get(path, i-1);
    matPower = gsl_matrix_power(h->transition_matrix, gsl_vector_get(locs,i) - gsl_vector_get(locs,i - 1));
    rand = gsl_rng_uniform(rn);
    sum = 0;
    for(j = 0; j < h->nstates; j++){
      sum += gsl_matrix_get(matPower, pastState, j);
      if (rand <= sum){
	gsl_vector_set(path,i,j);
	rand = 2.0;
      }
    }
    gsl_matrix_free(matPower);
  }
  return(path);
}

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
#ifndef GMM_MODEL_H
#define GMM_MODEL_H

#include "DoubleVector.h"
#include "GMM/GmmAggregateNewComp.h"
#include "Handle.h"
#include "LambdaCreationFunctions.h"
#include "Object.h"
#include "PDBVector.h"

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_vector.h>
#include <math.h>
#include <vector>

// By Tania, October 2017

using namespace pdb;

/* This class contains our current GMM model, that is:
 * - num of dimensions;
 * - weights: set of (prior) weights for each component
 * - means: mean for each component
 * - covars: covariance for each component
 * - inv_covars: inverse covariance for each component
 * - dets: determinant of each covariance matrix
 */
class GmmModel : public Object {

private:
  int ndim = 0;
  int k = 0;
  Handle<DoubleVector> weights;
  Handle<Vector<Handle<DoubleVector>>> means;      // k*ndim
  Handle<Vector<Handle<DoubleVector>>> covars;     // k*ndim*ndim
  Handle<Vector<Handle<DoubleVector>>> inv_covars; // k*ndim*ndim
  Handle<DoubleVector> dets;

public:
  ENABLE_DEEP_COPY

  GmmModel() {}

  GmmModel(int k, int ndim) {
    std::cout << "K= " << k << std::endl;
    std::cout << "ndim= " << ndim << std::endl;

    this->ndim = ndim;
    this->k = k;

    this->weights = pdb::makeObject<DoubleVector>(k);
    this->dets = pdb::makeObject<DoubleVector>(k);

    this->means = pdb::makeObject<Vector<Handle<DoubleVector>>>();
    this->covars = pdb::makeObject<Vector<Handle<DoubleVector>>>();
    this->inv_covars = pdb::makeObject<Vector<Handle<DoubleVector>>>();

    for (int i = 0; i < k; i++) {
      this->weights->setDouble(i, 1.0 / k);
      this->means->push_back(makeObject<DoubleVector>(ndim));
      this->covars->push_back(makeObject<DoubleVector>(ndim * ndim));
      this->inv_covars->push_back(makeObject<DoubleVector>(ndim * ndim));
    }

    std::cout << "Model created!" << std::endl;
  }

  // this prints the first 10 elements of a DoubleVector
  void print_vector(Handle<DoubleVector> v) {
    int MAX = 10;
    if (v->size < MAX)
      MAX = v->size;
    for (int i = 0; i < MAX; i++) {
      std::cout << i << ": " << v->getDouble(i) << "; ";
    }
    std::cout << std::endl;
  }

  // Print current model
  void print() {
    // Print mean
    std::cout << "**Current model**" << std::endl;
    std::cout << "Weights: ";
    print_vector(this->weights);
    for (int i = 0; i < k; i++) {
      std::cout << "**Component " << i << std::endl;
      std::cout << "Mean: ";
      print_vector((*this->means)[i]);
      std::cout << "Covar: ";
      print_vector((*this->covars)[i]);
      std::cout << "Inv covars: ";
      print_vector((*this->inv_covars)[i]);
    }
  }

  int getNumK() { return k; }

  int getNDim() { return ndim; }

  double getWeight(int i) { return this->weights->getDouble(i); }

  void updateWeights(std::vector<double> w) {
    for (int i = 0; i < k; i++) {
      this->weights->setDouble(i, w[i]);
    }
    std::cout << "updated weights";
    print_vector(this->weights);
  }

  void updateMeans(std::vector<std::vector<double>> &m) {
    for (int i = 0; i < k; i++) {
      std::cout << "updated mean for comp " << i << std::endl;

      double *mean = (*this->means)[i]->getRawData();
      for (int j = 0; j < ndim; j++) {
        mean[j] = m[i][j];
      }
      print_vector((*this->means)[i]);
    }
  }

  void updateCovars(std::vector<std::vector<double>> &c) {
    for (int i = 0; i < k; i++) {
      std::cout << "updated covar for comp " << i << std::endl;

      double *covar = (*this->covars)[i]->getRawData();

      for (int j = 0; j < (ndim * ndim); j++) {
        covar[j] = c[i][j];
      }
      print_vector((*this->covars)[i]);
    }

    // calcInvCovars();
  }

  // It calculates determinants and inverse covariances based on Cholesky
  // decomposition
  void calcInvCovars() {

    for (int i = 0; i < this->getNumK(); i++) {

      gsl_matrix_view covar =
          gsl_matrix_view_array((*covars)[i]->data->c_ptr(), ndim, ndim);
      gsl_matrix_view inv_covar =
          gsl_matrix_view_array((*inv_covars)[i]->data->c_ptr(), ndim, ndim);

      gsl_matrix_memcpy(&inv_covar.matrix, &covar.matrix);

      // gsl_permutation *p = gsl_permutation_alloc(ndim);

      gsl_linalg_cholesky_decomp(&inv_covar.matrix);

      // Calculate determinant -> logdet
      double ax = 0.0;
      for (int i = 0; i < ndim; i++) {
        ax += log(gsl_matrix_get(&inv_covar.matrix, i, i));
      }

      // Check for
      gsl_linalg_cholesky_invert(&inv_covar.matrix);

      dets->setDouble(i, ax);
    }
  }

  // It returns the rvalue, responsability, posterior probability of the
  // datapoint
  // for GMM component i
  double log_normpdf(int i, Handle<DoubleVector> inputData) {

    gsl_vector_view data =
        gsl_vector_view_array(inputData->data->c_ptr(), ndim);
    gsl_vector_view mean =
        gsl_vector_view_array((*means)[i]->data->c_ptr(), ndim);
    gsl_matrix_view covar =
        gsl_matrix_view_array((*covars)[i]->data->c_ptr(), ndim, ndim);
    gsl_matrix_view inv_covar =
        gsl_matrix_view_array((*inv_covars)[i]->data->c_ptr(), ndim, ndim);

    double ax, ay;
    int s;
    gsl_vector *ym;
    gsl_vector *datan = gsl_vector_alloc(ndim);

    ax = dets->getDouble(i);

    gsl_vector_memcpy(datan, &data.vector);
    gsl_vector_sub(datan, &mean.vector);

    ym = gsl_vector_alloc(ndim);
    gsl_blas_dsymv(CblasUpper, 1.0, &inv_covar.matrix, datan, 0.0, ym);
    gsl_blas_ddot(datan, ym, &ay);

    gsl_vector_free(ym);
    gsl_vector_free(datan);

    ay = -ax - 0.5 * ay;

    return ay;
  }

  // It calculates the sum of elements in vector v in log space
  // To do this, we extract the maximum factor (maxVal) from
  // all the values
  double logSumExp(double *v, size_t size) {

    double maxVal = v[0];
    double sum = 0;

    for (int i = 1; i < size; i++) {
      if (v[i] > maxVal) {
        maxVal = v[i];
      }
    }

    for (int i = 0; i < size; i++) {
      sum += exp(v[i] - maxVal);
    }

    return log(sum) + maxVal;
  }

  double logSumExp(DoubleVector &vector) {
    return logSumExp(vector.getRawData(), vector.size);
  }

  double logSumExp(Vector<double> &vector) {
    return logSumExp(vector.c_ptr(), vector.size());
  }

  // It uses the partial results from the Aggregate to update the model

  double updateModel(GmmAggregateNewComp update) {

    const double MIN_COVAR = 1e-6;

    double *sumWeights = update.getSumWeights().c_ptr();
    double totalSumR = 0.0;

    for (int i = 0; i < k; i++) {
      totalSumR += sumWeights[i];
    }

    // Update weights -> weights = sumR/totalSumR
    gsl_vector_view gSumWeights = gsl_vector_view_array(sumWeights, k);
    gsl_vector_view gweights =
        gsl_vector_view_array((*weights).data->c_ptr(), k);
    gsl_vector_memcpy(&gweights.vector, &gSumWeights.vector);
    gsl_vector_scale(&gweights.vector, 1.0 / totalSumR);

    // Update Mean and Covar for each component
    for (int i = 0; i < k; i++) {

      // Update means -> means = 1/sumR * weightedX
      gsl_vector_view mean =
          gsl_vector_view_array((*means)[i]->data->c_ptr(), ndim);
      gsl_vector_view gsumMean =
          gsl_vector_view_array(update.getSumMean(i).c_ptr(), ndim);
      gsl_vector_memcpy(&mean.vector, &gsumMean.vector);
      gsl_vector_scale(&mean.vector, 1.0 / sumWeights[i]);

      // Update covars -> gsumCovarv / sumR - mean*mean^T

      gsl_vector_view gsumCovarv =
          gsl_vector_view_array(update.getSumCovar(i).c_ptr(), ndim * ndim);
      gsl_vector_scale(&gsumCovarv.vector, 1.0 / sumWeights[i]);
      gsl_matrix_view gsumCovarm =
          gsl_matrix_view_array(update.getSumCovar(i).c_ptr(), ndim, ndim);
      gsl_blas_dsyr(CblasUpper, -1.0, &mean.vector, &gsumCovarm.matrix);

      // Add constant to diagonal and copy lower triangular
      double d;
      for (int row = 0; row < ndim; row++) {
        d = gsl_matrix_get(&gsumCovarm.matrix, row, row);
        gsl_matrix_set(&gsumCovarm.matrix, row, row, d + MIN_COVAR);

        for (int col = row + 1; col < ndim; col++) {
          // matrix[j][i] = matrix[i][j]
          d = gsl_matrix_get(&gsumCovarm.matrix, row, col);
          gsl_matrix_set(&gsumCovarm.matrix, col, row, d);
        }
      }

      gsl_vector_view covar =
          gsl_vector_view_array((*covars)[i]->data->c_ptr(), ndim * ndim);
      gsl_vector_memcpy(&covar.vector, &gsumCovarv.vector);
    }

    // Finally, calculate inverse covariances and determinants to be used in the
    // next iteration
    calcInvCovars();

    return update.getLogLikelihood();
  }

  ~GmmModel() {}
};

#endif


text:
#pragma once

/*** mathematics.h:

     Mathematical and statistical functions. ***/

/** Avoid including this file twice. **/

#ifndef MATHEMATICS_H
#define MATHEMATICS_H

/** Libraries. **/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/** Constants. **/

#define MATH_PI 3.1415926535897932385 /* Pi. */
#define MATH_2PI 6.2831853071795864770 /* 2 pi. */
#define MATH_PI_2 1.57079632679489661925 /* Pi / 2. */

/* Coordinates / axes. */

#define MATH_X 0 /* Castesian. */
#define MATH_Y 1
#define MATH_Z 2
#define MATH_RADIUS 0 /* Spherical. */
#define MATH_POLAR 1
#define MATH_AZIMUTH 2
#define MATH_MU 1 /* Pseudospherical (mu = cos(polar) replacing polar). */

/** Variables. **/

/* Random numbers. */

gsl_rng *RandomNumberGenerator; /* Random number generator. */
unsigned RandomSeed; /* Random seed. */
unsigned SortComponent = 0; /* Component used for sorting vectors. */

/** Function prototypes. **/

int SortAscending(const void *, const void *);
int SortIndexesAscending(const void *, const void *);
double ScalarProduct(double *, double *, unsigned);
double Norm(double *, unsigned);
double NormSquared(double *, unsigned);
double AngleCosine(double *, double *);

void SphericalToCartesian(double *, double *);
void PseudoSphericalToCartesian(double *, double *);
void CartesianToSpherical(double *, double *);
void CartesianToPseudoSpherical(double *, double *);
void Rotate(double *, double *, unsigned, double);
unsigned GetRandomSeed();
int InitialiseRandom(const gsl_rng_type *);
int FinaliseRandom();
double MetropolisHastings(double *, double *, double *, double *, unsigned, unsigned, unsigned, unsigned, double *, double *, double *,
                          void (*)(double *, double *, double *, double *, double *, unsigned), double (*)(double *, double *, double *, double *, double *, unsigned),
                          double (*)(double *));
double ParallelTempering(double *, double *, double *, double *, unsigned, unsigned, unsigned, unsigned, double *, double *, double *,
                         void (*)(double *, double *, double *, double *, double *, unsigned), double (*)(double *, double *, double *, double *, double *, unsigned),
                         double (*)(double *), unsigned, double *, double);
void TanWeights(double *, double *, double *, double *, unsigned, unsigned, unsigned, double *, double *, double *,
                double (*)(double *, double *, double *, double *, double *, unsigned), unsigned);
double TanNormalisation(double *, unsigned, double *);
double WeightAndSum(double *, double *, unsigned);
void UniformTransition(double *, double *, double *, double *, double *, unsigned);
double UniformTransitionPDF(double *, double *, double *, double *, double *, unsigned);
double PowerLawDeviate(double, double, double);
void IsotropicDirection(double, double *);
void StandardGaussianDeviates(double *);
void GaussianDeviates(double *, double, double);
complex ComplexGaussianDeviates(complex, complex);

/** SortAscending.

    - Function designed to be an argument of qsort().
    - Sorts the elements of a vector in ascending order. The vector elements are reorganised. **/

int SortAscending(const void *x, const void *y)
{
  double difference;

  difference = *((double *)x) - *((double *)y);
  if (difference > 0.0)
    return 1;
  else if (difference < 0.0)
    return -1;
  else
    return 0;
}

/** SortIndexesAscending.

    - Function designed to be an argument of qsort().
    - Sorts the indexes of a vector, in such a way that its elements are arranged in ascending order. The vector remains unchanged.
    - 'SortComponent' indicates the component used to sort, in the case that the vector elements have more than one dimension. **/

int SortIndexesAscending(const void *x, const void *y)
{
  double difference;

  difference = *(*(double **)x + SortComponent) - *(*(double **)y + SortComponent);
  if (difference > 0.0)
    return 1;
  else if (difference < 0.0)
    return -1;
  else
    return 0;
}

/** ScalarProduct.

    - Computes the scalar product between two vectors. **/

double ScalarProduct(double *vector_1, double *vector_2, unsigned vector_size)
{
  unsigned i;
  double product;

  /* Compute components product and sum. */

  product = 0.0;
  for (i = 0; i < vector_size; i++)
    product += vector_1[i] * vector_2[i];
  return product;
}

/** Norm.

    - Computes the norm of a vector. **/

double Norm(double *vector, unsigned vector_size)
{
  unsigned i;
  double *component_sq, max_component_sq, norm;

  /* Generate an auxiliary vector of component squares, in ascending order. */

  component_sq = malloc(vector_size * sizeof(double));
  for (i = 0; i < vector_size; i++)
    component_sq[i] = vector[i] * vector[i];
  qsort(component_sq, vector_size, sizeof(double), SortAscending);

  /* Divide component squares by the largest one to avoid overflow. */

  max_component_sq = component_sq[vector_size-1];
  if (max_component_sq <= 0.0) /* All components are null. */
    return 0.0;
  vector_size--;
  for (i = 0; i < vector_size; i++)
    component_sq[i] /= max_component_sq;

  /* Sum component squares in ascending order, to avoid loss of significance. */

  norm = 0.0;
  for (i = 0; i < vector_size; i++)
    norm += component_sq[i];
  norm += 1.0;

  /* Renormalise. */

  norm = sqrt(max_component_sq * norm);
  free(component_sq);
  return norm;
}

/** NormSquared.

    - Computes the square of the norm of a vector. **/

double NormSquared(double *vector, unsigned vector_size)
{
  unsigned i;
  double *component_sq, max_component_sq, norm;

  /* Generate an auxiliary vector of component squares, in ascending order. */

  component_sq = malloc(vector_size * sizeof(double));
  for (i = 0; i < vector_size; i++)
    component_sq[i] = vector[i] * vector[i];
  qsort(component_sq, vector_size, sizeof(double), SortAscending);

  /* Divide component squares by the largest one to avoid overflow. */

  max_component_sq = component_sq[vector_size-1];
  if (max_component_sq <= 0.0) /* All components are null. */
    return 0.0;
  vector_size--;
  for (i = 0; i < vector_size; i++)
    component_sq[i] /= max_component_sq;

  /* Sum component squares in ascending order, to avoid loss of significance. */

  norm = 0.0;
  for (i = 0; i < vector_size; i++)
    norm += component_sq[i];
  norm += 1.0;

  /* Renormalise. */

  norm *= max_component_sq;
  free(component_sq);
  return norm;
}

/** AngleCosine.

    - Computes the cosine of the angle between two 3D vectors. **/

double AngleCosine(double *vector_1, double *vector_2)
{
  double norm_1, norm_2, cosine;

  norm_1 = Norm(vector_1, 3);
  norm_2 = Norm(vector_2, 3);
  if (norm_1 <= 0.0 || norm_2 <= 0.0) /* At least one vector is null. */
    return 0.0;
  cosine = ScalarProduct(vector_1, vector_2, 3) / (norm_1 * norm_2);
  if (cosine > 1.0) /* Correct eventual round-off errors. */
    cosine = 1.0;
  else if (cosine < -1.0)
    cosine = -1.0;
  return cosine;
}

/** SphericalToCartesian.

    - Converts a 3D vector from spherical (radius, polar angle, azimuth) to Cartesian form.
    - Initial and final vectors may be the same. **/

void SphericalToCartesian(double *vec_in, double *vec_out)
{
  double sin_polar, cos_polar, radius, azimuth;

  radius = vec_in[MATH_RADIUS];
  sin_polar = sin(vec_in[MATH_POLAR]);
  cos_polar = cos(vec_in[MATH_POLAR]);
  azimuth = vec_in[MATH_AZIMUTH];
  vec_out[MATH_X] = radius * sin_polar * cos(azimuth);
  vec_out[MATH_Y] = radius * sin_polar * sin(azimuth);
  vec_out[MATH_Z] = radius * cos_polar;
}

/** PseudoSphericalToCartesian.

    - Converts a 3D vector from pseudo spherical (radius, polar angle cosine, azimuth) to Cartesian form.
    - Initial and final vectors may be the same. **/

void PseudoSphericalToCartesian(double *vec_in, double *vec_out)
{
  double cos_polar, sin_polar, radius, azimuth;

  radius = vec_in[MATH_RADIUS];
  azimuth = vec_in[MATH_AZIMUTH];
  cos_polar = vec_in[MATH_MU];
  sin_polar = 1.0 - cos_polar * cos_polar;
  if (sin_polar <= 0.0) /* Correct eventual round-off errors. */
    sin_polar = 0.0;
  else
    sin_polar = sqrt(sin_polar);
  vec_out[MATH_X] = radius * sin_polar * cos(azimuth);
  vec_out[MATH_Y] = radius * sin_polar * sin(azimuth);
  vec_out[MATH_Z] = radius * cos_polar;
}

/** CartesianToSpherical.

    - Converts a 3D vector from Cartesian to spherical (radius, polar angle, azimuth) form.
    - Initial and final vectors may be the same. **/

void CartesianToSpherical(double *vec_in, double *vec_out)
{
  double radius, azimuth, polar;

  radius = Norm(vec_in, 3);
  if (radius > 0.0)
    polar = acos(vec_in[MATH_Z] / radius);
  else
    polar = 0.0;
  azimuth = atan2(vec_in[MATH_Y], vec_in[MATH_X]);
  vec_out[MATH_RADIUS] = radius;
  vec_out[MATH_POLAR] = polar;
  vec_out[MATH_AZIMUTH] = azimuth;
}

/** CartesianToPseudoSpherical.

    - Converts a 3D vector from Cartesian to pseudo spherical (radius, polar angle cosine, azimuth) form.
    - Initial and final vectors may be the same. **/

void CartesianToPseudoSpherical(double *vec_in, double *vec_out)
{
  double radius, cos_polar, azimuth;

  radius = Norm(vec_in, 3);
  if (radius > 0.0)
    cos_polar = vec_in[MATH_Z] / radius;
  else
    cos_polar = 1.0;
  azimuth = atan2(vec_in[MATH_Y], vec_in[MATH_X]);
  vec_out[MATH_RADIUS] = radius;
  vec_out[MATH_MU] = cos_polar;
  vec_out[MATH_AZIMUTH] = azimuth;
}

/** Rotate.

    - Rotates a vector.
    - Vector is rotated counterclockwise (or axes clockwise) as seen from the tip of the rotation axis specified by argument 'axis'.
    - Initial and final vectors may be the same. **/

void Rotate(double *vec_in, double *vec_out, unsigned axis, double angle)
{
  unsigned i, j;
  double cos_angle, sin_angle, x, y;

  cos_angle = cos(angle);
  sin_angle = sin(angle);
  axis %= 3;
  vec_out[axis] = vec_in[axis];
  i = (axis + 1) % 3;
  j = (axis + 2) % 3;
  x = vec_in[i];
  y = vec_in[j];
  vec_out[i] = x * cos_angle - y * sin_angle;
  vec_out[j] = x * sin_angle + y * cos_angle;
}

/** GetRandomSeed.

    - Provides a random seed by reading garbage from /dev/urandom. **/

unsigned GetRandomSeed()
{
  FILE *random_file;
  unsigned seed;
  int error;

  seed = 0;
  random_file = fopen("/dev/urandom", "rb");
  if (!random_file)
    return 0;
  do
  {
    error = fread(&seed, 1, sizeof(unsigned int), random_file);
    if (error != sizeof(unsigned int))
      return 0;
  }
  while (!seed);
  fclose(random_file);
  return seed;
}

/** InitialiseRandom.

    - Initialises the random number generator. **/

int InitialiseRandom(const gsl_rng_type *random_generator)
{
  RandomSeed = GetRandomSeed();
  if (!RandomSeed)
    return EXIT_FAILURE;
  RandomNumberGenerator = gsl_rng_alloc(random_generator);
  if (!RandomNumberGenerator)
    return EXIT_FAILURE;
  gsl_rng_set(RandomNumberGenerator, RandomSeed);
  return EXIT_SUCCESS;
}

/** FinaliseRandom.

    - Finalises the random number generator. **/

int FinaliseRandom()
{
  /* Release random number generator memory. */

  if (RandomNumberGenerator)
  {
    gsl_rng_free(RandomNumberGenerator);
    return EXIT_SUCCESS;
  }
  return EXIT_FAILURE;
}

/** MetropolisHastings.

    - Constructs a Markov chain using Metropolis-Hastings sampling from a given PDF ('sampling_pdf').
    - 'states' contains the actual Markov chain, comprising 'chain_size' elements. Each element contains 'state_size' doubles. The initial state must be stored in
      the first element of 'states'. 'states' must be allocated by the calling function.
    - 'proposals' is the chain of proposals derived from each state (same length and element size as 'states'). It must be allocated by the calling function.
    - 'lower_bound' and 'upper_bound' must contain the absolute bounds for the variables defining the state.
    - 'step' must contain the parameter defining the size of the transition step.
    - 'sampling_pdf' must point to a function computing the PDF sampled by the Markov chain. Its only argument is a state of the chain; any other parameter must be
      passed through global variables.
    - 'transition_pdf' must point to a function computing the PDF of a Metropolis-Hastings transition, given a Markov chain state. Its arguments are the final and
      initial states of the transition, in this order. Any other parameter must be passed through global variables.
    - 'transition' must point to a function computing the proposed final state of a Metropolis-Hastings transition, given a Markov chain state. Its arguments are the
      final and initial states of the transition, in this order. Any other parameter must be passed through global variables.
    - 'proposals_PDF' is an array of doubles of the same length of the Markov chain. It must be allocated by the calling function if storage of sampling PDF values
      for the proposals is desired, or set to NULL otherwise.
    - 'states_PDF' is an array of doubles of the same length of the Markov chain. It must be allocated by the calling function if storage of sampling PDF values
      for the states is desired, or set to NULL otherwise.
    - If the initial state set in 'states' has a null value of the sampling PDF, the first transition probability is set to unity. Once a state with a non-null value
      of the sampling PDF is reached, the sampling PDF values do not vanish over the rest of the Markov chain by construction.
    - The order of the state and proposal chains is the following: any element of 'proposals' is sampled from the corresponding element of 'states', whereas any
      element of 'states' (except the first one) is the previous element of either 'states' or 'proposals'.
    - If 'save_states' is set to null, only the final state and proposal are given. Otherwise, the full chains are saved.
    - If 'log_flag' is non null, 'sampling_pdf' is assumed to give the natural logarithm of the PDF.
    - Returns the acceptance ratio of the chain. **/

double MetropolisHastings(double *states, double *proposals, double *proposals_pdf, double *states_pdf, unsigned state_size, unsigned chain_size, unsigned save_states,
                          unsigned log_flag, double *lower_bound, double *upper_bound, double *step,
			  void (*transition)(double *, double *, double *, double *, double *, unsigned),
                          double (*transition_pdf)(double *, double *, double *, double *, double *, unsigned), double (*sampling_pdf)(double *))
{
  unsigned i, accepted_transitions, state_size_bytes;
  double direct_transition_pdf, inverse_transition_pdf, pdf_ratio, state_pdf, proposal_pdf, *proposal, *state, *old_state, acceptance_ratio;

  /* Set initial stuff. */

  state_size_bytes = state_size * sizeof(double);
  proposal = proposals;
  state = states;
  state_pdf = sampling_pdf(state);
  if (states_pdf)
    states_pdf[0] = state_pdf;

  /* Construct the Markov chain. */

  accepted_transitions = 0;
  chain_size--;
  for (i = 0; i < chain_size; i++)
  {
    /* Set new proposal. */

    transition(proposal, state, lower_bound, upper_bound, step, state_size);
    proposal_pdf = sampling_pdf(proposal);
    if (proposals_pdf && save_states)
      proposals_pdf[i] = proposal_pdf;

    /* Compute transition probabilities. */

    direct_transition_pdf = transition_pdf(proposal, state, lower_bound, upper_bound, step, state_size);
    inverse_transition_pdf = transition_pdf(state, proposal, lower_bound, upper_bound, step, state_size);
    if (!log_flag && state_pdf <= 0.0)
      pdf_ratio = 1.1;
    else if (direct_transition_pdf <= 0.0)
      pdf_ratio = -0.1;
    else if (log_flag)
      pdf_ratio = exp(proposal_pdf - state_pdf) * inverse_transition_pdf / direct_transition_pdf;
    else
      pdf_ratio = proposal_pdf * inverse_transition_pdf / (state_pdf * direct_transition_pdf);

    /* Accept of reject transition. */

    if (save_states)
    {
      old_state = state;
      state += state_size;
    }
    if (gsl_rng_uniform(RandomNumberGenerator) < pdf_ratio)
    {
      memcpy(state, proposal, state_size_bytes);
      state_pdf = proposal_pdf;
      accepted_transitions++;
    }
    else if (save_states)
      memcpy(state, old_state, state_size_bytes);
    if (save_states)
    {
      proposal += state_size;
      if (states_pdf)
        states_pdf[i+1] = state_pdf;
    }
  }

  /* Compute final stuff. */

  transition(proposal, state, lower_bound, upper_bound, step, state_size);
  acceptance_ratio = accepted_transitions / (double) chain_size;
  if (proposals_pdf && save_states)
    proposals_pdf[chain_size] = sampling_pdf(proposal);
  return acceptance_ratio;
}

/** ParallelTempering.

    - Constructs a Markov chain using Metropolis-Hastings sampling from a given PDF ('sampling_pdf'), with parallel tempering.
    - 'states' contains the cold Markov chain, comprising 'chain_size' elements. Each element contains 'state_size' doubles. The initial state must be stored in
      the first element of 'states'. 'states' must be allocated by the calling function.
    - 'proposals' is the chain of proposals derived from each state (same length and element size as 'states'). It must be allocated by the calling function.
    - 'lower_bound' and 'upper_bound' must contain the absolute bounds for the variables defining the state.
    - 'step' must contain the parameter defining the size of the transition step.
    - 'sampling_pdf' must point to a function computing the PDF sampled by the cold Markov chain. Its only argument is a state of the chain; any other parameter must
      be passed through global variables.
    - 'transition_pdf' must point to a function computing the PDF of a Metropolis-Hastings transition, given a Markov chain state. Its arguments are the final and
      initial states of the transition, in this order. Any other parameter must be passed through global variables.
    - 'transition' must point to a function computing the proposed final state of a Metropolis-Hastings transition, given a Markov chain state. Its arguments are the
      final and initial states of the transition, in this order. Any other parameter must be passed through global variables.
    - 'proposals_PDF' is an array of doubles of the same length of the Markov chain. It must be allocated by the calling function if storage of sampling PDF values
      for the proposals is desired, or set to NULL otherwise.
    - 'states_PDF' is an array of doubles of the same length of the Markov chain. It must be allocated by the calling function if storage of sampling PDF values
      for the states is desired, or set to NULL otherwise.
    - If the initial state set in 'states' has a null value of the sampling PDF, the first transition probability is set to unity. Once a state with a non-null value
      of the sampling PDF is reached, the sampling PDF values do not vanish over the rest of the Markov chain by construction.
    - The order of the state and proposal chains is the following: any element of 'proposals' is sampled from the corresponding element of 'states', whereas any
      element of 'states' (except the first one) is the previous element of either 'states' or 'proposals'.
    - If 'save_states' is set to null, only the final state and proposal are given. Otherwise, the full chains are saved.
    - If 'log_flag' is non null, 'sampling_pdf' is assumed to give the natural logarithm of the PDF.
    - The number of hot chains in the parallel-tempering scheme is 'hot_chains_count', which have tempering factors given by 'tempering_factors'.
    - 'exchange_probability' gives the probability that two chains are interchanged in a given step.
    - Returns the acceptance ratio of the chain. **/

double ParallelTempering(double *states, double *proposals, double *proposals_pdf, double *states_pdf, unsigned state_size, unsigned chain_size, unsigned save_states,
                         unsigned log_flag, double *lower_bound, double *upper_bound, double *step,
                         void (*transition)(double *, double *, double *, double *, double *, unsigned),
                         double (*transition_pdf)(double *, double *, double *, double *, double *, unsigned), double (*sampling_pdf)(double *),
                         unsigned hot_chains_count, double *tempering_factors, double exchange_probability)
{
  unsigned i, j, accepted_transitions, state_size_bytes;
  double direct_transition_pdf, inverse_transition_pdf, pdf_ratio, state_pdf, proposal_pdf, *proposal, *state, *old_state, acceptance_ratio, x, *aux_state;
  double *cool_state, hot_proposals_pdf, factor, *hot_proposals, *hot_states, *hot_states_pdf, *cool_state_pdf;

  /* Set initial stuff. */

  state_size_bytes = state_size * sizeof(double);
  proposal = proposals;
  state = states;
  state_pdf = sampling_pdf(state);
  if (states_pdf)
    states_pdf[0] = state_pdf;
  hot_proposals = malloc(hot_chains_count * state_size_bytes);
  hot_states = malloc(hot_chains_count * state_size_bytes);
  aux_state = malloc(state_size_bytes);
  for (i = 0; i < hot_chains_count; i++)
    memcpy(hot_states + i * state_size, state, state_size_bytes); /* Copy initial state of cold chain to all hot chains. */
  hot_states_pdf = malloc(hot_chains_count * sizeof(double));
  for (i = 0; i < hot_chains_count; i++)
  {
    if (log_flag)
      hot_states_pdf[i] = state_pdf * tempering_factors[i];
    else
      hot_states_pdf[i] = pow(state_pdf, tempering_factors[i]);
  }

  /* Construct the Markov chains. */

  accepted_transitions = 0;
  chain_size--;
  for (i = 0; i < chain_size; i++)
  {
    /* Set new proposal for cold chain. */

    transition(proposal, state, lower_bound, upper_bound, step, state_size);
    proposal_pdf = sampling_pdf(proposal);
    if (proposals_pdf && save_states)
      proposals_pdf[i] = proposal_pdf;

    /* Compute transition probabilities for cold chain. */

    direct_transition_pdf = transition_pdf(proposal, state, lower_bound, upper_bound, step, state_size);
    inverse_transition_pdf = transition_pdf(state, proposal, lower_bound, upper_bound, step, state_size);
    if (!log_flag && state_pdf <= 0.0)
      pdf_ratio = 1.1;
    else if (direct_transition_pdf <= 0.0)
      pdf_ratio = -0.1;
    else if (log_flag)
      pdf_ratio = exp(proposal_pdf - state_pdf) * inverse_transition_pdf / direct_transition_pdf;
    else
      pdf_ratio = proposal_pdf * inverse_transition_pdf / (state_pdf * direct_transition_pdf);

    /* Accept of reject transition for cold chain. */

    if (save_states)
    {
      old_state = state;
      state += state_size;
    }
    if (gsl_rng_uniform(RandomNumberGenerator) < pdf_ratio)
    {
      memcpy(state, proposal, state_size_bytes);
      state_pdf = proposal_pdf;
      accepted_transitions++;
    }
    else if (save_states)
      memcpy(state, old_state, state_size_bytes);
    if (save_states)
    {
      proposal += state_size;
      if (states_pdf)
        states_pdf[i+1] = state_pdf;
    }

    /* Hot chains. */

    for (j = 0; j < hot_chains_count; j++)
    {
      /* Set new proposal. */

      transition(hot_proposals + j * state_size, hot_states + j * state_size, lower_bound, upper_bound, step, state_size);
      hot_proposals_pdf = sampling_pdf(hot_proposals + j * state_size);
      if (log_flag)
        hot_proposals_pdf *= tempering_factors[j];
      else
        hot_proposals_pdf = pow(hot_proposals_pdf, tempering_factors[j]);

      /* Compute transition probabilities. */

      direct_transition_pdf = transition_pdf(hot_proposals + j * state_size, hot_states + j * state_size, lower_bound, upper_bound, step, state_size);
      inverse_transition_pdf = transition_pdf(hot_states + j * state_size, hot_proposals + j * state_size, lower_bound, upper_bound, step, state_size);
      if (!log_flag && hot_states_pdf[j] <= 0.0)
        pdf_ratio = 1.1;
      else if (direct_transition_pdf <= 0.0)
        pdf_ratio = -0.1;
      else if (log_flag)
        pdf_ratio = exp(hot_proposals_pdf - hot_states_pdf[j]) * inverse_transition_pdf / direct_transition_pdf;
      else
        pdf_ratio = hot_proposals_pdf * inverse_transition_pdf / (hot_states_pdf[j] * direct_transition_pdf);

      /* Accept of reject transition. */

      if (gsl_rng_uniform(RandomNumberGenerator) < pdf_ratio)
      {
        memcpy(hot_states + j * state_size, hot_proposals + j * state_size, state_size_bytes);
        hot_states_pdf[j] = hot_proposals_pdf;
      }
    }

    /* Exchange chain states. */

    if (gsl_rng_uniform(RandomNumberGenerator) < exchange_probability)
    {
      j = gsl_rng_uniform_int(RandomNumberGenerator, hot_chains_count);
      factor = tempering_factors[j];
      if (j)
      {
        factor /= tempering_factors[j-1];
        cool_state_pdf = hot_states_pdf + j - 1;
        cool_state = hot_states + (j - 1) * state_size;
      }
      else
      {
        cool_state_pdf = &state_pdf;
        cool_state = state;
      }
      if (log_flag)
        pdf_ratio = exp(*cool_state_pdf * (factor - 1.0) + hot_states_pdf[j] * (1.0 / factor - 1.0));
      else
        pdf_ratio = pow(*cool_state_pdf, factor - 1.0) * pow(hot_states_pdf[j], 1.0 / factor - 1.0);
      if (gsl_rng_uniform(RandomNumberGenerator) < pdf_ratio)
      {
        if (log_flag)
        {
          x = hot_states_pdf[j];
          hot_states_pdf[j] = *cool_state_pdf * factor;
          *cool_state_pdf = x / factor;
        }
        else
        {
          x = hot_states_pdf[j];
          hot_states_pdf[j] = pow(*cool_state_pdf, factor);
          *cool_state_pdf = pow(x, 1.0 / factor);
        }
        memcpy(aux_state, hot_states + j * state_size, state_size_bytes);
        memcpy(hot_states + j * state_size, cool_state, state_size_bytes);
        memcpy(cool_state, aux_state, state_size_bytes);
      }
    }
  }
  transition(proposal, state, lower_bound, upper_bound, step, state_size);
  acceptance_ratio = accepted_transitions / (double) chain_size;
  if (proposals_pdf && save_states)
    proposals_pdf[chain_size] = sampling_pdf(proposal);
  free(hot_proposals);
  free(hot_states);
  free(hot_states_pdf);
  free(aux_state);
  return acceptance_ratio;
}

/** TanWeights.

    - Computes the weights proposed by Tan (2006, J.Comp.Graph.Stat., 15, 735) to use Metropolis-Hastings sampling for the evaluation of integrals via the Monte
      Carlo method.
    - 'states' contains the actual Markov chain, comprising 'chain_size' elements. Each element contains 'state_size' doubles.
    - 'proposals' contains the chain of proposals derived from each state (same length and element size as 'states').
    - 'proposals_pdf' is an array of doubles of the same length of the Markov chain, containing the sampling PDF values for the proposals.
    - If the present chain is an extension of an old one, 'old_chain_size' is the size of the latter.
    - 'weights' is an array of 'chain_size' elements for storing the resulting weights. It must be allocated by the calling function. Elements of 'weights' between
      'old_chain_size' and 'chain_size' are erased before computing.
    - 'lower_bound' and 'upper_bound' must contain the absolute bounds for the variables defining the state.
    - 'step' must contain the parameter defining the size of the transition step.
    - 'sort_component' indicates the component of the state used for sorting them, to reduce the computational cost if the PDF is bound.
    - 'transition_pdf' must point to a function computing the PDF of a Metropolis-Hastings transition, given a Markov chain state. Its arguments are the final and
      initial states of the transition, in this order. Any other parameter must be passed through global variables. **/

void TanWeights(double *states, double *proposals, double *proposals_pdf, double *weights, unsigned state_size, unsigned chain_size, unsigned old_chain_size,
                double *lower_bound, double *upper_bound, double *step, double (*transition_pdf)(double *, double *, double *, double *, double *, unsigned),
                unsigned sort_component)
{
  unsigned i, j, k, m, n;
  double *proposal, **pointers, pdf;

  if (sort_component >= state_size)
    SortComponent = 0;
  else
    SortComponent = sort_component;
  memset(weights + old_chain_size, 0, (chain_size - old_chain_size) * sizeof(double));
  pointers = malloc(chain_size * sizeof(double *));
  for (i = 0; i < chain_size; i++)
    pointers[i] = states + i * state_size;
  qsort(pointers, chain_size, sizeof(double *), SortIndexesAscending);
  for (i = 0; i < chain_size; i++)
  {
    k = pointers[i] - states;
    m = k / state_size;
    proposal = proposals + k;
    for (j = i; j < chain_size; j++)
    {
      n = (pointers[j] - states) / state_size;
      if (m < old_chain_size && n < old_chain_size)
        continue; /* Avoid re-computation of weights already computed. */
      pdf = transition_pdf(proposal, pointers[j], lower_bound, upper_bound, step, state_size);
      if (pdf <= 0.0)
        break;
      weights[m] += pdf;
    }
    for (j = 1; j <= i; j++)
    {
      n = (pointers[i-j] - states) / state_size;
      if (m < old_chain_size && n < old_chain_size)
        continue;
      pdf = transition_pdf(proposal, pointers[i-j], lower_bound, upper_bound, step, state_size);
      if (pdf <= 0.0)
        break;
      weights[m] += pdf;
    }
    if (weights[m] > 0.0)
      weights[m] = proposals_pdf[m] / weights[m];
  }
  free(pointers);
}

/** TanNormalisation.

    - Computes the PDF normalisation proposed by Tan (2006, J.Comp.Graph.Stat., 15, 735), together with an estimate of its undertainty ('error').
    - 'weights' is the array of weights (of size 'element_count') computed by 'TanWeights'. **/

double TanNormalisation(double *weights, unsigned element_count, double *error)
{
  unsigned i;
  double norm, **pointers, x;

  /* Perform weighted sum. */

  pointers = malloc(element_count * sizeof(double *));
  for (i = 0; i < element_count; i++)
    pointers[i] = weights + i;
  qsort(pointers, element_count, sizeof(double *), SortIndexesAscending);
  norm = 0.0;
  *error = 0.0;
  for (i = 0; i < element_count; i++)
  {
    x = pointers[i][0];
    norm += x;
    *error += x * x;
  }
  *error = sqrt(*error - norm * norm / element_count);
  free(pointers);
  return norm;
}

/** TanMean.

    - Computes the mean of a function as proposed by Tan (2006, J.Comp.Graph.Stat., 15, 735), together with an estimate of its undertainty ('error').
    - 'weights' is the array of weights (of size 'element_count') computed by 'TanWeights'.
    - 'norm' is the PDF normalisation computed by 'TanNormalisation'.
    - 'function' is the array of function values (of size 'element_count') at the proposal states. **/

double TanMean(double *weights, double *function, double norm, unsigned element_count, double *error)
{
  unsigned i;
  double mean, x;

  /* Perform weighted sum. */

  mean = 0.0;
  *error = 0.0;
  for (i = 0; i < element_count; i++)
    mean += weights[i] * function[i];
  mean /= norm;
  for (i = 0; i < element_count; i++)
  {
    x = (function[i] - mean) * weights[i];
    *error += x * x;
  }
  *error = sqrt(*error) / norm;
  return mean;
}

/** UniformTransition.

    - Samples an n-dimensional transition in a Markov chain, from a uniform PDF.
    - 'state' contains the present state of the Markov chain, comprising 'state_size' doubles.
    - 'proposal' will store the resulting proposed state (same size as 'state').
    - 'lower_bound' and 'upper_bound' must contain the absolute bounds for the variables defining the state.
    - 'step' must contain the parameter defining the size of the transition step. **/

void UniformTransition(double *proposal, double *state, double *lower_bound, double *upper_bound, double *step, unsigned state_size)
{
  unsigned i;
  double min, max;

  for (i = 0; i < state_size; i++)
  {
    if (step[i] <= 0.0)
      proposal[i] = state[i];
    else
    {
      min = state[i] - 0.5 * step[i];
      max = min + step[i];
      if (min < lower_bound[i])
        min = lower_bound[i];
      if (max > upper_bound[i])
        max = upper_bound[i];
      proposal[i] = min + (max - min) * gsl_rng_uniform(RandomNumberGenerator);
    }
  }
}

/** UniformTransitionPDF.

    - Computes the n-dimensional PDF value of a Markov chain transition, using a uniform PDF.
    - 'state' contains the state of the Markov chain from which the transition occurs, comprising 'state_size' doubles.
    - 'proposal' contains the proposed state of the chain after the transition (same size as 'state').
    - 'lower_bound' and 'upper_bound' must contain the absolute bounds for the variables defining the state.
    - 'step' must contain the parameter defining the size of the transition step. **/

double UniformTransitionPDF(double *proposal, double *state, double *lower_bound, double *upper_bound, double *step, unsigned state_size)
{
  unsigned i;
  double min, max, density;

  density = 1.0;
  for (i = 0; i < state_size; i++)
  {
    if (step[i] <= 0.0)
      continue;
    min = state[i] - 0.5 * step[i];
    max = min + step[i];
    if (min < lower_bound[i])
      min = lower_bound[i];
    if (max > upper_bound[i])
      max = upper_bound[i];
    if (proposal[i] < min || proposal[i] > max)
      return 0.0;
    density /= max - min;
  }
  return density;
}

/** PowerLawDeviate.

    - Generates a power-law deviate with an index 'index' in any closed interval ['min', 'max']. **/

double PowerLawDeviate(double min, double max, double index)
{
  double deviate, std_deviate;

  std_deviate = gsl_rng_uniform(RandomNumberGenerator);
  if (index == -1.0)
    deviate = min * exp(std_deviate * log(max / min));
  else if (index == 0.0)
    deviate = min + (max - min) * std_deviate;
  else
  {
    index += 1.0;
    min = pow(min, index);
    max = pow(max, index);
    deviate = pow(min + (max - min) * std_deviate, 1.0 / index);
  }
  return deviate;
}

/** IsotropicDirection.

    - Generates a random direction isotropically distributed within a cone of a given semi-aperture 'semi_aperture'.
    - The output direction has unit norm, by construction.
    - Relies on the fact that, for an isotropic distribution of directions in three dimensions, the projection onto any axis is distributed uniformly in [-1,1]. **/

void IsotropicDirection(double semi_aperture, double *direction)
{
  double cos_polar, sin_polar, azimuth, cos_polar_max;

  /* Sample polar angle cosine with uniform PDF in [cos_polar_max, 1]. */

  cos_polar_max = cos(semi_aperture);
  cos_polar = cos_polar_max + gsl_rng_uniform(RandomNumberGenerator) * (1.0 - cos_polar_max);
  sin_polar = 1.0 - cos_polar * cos_polar;
  if (sin_polar > 0.0) /* Correct eventual round-off errors. */
    sin_polar = sqrt(sin_polar);
  else
    sin_polar = 0.0;

  /* Sample azimuth with uniform PDF in [0, 2 pi). */

  azimuth = MATH_2PI * gsl_rng_uniform(RandomNumberGenerator);

  /* Compute direction. */

  direction[MATH_X] = sin_polar * cos(azimuth);
  direction[MATH_Y] = sin_polar * sin(azimuth);
  direction[MATH_Z] = cos_polar;
}

/** StandardGaussianDeviates.

    - Generates a couple of standard gaussian deviates.
    - Deviates are stored in 'deviates', which must be allocated by the calling program. **/

void StandardGaussianDeviates(double *deviates)
{
  double x, y, r;

  y = 2.0 * gsl_rng_uniform(RandomNumberGenerator) - 1.0;
  do
  {
    x = y;
    y = 2.0 * gsl_rng_uniform(RandomNumberGenerator) - 1.0;
    r = x * x + y * y;
  }
  while (r == 0.0 || r >= 1.0);
  r = sqrt(-2.0 * log(r) / r);
  deviates[0] = x * r;
  deviates[1] = y * r;
}

/** GaussianDeviates.

    - Generates a couple of gaussian deviates with mean 'mean' and standard deviation 'std_dev'.
    - Deviates are stored in 'deviates', which must be allocated by the calling program. **/

void GaussianDeviates(double *deviates, double mean, double std_dev)
{
  double x, y, r;

  y = 2.0 * gsl_rng_uniform(RandomNumberGenerator) - 1.0;
  do
  {
    x = y;
    y = 2.0 * gsl_rng_uniform(RandomNumberGenerator) - 1.0;
    r = x * x + y * y;
  }
  while (r == 0.0 || r >= 1.0);
  r = sqrt(-2.0 * log(r) / r) * std_dev;
  deviates[0] = x * r + mean;
  deviates[1] = y * r + mean;
}

/** ComplexGaussianDeviates.

    - Generates a gaussian deviate with mean 'mean' and standard deviation 'std_dev' in the complex plane. **/

complex ComplexGaussianDeviates(complex mean, complex std_dev)
{
  double x, y, r;
  complex deviate;

  y = 2.0 * gsl_rng_uniform(RandomNumberGenerator) - 1.0;
  do
  {
    x = y;
    y = 2.0 * gsl_rng_uniform(RandomNumberGenerator) - 1.0;
    r = x * x + y * y;
  }
  while (r == 0.0 || r >= 1.0);
  r = sqrt(-2.0 * log(r) / r);
  deviate = mean + r * (x * creal(std_dev) + I * y * cimag(std_dev));
  return deviate;
}

#endif

/*** End of 'mathematics.h'. ***/


text:
/***************************************************************************
 *   Copyright (C) 2006 by Reed A. Cartwright                              *
 *   reed@scit.us                                                          *
 *                                                                         *
 *   Permission is hereby granted, free of charge, to any person obtaining *
 *   a copy of this software and associated documentation files (the       *
 *   "Software"), to deal in the Software without restriction, including   *
 *   without limitation the rights to use, copy, modify, merge, publish,   *
 *   distribute, sublicense, and/or sell copies of the Software, and to    *
 *   permit persons to whom the Software is furnished to do so, subject to *
 *   the following conditions:                                             *
 *                                                                         *
 *   The above copyright notice and this permission notice shall be        *
 *   included in all copies or substantial portions of the Software.       *
 *                                                                         *
 *   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,       *
 *   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF    *
 *   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.*
 *   IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR     *
 *   OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
 *   ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR *
 *   OTHER DEALINGS IN THE SOFTWARE.                                       *
 ***************************************************************************/

#ifndef EMDEL_H
#define EMDEL_H

#include <boost/config.hpp>
#include <boost/program_options.hpp>
namespace po = boost::program_options;

#include "em.h"

class emdel_app  {
public:
	emdel_app(int argc, char *argv[]);
	virtual ~emdel_app() { }

	virtual int run();

	std::string args_comment() const;

	po::options_description desc;

	// use X-Macros to specify argument variables
	struct arg_t {
#	define XCMD(lname, sname, desc, type, def) type V(lname) ;
#	include "emdel.cmds"
#	undef XCMD
	};
protected:
	arg_t arg;

};

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

namespace std {
template<typename _Tp, typename _CharT, typename _Traits>
basic_ostream<_CharT, _Traits>&
operator<<(basic_ostream<_CharT, _Traits>& os, const std::vector<_Tp> &v)
{
	if(v.size() == 1)
	{
		os << v.front();
	}
	else if(v.size() > 1)
	{
		std::copy(v.begin(), v.end()-1, std::ostream_iterator<_Tp, _CharT, _Traits>(os, " "));
		os << v.back();
	}
	return os;
}
}

#ifdef BOOST_WINDOWS
#	include <process.h>
#	define getpid _getpid
#endif

inline unsigned int rand_seed_start()
{
	unsigned int u = static_cast<unsigned int>(time(NULL));
	unsigned int v = static_cast<unsigned int>(getpid());
	v += (v << 15) + (v >> 3); // Spread 5-decimal PID over 32-bit number
	return u ^ (v + 0x9e3779b9u + (u<<6) + (u>>2));
}

inline unsigned int rand_seed()
{
	static unsigned int u = rand_seed_start();
	return (u = u*1664525u + 1013904223u);
}

class gslrand
{
public:
	typedef unsigned int int_type;
	typedef int_type seed_type;
	typedef double real_type;
	gslrand() {
		//gsl_rng_env_setup();
		T = gsl_rng_mt19937;
		r = gsl_rng_alloc(T);
	}
	virtual ~gslrand() {
		gsl_rng_free(r);
	}

	inline void seed(seed_type s) { gsl_rng_set(r, s); }
	inline int_type get() { return gsl_rng_get(r); }
	inline int_type operator()() { return get(); }

	inline int_type uniform() { return get(); }
	inline int_type uniform(int_type max) { return gsl_rng_uniform_int(r, max); }
	inline real_type uniform01() { return gsl_rng_uniform(r); }
	inline int_type geometric(real_type p) { return gsl_ran_geometric(r,p); }

private:
	const gsl_rng_type *T;
	gsl_rng *r;
};

extern gslrand myrand;

void set_rand_seed(unsigned int u);

#endif

max_stars_repo_name: reedacartwright/emdel

max_issues_repo_name: reedacartwright/emdel

max_forks_repo_name: reedacartwright/emdel