text:
/*
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
//#include "gsl-sprng.h"
#include <mpi.h>
#include <time.h>
#include <unistd.h>

const int SAMPLENUM = 1000;
int obsnum;
gsl_rng * r;
int procnum;

typedef struct {
	int prey;
	int predator0;
	int predator1;
} lvstate;

typedef struct {
	double time;
	double rawdata;
} obsdata;

void simPrior(lvstate *simData, MPI_Datatype datatype, MPI_Comm comm);
void rowsample(int *rows, double *w);
void stepLV(lvstate *state, double *t0p, double *dtp, double *lvParam, MPI_Datatype datatype, MPI_Comm comm);
void peturb(double *lvParam, MPI_Datatype datatype, MPI_Comm comm);

double obsLik(lvstate *mystate, double obs)
{
		const double SIGMA = 20.0;
		return log(gsl_ran_gaussian_pdf((obs - mystate->prey), SIGMA));
}

void mpiStepLV(lvstate *simData, double *t0p, double *dtp, double *lvParam, MPI_Datatype datatype, MPI_Comm comm)
{
	//printf("...callStepLV...");
	int j, offset, destNode, sourceNode, rank;
	int procnum, avgsize, samplemod;
	double t0, dt;
	MPI_Status status;

	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD, &procnum);
	//const unsigned int MAX_MPI_BUFFER = 524300000;
	//printf("rank in mpiStepLV %i \n", rank);
	//for (j = 0; j < SAMPLENUM; j++) {
		//MPI_Bcast(&simData[j].predator, 1, MPI_INT, 0, comm);
		//MPI_Bcast(&simData[j].prey, 1, MPI_INT, 0, comm);
	//}
	MPI_Bcast(lvParam, 5, MPI_DOUBLE, 0, comm);
	MPI_Bcast(t0p, 1, MPI_DOUBLE, 0, comm);
	MPI_Bcast(dtp, 1, MPI_DOUBLE, 0, comm);

	avgsize = SAMPLENUM/(procnum-1);
	samplemod = SAMPLENUM%(procnum-1);

	t0 = *t0p;
	dt = *dtp;

	if (rank==0)
	{
		//send tasks to other nodes

		offset = 0;

		for (destNode=1; destNode<procnum; destNode++){
			//MPI_Send(&lvParam, 3, MPI_DOUBLE, destNode, 0, MPI_COMM_WORLD);
			MPI_Send(&offset, 1, MPI_INT, destNode, 0, MPI_COMM_WORLD);
			for (j=0;j<avgsize;j++)
			{
				MPI_Send(&simData[offset+j].prey, 1, MPI_INT, destNode, 0, MPI_COMM_WORLD);
				MPI_Send(&simData[offset+j].predator0, 1, MPI_INT, destNode, 0, MPI_COMM_WORLD);
				MPI_Send(&simData[offset+j].predator1, 1, MPI_INT, destNode, 0, MPI_COMM_WORLD);
			}
			offset = offset + avgsize;
		}

		for (j=0;j<samplemod;j++){
			stepLV(simData+offset+j, &t0, &dt, lvParam, MPI_FLOAT, MPI_COMM_WORLD);
		}

		// Wait to receive results from each node
		for (sourceNode=1; sourceNode<procnum; sourceNode++){

			MPI_Recv(&offset, 1, MPI_INT, sourceNode, 0, MPI_COMM_WORLD, &status);
			for (j=0;j<avgsize;j++)
			{
				MPI_Recv(&simData[offset+j].prey, 1, MPI_INT, sourceNode, 0, MPI_COMM_WORLD, &status);
				MPI_Recv(&simData[offset+j].predator0, 1, MPI_INT, sourceNode, 0, MPI_COMM_WORLD, &status);
				MPI_Recv(&simData[offset+j].predator1, 1, MPI_INT, sourceNode, 0, MPI_COMM_WORLD, &status);
			}
		}
	}
	else if(rank>0 && rank<procnum)
	{
		//Receive data distributed from master node
		MPI_Recv(&offset, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
		for (j=0;j<avgsize;j++)
		{
			MPI_Recv(&simData[offset+j].prey, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
			MPI_Recv(&simData[offset+j].predator0, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
			MPI_Recv(&simData[offset+j].predator1, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
		}

		//doing stepLV for every slave processor/node
		for (j=0;j<avgsize;j++)
		{
			stepLV(simData+offset+j, &t0, &dt, lvParam, MPI_FLOAT, MPI_COMM_WORLD);
		}
		//printf("---5---");
		MPI_Send(&offset, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
		for (j=0;j<avgsize;j++)
		{
			MPI_Send(&simData[offset+j].prey, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
			MPI_Send(&simData[offset+j].predator0, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
			MPI_Send(&simData[offset+j].predator1, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
		}

	}
	//MPI_Buffer_detach( &buff, &b_size );
	return;
}

//particle filter: output new estimated lvstate and its likehood
void pfPropPar(lvstate ppstate[16], obsdata *myobsdata, double *lvParam, double *ll, MPI_Datatype datatype, MPI_Comm comm)
{
	int rank, procnum;
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD, &procnum);

	typedef struct {
		lvstate *vpptuple;
	} stateVec;

	int i, j, k, l;

	double timeUnit = 2.0;
	int deltaTimeUnit = 1;
	double delTime = timeUnit * deltaTimeUnit;
	double curTime = 0.0;

	double lw[SAMPLENUM];
	double w[SAMPLENUM];
	double wSum;
	double maxLw;
	double likehood = 0.0;
	unsigned int rows[SAMPLENUM];

	lvstate *simData = malloc(SAMPLENUM * sizeof(lvstate));
	if (simData == NULL) {
		printf("Fail to initiate simData...EXIT!");
		goto freeMem;
	}
	else {
		memset(simData, 0, SAMPLENUM * sizeof(lvstate));
	}
	//initiate state: get simulated prey-predator tuple
		stateVec *myStateVec = malloc(obsnum * sizeof(stateVec));
		if (myStateVec == NULL) {
			printf("Fail to initiate myStateVec...Exit");
			goto freeMem;
		}
		memset(myStateVec, 0, obsnum * sizeof(stateVec));

		stateVec *tempStateVec = malloc(obsnum * sizeof(stateVec));
		if (tempStateVec == NULL) {
			printf("Fail to initiate tempStateVec...Exit");
			goto freeMem;
		}
		memset(tempStateVec, 0, obsnum * sizeof(stateVec));

	if (rank==0)
	{
		simPrior(simData, MPI_FLOAT, MPI_COMM_WORLD);
		*ll = 0.0;

		if ((tempStateVec + 0)->vpptuple == NULL) {
			(tempStateVec + 0)->vpptuple = (lvstate*) malloc(
					SAMPLENUM * sizeof(lvstate));
			if ((tempStateVec + 0)->vpptuple == NULL) {
				printf("Fail to initiate myStateVec->vpptuple...Exit");
				goto freeMem;
			}
		}

		memset(tempStateVec[0].vpptuple, 0, SAMPLENUM * sizeof(lvstate));
		for (j = 0; j < SAMPLENUM; j++) {
			((tempStateVec + 0)->vpptuple + j)->predator0 = (simData + j)->predator0;
			((tempStateVec + 0)->vpptuple + j)->predator1 = (simData + j)->predator1;
			((tempStateVec + 0)->vpptuple + j)->prey = (simData + j)->prey;
		}
	}

	for (i = 1; i < obsnum; i++)
	{
		curTime = myobsdata[i-1].time;
		wSum = 0.0;

		if (rank==0)
		{
			if ((tempStateVec + i)->vpptuple == NULL) {
				(tempStateVec + i)->vpptuple = (lvstate*) malloc(
						SAMPLENUM * sizeof(lvstate));
				if ((tempStateVec + i)->vpptuple == NULL) {
					printf("Fail to initiate myStateVec->vpptuple...Exit");
					goto freeMem;
				}
				memset((tempStateVec + i)->vpptuple, 0,
						SAMPLENUM * sizeof(lvstate));
			}
		}

		mpiStepLV(simData, &curTime, &delTime, lvParam, MPI_FLOAT, MPI_COMM_WORLD);

		if (rank==0)
		{
			for (j = 0; j < SAMPLENUM; j++)
			{
				((tempStateVec + i)->vpptuple + j)->predator0 = (simData + j)->predator0;
				((tempStateVec + i)->vpptuple + j)->predator1 = (simData + j)->predator1;
				((tempStateVec + i)->vpptuple + j)->prey = (simData + j)->prey;
				lw[j] = obsLik(((tempStateVec + i)->vpptuple + j), myobsdata[i].rawdata);
			}

			//get the max line weight
			for (j = 0; j < SAMPLENUM; j++) {
				if (j == 0) {
					maxLw = lw[0];
				} else {
					if (lw[j] > maxLw) {
						maxLw = lw[j];
					}
				}
			}
			//normalize the line wight and output list w[]
			for (l = 0; l < SAMPLENUM; l++) {
				w[l] = exp(lw[l] - maxLw);
				wSum = wSum + w[l];
			}
			rowsample(rows, w);

			//reorganize the vector with only keeping the rows in 'rows'
			for (k = 0; k <= i; k++)
			{
				if ((myStateVec + k)->vpptuple == NULL) {
					//printf("init mystatevec  \n");
					(myStateVec + k)->vpptuple = (lvstate*) malloc(
							SAMPLENUM * sizeof(lvstate));
					if ((myStateVec + k)->vpptuple == NULL) {
						printf("Fail to initiate myStateVec->vpptuple...Exit");
						goto freeMem;
					}
				}

				for (j = 0; j < SAMPLENUM; j++) {
					((myStateVec + k)->vpptuple + j)->predator0 =
							((tempStateVec + k)->vpptuple + rows[j])->predator0;
					((myStateVec + k)->vpptuple + j)->predator1 =
							((tempStateVec + k)->vpptuple + rows[j])->predator1;
					((myStateVec + k)->vpptuple + j)->prey =
							((tempStateVec + k)->vpptuple + rows[j])->prey;
				}
				for (j = 0; j < SAMPLENUM; j++) {
					((tempStateVec + k)->vpptuple + j)->predator0 =
							((myStateVec + k)->vpptuple + j)->predator0;
					((tempStateVec + k)->vpptuple + j)->predator1 =
							((myStateVec + k)->vpptuple + j)->predator1;
					((tempStateVec + k)->vpptuple + j)->prey =
							((myStateVec + k)->vpptuple + j)->prey;
				}
			}
			for (j = 0; j < SAMPLENUM; j++) {
				simData[j].predator0 = ((myStateVec + i)->vpptuple + j)->predator0;
				simData[j].predator1 = ((myStateVec + i)->vpptuple + j)->predator1;
				simData[j].prey = ((myStateVec + i)->vpptuple + j)->prey;
			}
			likehood = likehood + maxLw + log(wSum/SAMPLENUM);
			//printf("likehood: %f, maxLW: %f, mean: %f, wSum: %f\n", likehood, maxLw, log(wSum/SAMPLENUM), wSum);
			*ll = likehood;
			rows[SAMPLENUM] = 0;
		}

	}
	if (rank==0)
	{
		for (i = 0; i < obsnum; i++) {
			ppstate[i].prey = ((myStateVec + i)->vpptuple + 0)->prey;
			ppstate[i].predator0 = ((myStateVec + i)->vpptuple + 0)->predator0;
			ppstate[i].predator1 = ((myStateVec + i)->vpptuple + 0)->predator1;
		}
	}

freeMem:

	if (myStateVec != NULL) {
		for (i = 0; ((myStateVec + i)->vpptuple != NULL) && i < obsnum; i++) {
			free((myStateVec + i)->vpptuple);
			(myStateVec + i)->vpptuple = NULL;
		}
		free(myStateVec);
		myStateVec = NULL;
	}
	if (tempStateVec != NULL) {
		for (i = 0; NULL != (tempStateVec + i)->vpptuple && i < obsnum; i++) {
			free((tempStateVec + i)->vpptuple);
			(tempStateVec + i)->vpptuple = NULL;
		}
		free(tempStateVec);
		tempStateVec = NULL;
	}
	if (simData != NULL) {
		free(simData);
		simData = NULL;
	}
}

void stepLV(lvstate *state, double *t0p, double *dtp, double *lvParam, MPI_Datatype datatype, MPI_Comm comm) {
	int rank;
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD, &procnum);

	register double t0 = *t0p, dt = *dtp, t;
	register double h0, h1, h2, h3, h4, h5, u, l0, l1, l2, l3, l4;
	register int prey, predator0, predator1;
	prey = state->prey;
	predator0 = state->predator0;
	predator1 = state->predator1;
	l0 = lvParam[0];
	l1 = lvParam[1];
	l2 = lvParam[2];
	l3 = lvParam[3];
	l4 = lvParam[4];

	while (dt > 0) {
		h1 = l0 * prey;
		h2 = l1 * prey * predator0;
		h3 = l2 * prey * predator1;
		h4 = l3 * predator0;
		h5 = l4 * predator1;
		h0 = h1 + h2 + h3 + h4 + h5;

		if ((h0 < (1e-10)) || (prey >= 1000000))
			t = 1e99;
		else {
			t = gsl_ran_exponential(r, 1/h0);
		}
		if (t > dt) {
			state->prey = prey;
			state->predator0 = predator0;
			state->predator1 = predator1;
			return;
		} else {
			u = gsl_rng_uniform(r);
			if (u < (h1 / h0)) {
				prey = prey + 1;
			} else if (u < ((h1 + h2) / h0)) {
				prey = prey - 1;
				predator0 = predator0 + 1;
			} else if (u < ((h1 + h2 + h3) / h0)) {
				prey = prey - 1;
				predator1 = predator1 + 1;
			} else if (u < ((h1 + h2 + h3 + h4) / h0)) {
				predator0 = predator0 - 1;
			} else {
				predator1 = predator1 - 1;
			}
			dt = dt - t;

		}
	}
	return;
}

//mcmc process
void runPmmhPath(int its, double *lvParam, double *obslik, lvstate ppstate[16],
		obsdata *myobsdata, MPI_Datatype datatype, MPI_Comm comm)
{
	int rank;
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD, &procnum);

	int i, j;
	double propParam[5];
	double curParam[5];
	double ll;
	double propMll, curMll = -1e99;
	lvstate curPath[obsnum];
	lvstate propPath[obsnum];
	FILE *fp = fopen("mcmc-out.txt", "wab+");
	if (fp == NULL) {
		printf("file is null.\n");
		return;
	}

	if (rank==0)
	{

		fprintf(fp, "th1,th2,th3,th4,th5,");
		for (i = 0; i <= 30; i++) {
			if ((i % 2) == 0) {
				if (i == 30) {
					fprintf(fp, "x%i,y0_%i,y1_%i", i, i, i);
				} else {
					fprintf(fp, "x%i,y0_%i,y1_%i,", i, i, i);
				}
			}
		}
		fprintf(fp, "\n");
		memcpy(curPath, ppstate, sizeof(lvstate) * obsnum);
		memcpy(curParam, lvParam, sizeof(double) * 5);
	}

	for (i = its; i > 0; i--) {
		memcpy(propParam, curParam, sizeof(double) * 5);
		peturb(propParam, MPI_FLOAT, MPI_COMM_WORLD);

		pfPropPar(propPath, myobsdata, propParam, &ll, MPI_FLOAT, MPI_COMM_WORLD);

		if (rank==0)
		{
			//printf("likelihood %f, parameter %f, %f, %f, %f, %f:\n", ll, propParam[0], propParam[1], propParam[2], propParam[3], propParam[4]);
			propMll = ll;
			if (log(gsl_ran_flat(r, 0.0, 1.0)) < (propMll - curMll)) {
				curMll = propMll;
				memcpy(curParam, propParam, sizeof(double) * 5);
				memcpy(curPath, propPath, sizeof(lvstate) * obsnum);
			}
		}

		if (rank==0)
		{
			printf("run for the %ith iteration.\n", i);
			//write current parameters and state sequence into output
			fprintf(fp, "%f,%f,%f,%f,%f,", curParam[0], curParam[1], curParam[2], curParam[3], curParam[4]);

			for (j = 0; j < obsnum; j++) {
				if (j == (obsnum - 1)) {
					fprintf(fp, "%i,%i,%i", curPath[j].prey, curPath[j].predator0, curPath[j].predator1);
				} else {
					fprintf(fp, "%i,%i,%i,", curPath[j].prey, curPath[j].predator0, curPath[j].predator1);
				}
			}
			fprintf(fp, "\n");
			printf("end for the %ith iteration.\n", i);
		}

	}
	fclose(fp);

}

void rowsample(int *rows, double *w)
{
	gsl_ran_discrete_t * grdp;
	int row;
	int i;
	grdp = gsl_ran_discrete_preproc(SAMPLENUM, w);

	for (i = 0; i < SAMPLENUM; i++) {
		row = (int) gsl_ran_discrete(r, grdp);
		rows[i] = row;
	}

	gsl_ran_discrete_free(grdp);
	return;

}

void peturb(double *lvParam, MPI_Datatype datatype, MPI_Comm comm) {
	int rank;
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	if (rank==0)
	{
		const double SIGMA = 0.035;
		lvParam[0] = lvParam[0] * exp(gsl_ran_gaussian(r, SIGMA));
		lvParam[1] = lvParam[1] * exp(gsl_ran_gaussian(r, SIGMA));
		lvParam[2] = lvParam[2] * exp(gsl_ran_gaussian(r, SIGMA));
		lvParam[3] = lvParam[3] * exp(gsl_ran_gaussian(r, SIGMA));
		lvParam[4] = lvParam[4] * exp(gsl_ran_gaussian(r, SIGMA));
	}
	return;
}

//to get the prior simulate value of prey and predator
void simPrior(lvstate *simData, MPI_Datatype datatype, MPI_Comm comm)
{
	int rank;
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);

	int i;

	const double PREY_MEAN = 800.0;
	const double PRED0_MEAN = 500.0;
	const double PRED1_MEAN = 600.0;

	for (i = 0; i < SAMPLENUM; i++) {
		simData[i].prey = gsl_ran_poisson(r, PREY_MEAN);
		simData[i].predator0 = gsl_ran_poisson(r, PRED0_MEAN);
		simData[i].predator1 = gsl_ran_poisson(r, PRED1_MEAN);
	}
	return;
}

void runModel(int its, MPI_Datatype datatype, MPI_Comm comm) {
	int rank;
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);

	obsnum = 16;
	obsdata myobsdata[obsnum];
	lvstate mylvstate[obsnum];
	double ll = 0.0;
	double lvParam[5] = { 10.0, 0.005, 0.0025, 6.0, 3.0 };

	//produce time list(with even number in)
	if (rank==0) {
		int i;
		int j = 0;

		for (i = 0; i <= 30; i++) {
			if ((i % 2) == 0) {
				myobsdata[j].time = (double) i;
				j++;
			}
		}

		//read external txt file
		FILE *file = fopen("LV12.txt", "r");
		char line[1024];
		if (file == NULL) {
			printf("file is null.\n");
			exit(EXIT_FAILURE);
		}
		if (file != NULL) {
			i = 0;
			while (fgets(line, 1024, file)) {
				myobsdata[i].rawdata = atof(line);
				i++;
			}
			fclose(file);
		}
	}

	runPmmhPath(its, lvParam, &ll, mylvstate, myobsdata, MPI_FLOAT, MPI_COMM_WORLD);
	return;
}

int main(int argc,char *argv[])
{
	int its, i, rank;
	double x, can, a, alpha, Nsum, mytime, maxtime, mintime;
	clock_t begin, end;
	double time_spent;

	begin = clock();
	const gsl_rng_type * T;
	long seed;
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	seed = time(NULL) * getpid();    // set seed
	gsl_rng_set(r, seed);

	MPI_Init(NULL, NULL);
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD, &procnum);

	if (argc == 1) {
		its = 50;
	}
	else {
		its = atoi(argv[1]);
	}

	runModel(its, MPI_FLOAT, MPI_COMM_WORLD);

	if (rank==0) {
		printf("Starting lv-pmcmc ...\n");
		printf("Running for %i", its);
		printf(" iterations, Done.\n");
	}

	gsl_rng_free (r);

	// here, do your time-consuming job
	end = clock();
	time_spent = (double)(end - begin) / CLOCKS_PER_SEC;

	if (rank==0){
		printf("Time consuming in total is %f\n", time_spent);
	}

	MPI_Finalize();
    return(EXIT_SUCCESS);
}
*/

max_stars_repo_name: scania/scania-pmcmc

max_issues_repo_name: scania/scania-pmcmc

max_forks_repo_name: scania/scania-pmcmc


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

#pragma once

#include <gsl/gsl-lite.hpp>
#include <units/bits/external/hacks.h>
#include <utility>

inline constexpr struct validated_tag {
} validated;

template<std::movable T, std::predicate<T> Validator>
class validated_type {
  T value_;
public:
  using value_type = T;

  static constexpr bool validate(const T& value) { return Validator()(value); }

  constexpr explicit validated_type(const T& value) noexcept(std::is_nothrow_copy_constructible_v<T>)
    requires std::copyable<T>
  : value_(value)
  {
    gsl_Expects(validate(value_));
  }

  constexpr explicit validated_type(T&& value) noexcept(std::is_nothrow_move_constructible_v<T>) :
      value_(std::move(value))
  {
    gsl_Expects(validate(value_));
  }

  constexpr validated_type(const T& value, validated_tag) noexcept(std::is_nothrow_copy_constructible_v<T>)
    requires std::copyable<T>
  : value_(value)
  {
  }

  constexpr validated_type(T&& value, validated_tag) noexcept(std::is_nothrow_move_constructible_v<T>) :
      value_(std::move(value))
  {
  }

  constexpr explicit(false) operator T() const noexcept(std::is_nothrow_copy_constructible_v<T>)
    requires std::copyable<T>
  {
    return value_;
  }

  constexpr T& value() & noexcept = delete;
  constexpr const T& value() const& noexcept { return value_; }
  constexpr T&& value() && noexcept { return std::move(value_); }
  constexpr const T&& value() const&& noexcept { return std::move(value_); }

  bool operator==(const validated_type&) const
    requires std::equality_comparable<T>
  = default;
  auto operator<=>(const validated_type&) const
    requires std::three_way_comparable<T>
  = default;
};

ext: h

max_stars_repo_name: chiphogg/units

max_issues_repo_name: chiphogg/units

max_forks_repo_name: chiphogg/units


text:
#ifndef __UTIL_H
#define __UTIL_H

#include <glib.h>
#include <gsl/gsl_matrix.h>

struct _CLArgumentEntry {
  GString *block;
  char *label;
};

struct _LabeledMatrix {
  gsl_matrix *mat;
  char **labels1;
  char **labels2;
};

typedef struct _LabeledMatrix LabeledMatrix ;
typedef struct _CLArgumentEntry CLArgumentEntry ;

char *complearn_make_temp_dir_name(void);
char *complearn_make_temp_dir(void);
int make_directory(const char *name);
int complearn_make_directory_if_necessary(const char *name);
GString *complearn_read_whole_file(const char *fname);
int complearn_count_strings(const char * const * str);
int complearn_remove_directory_recursively(const char *dirname, int f_err_out);
GString *complearn_read_whole_file_ptr(FILE *fp);
GString *complearn_read_whole_file_io(GIOChannel *gio);
GSList *complearn_read_directory_of_files(const char *dirname);
GSList *complearn_read_list_of_files(const char *filename);
GSList *complearn_read_list_of_strings(const char *filename);
gsl_matrix *complearn_average_matrix(gsl_matrix *a);
gsl_matrix *complearn_svd_project(gsl_matrix *a);
CLArgumentEntry *complearn_new_arg(GString *block, const char *label);
int complearn_get_pid(void);
char *complearn_get_hostname();
char ** complearn_dupe_strings(const char * const * str);
char ** complearn_fix_labels(char **inp);
char *complearn_make_temp_file_name(void);
char *complearn_chdir(const char *newdir);
int complearn_write_file(const char *fname, const GString *f);
void complearn_mrbreaker(void);

#endif

max_stars_repo_name: rudi-cilibrasi/classic-complearn

max_issues_repo_name: rudi-cilibrasi/classic-complearn

max_forks_repo_name: rudi-cilibrasi/classic-complearn


text:
#pragma once
#include <array>
#include <gsl/span>
#include <sam.h>

namespace mcu {

#define EEPROM [[gnu::section(".eeprom")]]

class FlashEEPROM
{
public:
	FlashEEPROM() = delete;

	template <typename T>
	static void write(const T* flashAddress, const T& src) noexcept
	{
		write(flashAddress, gsl::span(&src, 1));
	}
	template <typename T>
	static T read(const T* flashAddress) noexcept
	{
		T result;
		read(flashAddress, gsl::span(&result, 1));
		return result;
	}

	template <typename T>
	static void write(const T* flashAddress, gsl::span<const T> src) noexcept
	{
		write(reinterpret_cast<const std::byte*>(flashAddress), gsl::as_bytes(src));
	}
	template <typename T>
	static void read(const T* flashAddress, gsl::span<T> dst) noexcept
	{
		read(reinterpret_cast<const std::byte*>(flashAddress), gsl::as_writable_bytes(dst));
	}

	static void write(const std::byte* flashAddress, gsl::span<const std::byte> src) noexcept;
	static void read(const std::byte* flashAddress, gsl::span<std::byte> dst) noexcept;

	static void commit() noexcept;
	static bool needsCommit() noexcept { return cacheDirty; }

private:
	using Page = std::array<unsigned, FLASH_PAGE_SIZE / sizeof(unsigned)>;
	using Row = std::array<Page, 4>;

	static constexpr size_t RowSize = sizeof(Row);
	static constexpr uintptr_t RowOffsetMask = sizeof(Row) - 1;

	union Cache {
		std::array<std::byte, sizeof(Row)> bytes;
		Row pages;
	};
	static Cache cache;
	static Row* cacheTag;
	static bool cacheDirty;

	static void cacheRow(Row* row) noexcept;
	static void copyPage(const Page& src, Page& dst) noexcept;
};

} // namespace mcu


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

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

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

struct TensorResizePara
{
    MLayer *prev;

    int height;
    int width;

    int res_valid;
};

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

    para->prev = mNetworkLayer(ini,mINIRead(ini,name,"prev"));
    mException((para->prev == NULL),EXIT,"invalid prev");
    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    para->height= DFLT; mINIRead(ini,name,"height","%d",&(para->height));
    para->width = DFLT; mINIRead(ini,name,"width" ,"%d",&(para->width ));

    return para;
}

struct HandleTensorResize
{
    int *lx;
    int *ly;
    float *wx;
    float *wy;
};
void endTensorResize(void *info)
{
    struct HandleTensorResize *handle = (struct HandleTensorResize *)info;
    if(handle->lx !=NULL) mFree(handle->lx);
    if(handle->ly !=NULL) mFree(handle->ly);
    if(handle->wx !=NULL) mFree(handle->wx);
    if(handle->wy !=NULL) mFree(handle->wy);
}
#define HASH_TensorResize 0xb6f42d3a

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

    MHandle *hdl=mHandle(out,TensorResize);
    struct HandleTensorResize *handle = (struct HandleTensorResize *)(hdl->handle);

    if(para->height<=0) para->height= in->height;
    if(para->width <=0) para->width = in->width;

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

    if(handle->lx !=NULL) {mFree(handle->lx);} handle->lx=(int   *)mMalloc(para->width *sizeof(int));
    if(handle->ly !=NULL) {mFree(handle->ly);} handle->ly=(int   *)mMalloc(para->height*sizeof(int));
    if(handle->wx !=NULL) {mFree(handle->wx);} handle->wx=(float *)mMalloc(para->width *sizeof(float));
    if(handle->wy !=NULL) {mFree(handle->wy);} handle->wy=(float *)mMalloc(para->height*sizeof(float));

    float kx = (float)(in->width)/(float)(out->width);
    for(int i=0;i<out->width;i++)
    {
        float l=kx*((float)i-0.5);
        if(l<0) l=0.0f;else if(l>=in->width-1) l=(float)(in->width)-1.00001;
        handle->lx[i] = (int)l;
        handle->wx[i] = (float)(handle->lx[i]+1)-l;
    }

    float ky = (float)(in->height)/(float)(out->height);
    for(int i=0;i<out->height;i++)
    {
        float l=ky*((float)i-0.5);
        if(l<0) l=0.0f;else if(l>=in->height-1) l=(float)(in->height)-1.00001;
        handle->ly[i] = (int)l;
        handle->wy[i] = (float)(handle->ly[i]+1)-l;
    }

    hdl->valid = 1;
}

void mTensorResizeForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Resize",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorResizePara *para = (struct TensorResizePara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *out=layer->tns;

    TensorResizeSet(layer);

    MHandle *hdl=mHandle(out,TensorResize);
    struct HandleTensorResize *handle = (struct HandleTensorResize *)(hdl->handle);

    int  in_size = in->width* in->height;
    int out_size =out->width*out->height;

    for(int h=0;h<out->height;h++)for(int w=0;w<out->width;w++)
    {
        int out_idx = h*out->width+w;
        int in_idx1 =(handle->ly[h]  )*in->width + (handle->lx[w]  );
        int in_idx2 =(handle->ly[h]  )*in->width + (handle->lx[w]+1);
        int in_idx3 =(handle->ly[h]+1)*in->width + (handle->lx[w]  );
        int in_idx4 =(handle->ly[h]+1)*in->width + (handle->lx[w]+1);
        float w1 = (    handle->wy[h])*(    handle->wx[w]);
        float w2 = (    handle->wy[h])*(1.0-handle->wx[w]);
        float w3 = (1.0-handle->wy[h])*(    handle->wx[w]);
        float w4 = (1.0-handle->wy[h])*(1.0-handle->wx[w]);

        for(int b=0;b<in->batch;b++)for(int c=0;c<in->channel;c++)
        {
            out->data[b][c*out_size+out_idx] = in->data[b][c*in_size+in_idx1]*w1
                                              +in->data[b][c*in_size+in_idx2]*w2
                                              +in->data[b][c*in_size+in_idx3]*w3
                                              +in->data[b][c*in_size+in_idx4]*w4;
        }
    }

    layer->state = MORN_FORWARD;
}

void mTensorResizeBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Resize",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorResizePara *para = (struct TensorResizePara *)(layer->para);
    if(para->res_valid==0) return;
    MTensor *res= para->prev->res;
    MTensor *out= layer->res;

    MHandle *hdl=mHandle(layer->tns,TensorResize);
    struct HandleTensorResize *handle = (struct HandleTensorResize *)(hdl->handle);
    mException((hdl->valid == 0),EXIT,"no forward operate");

    int  in_size =res->width*res->height;
    int out_size =out->width*out->height;

    if(para->prev->state == MORN_FORWARD)
    {
        for(int b=0;b<res->batch;b++) memset(res->data[b],0,in_size*res->channel*sizeof(float));
        para->prev->state = MORN_BACKWARD;
    }

    for(int h=0;h<out->height;h++)for(int w=0;w<out->width;w++)
    {
        int out_idx = h*out->width+w;
        int in_idx1 =(handle->ly[h]  )*res->width + (handle->lx[w]  );
        int in_idx2 =(handle->ly[h]  )*res->width + (handle->lx[w]+1);
        int in_idx3 =(handle->ly[h]+1)*res->width + (handle->lx[w]  );
        int in_idx4 =(handle->ly[h]+1)*res->width + (handle->lx[w]+1);
        float w1 = (    handle->wy[h])*(    handle->wx[w]);
        float w2 = (    handle->wy[h])*(1.0-handle->wx[w]);
        float w3 = (1.0-handle->wy[h])*(    handle->wx[w]);
        float w4 = (1.0-handle->wy[h])*(1.0-handle->wx[w]);

        for(int b=0;b<res->batch;b++)for(int c=0;c<res->channel;c++)
        {
            res->data[b][c*in_size+in_idx1] += out->data[b][c*out_size+out_idx]*w1;
            res->data[b][c*in_size+in_idx2] += out->data[b][c*out_size+out_idx]*w2;
            res->data[b][c*in_size+in_idx3] += out->data[b][c*out_size+out_idx]*w3;
            res->data[b][c*in_size+in_idx4] += out->data[b][c*out_size+out_idx]*w4;
        }
    }
}

ext: c

max_stars_repo_name: ishine/Morn

max_issues_repo_name: ishine/Morn

max_forks_repo_name: ishine/Morn


text:
#include "ridc.h"
#include <stdlib.h>
#include <stdio.h>
#include <cmath>
#include <gsl/gsl_linalg.h>

#ifndef _BRUSSELATOR_H_
#define _BRUSSELATOR_H_

class Brusselator_GSL : public ODE {
public:
  Brusselator_GSL(int my_neq, int my_nt, double my_ti, double my_tf, double my_dt) {
    neq = my_neq;
    nt = my_nt;
    ti = my_ti;
    tf = my_tf;
    dt = my_dt;
  }

  void rhs(double t, double *u, double *f) {
    /** user implemented rhs function, u'=rhs(t,u)
       @return (by reference) f: rhs(t,u)
       @param t current time step
       @param u solution u at time t
       @param f rhs(t,u)
    */
    double A = 1.0;
    double B = 3.0;
    double alpha = 0.02;
    double dx = 1.0/(neq+1);
    int Nx = neq/2;

    f[0] = A+u[0]*u[0]*u[Nx] -(B+1.0)*u[0] +alpha/dx/dx*(u[1]-2*u[0]+1.0);
    f[Nx-1] = A+u[Nx-1]*u[Nx-1]*u[2*Nx-1] -(B+1.0)*u[Nx-1]
      + alpha/dx/dx*(1.0-2*u[Nx-1]+u[Nx-2]);

    f[Nx] = B*u[0]- u[0]*u[0]*u[Nx] + alpha/dx/dx*(u[Nx+1]-2*u[Nx]+3.0);
    f[2*Nx-1] = B*u[Nx-1]- u[Nx-1]*u[Nx-1]*u[2*Nx-1]
      + alpha/dx/dx*(3.0-2*u[2*Nx-1]+u[2*Nx-2]);

    for (int i=1; i<Nx-1; i++) {
      f[i] = A + u[i]*u[i]*u[Nx+i] -(B+1.0)*u[i] +alpha/dx/dx*(u[i+1]-2*u[i]+u[i-1]);
      f[Nx+i] = B*u[i]- u[i]*u[i]*u[Nx+i] + alpha/dx/dx*(u[Nx+i+1]-2*u[Nx+i]+u[Nx+i-1]);
    }
  }

  void step(double t, double * u, double * unew) {
    /** user implemented step function, for advancing the solution from t to t+dt
       @return (by reference) unew: solution at time t+dt
       @param t current time step
       @param u solution u at time t
       @param unew solution at time t+dt
    */
    double tnew = t + dt;

    double NEWTON_TOL = 1.0e-14;
    int NEWTON_MAXSTEP = 1000;

    double * J;
    double * stepsize;
    double * uguess;

    stepsize = new double[neq];
    uguess = new double[neq];

    J = new double[neq*neq];

    // GSL
    gsl_matrix_view m;
    gsl_vector_view b;
    gsl_vector *x = gsl_vector_alloc(neq);
    gsl_permutation *p = gsl_permutation_alloc(neq);

    // set initial condition
    for (int i=0; i<neq; i++) {
      uguess[i] = u[i];
    }

    double maxstepsize;

    int counter=0;
    int * pivot = new int[neq];

    while (1) {

      // create rhs
      newt(tnew,u,uguess,stepsize);

      // compute numerical Jacobian
      jac(tnew,uguess,J);

      m = gsl_matrix_view_array(J,neq,neq);
      b = gsl_vector_view_array(stepsize,neq);

      int s; // presumably, this is some flag info

      gsl_linalg_LU_decomp (&m.matrix,p,&s);
      gsl_linalg_LU_solve (&m.matrix,p,&b.vector, x);

      // check for convergence
      maxstepsize = 0.0;
      for (int i = 0; i<neq; i++) {
	uguess[i] = uguess[i] - x->data[i];
	if (std::abs(stepsize[i])>maxstepsize) {
	  maxstepsize = std::abs(x->data[i]);
	}
      }

      // if update sufficiently small enough
      if ( maxstepsize < NEWTON_TOL) {
	break;
      }

      counter++;
      //error if too many steps
      if (counter > NEWTON_MAXSTEP) {
	fprintf(stderr,"max Newton iterations reached\n");
	exit(42);
      }
    } // end Newton iteration

    for (int i=0; i<neq; i++) {
      unew[i] = uguess[i];
    }

    delete [] pivot;
    delete [] uguess;
    delete [] stepsize;
    delete [] J;

    gsl_permutation_free (p);
    gsl_vector_free (x);
  }

 private:
  void newt(double t, double *uprev, double *uguess,
	    double *g){
    /** Helper function to compute the next Newton step
       for solving a system of equations

       @return (by reference) g how far from zero we are
       @param t current time step
       @param uguess current solution guess
       @param uprev solution at previous time step
       @param g how far from zero we are, returned by reference
    */
    rhs(t,uguess,g);
    for (int i =0; i<neq; i++) {
      g[i] = uguess[i]-uprev[i]-dt*g[i];
    }
  }

  void jac(double t, double *u, double *J){
    /** Helper function to the Jacobian matrix (using finite differences)
       for advancing the solution from time t(n) to t(n+1) using an
       implicit Euler step on a system of equations

       @return (by reference) J the Jacobian for the Newton step
       @param t current time step
       @param u function value at the current time step
       @param J Jacobian, returned by reference
    */
    double d = 1e-5; // numerical Jacobian approximation
    double *u1;
    double *f1;
    double *f;

    // need to allocate memory
    u1 = new double[neq];
    f1 = new double[neq];
    f = new double[neq];

    // note, instead of using newt, use rhs for efficiency
    rhs(t,u,f);

    for (int i = 0; i<neq; i++) {
      for (int j = 0; j<neq; j++) {
	u1[j] = u[j];
      }
      u1[i] = u1[i] + d;

      rhs(t,u1,f1);

      for (int j = 0; j<neq; j++) {
	J[j*neq+i] = -dt*(f1[j]-f[j])/d;
      }
      J[i*neq+i] = 1.0+J[i*neq+i];
    }

    // need to delete memory
    delete [] u1;
    delete [] f1;
    delete [] f;
  }

};

#endif // _BRUSSELATOR_H_

max_stars_repo_name: ongbw/ridc

max_issues_repo_name: ongbw/ridc

max_forks_repo_name: ongbw/ridc


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

#ifndef __GSL_SPLINALG_H__
#define __GSL_SPLINALG_H__

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_types.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

/* iteration solver type */
typedef struct
{
  const char *name;
  void * (*alloc) (const size_t n, const size_t m);
  int (*iterate) (const gsl_spmatrix *A, const gsl_vector *b,
                  const double tol, gsl_vector *x, void *);
  double (*normr)(const void *);
  void (*free) (void *);
} gsl_splinalg_itersolve_type;

typedef struct
{
  const gsl_splinalg_itersolve_type * type;
  double normr; /* current residual norm || b - A x || */
  void * state;
} gsl_splinalg_itersolve;

/* available types */
GSL_VAR const gsl_splinalg_itersolve_type * gsl_splinalg_itersolve_gmres;

/*
 * Prototypes
 */
gsl_splinalg_itersolve *
gsl_splinalg_itersolve_alloc(const gsl_splinalg_itersolve_type *T,
                             const size_t n, const size_t m);
void gsl_splinalg_itersolve_free(gsl_splinalg_itersolve *w);
const char *gsl_splinalg_itersolve_name(const gsl_splinalg_itersolve *w);
int gsl_splinalg_itersolve_iterate(const gsl_spmatrix *A,
                                   const gsl_vector *b,
                                   const double tol, gsl_vector *x,
                                   gsl_splinalg_itersolve *w);
double gsl_splinalg_itersolve_normr(const gsl_splinalg_itersolve *w);

__END_DECLS

#endif /* __GSL_SPLINALG_H__ */

max_stars_repo_name: shi-bash-cmd/qtTest

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
// Copyright [1999-2017] EMBL-European Bioinformatics Institute
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_cdf.h>

#include "multiSet.h"

typedef struct setComparisonData_st {
	Multiset * multi;
} SetComparisonData;

void SetComparisonSeek(WiggleIterator * iter, const char * chrom, int start, int finish) {
	SetComparisonData * data = (SetComparisonData* ) iter->data;
	seekMultiset(data->multi, chrom, start, finish);
	pop(iter);
}

////////////////////////////////////////////////////////
// T-test
////////////////////////////////////////////////////////

void TTestReductionPop(WiggleIterator * wi) {
	if (wi->done)
		return;

	SetComparisonData * data = (SetComparisonData *) wi->data;
	Multiset * multi = data->multi;

	if (multi->done) {
		wi->done = true;
		return;
	}

	// Go to first position where both of the sets have at least one value
	while (!multi->inplay[0] || !multi->inplay[1]) {
		popMultiset(multi);
		if (multi->done) {
			wi->done = true;
			return;
		}
	}
	wi->chrom = multi->chrom;
	wi->start = multi->start;
	wi->finish = multi->finish;

	// Compute measurements
	double sum1, sum2, sumSq1, sumSq2;
	int count1, count2;
	int index;

	sum1 = sum2 = 0;
	sumSq1 = sumSq2 = 0;
	count1 = multi->multis[0]->count;
	count2 = multi->multis[1]->count;

	for (index = 0; index < multi->multis[0]->count; index++) {
		if (multi->multis[0]->inplay[index]) {
			sum1 += multi->values[0][index];
			sumSq1 += multi->values[0][index] * multi->values[0][index];
		}
	}

	for (index = 0; index < multi->multis[1]->count; index++) {
		if (multi->multis[1]->inplay[index]) {
			sum2 += multi->values[1][index];
			sumSq2 += multi->values[1][index] * multi->values[1][index];
		}
	}

	// To avoid divisions by 0:
	if (count1 == 0 || count2 == 0) {
		wi->value = NAN;
		popMultiset(multi);
		return;
	}

	double mean1 = sum1 / count1;
	double mean2 = sum2 / count2;
	double meanSq1 = sumSq1 / count1;
	double meanSq2 = sumSq2 / count2;
	double var1 = meanSq1 - mean1 * mean1;
	double var2 = meanSq2 - mean2 * mean2;

	// To avoid divisions by 0:
	if (var1 + var2 == 0) {
		wi->value = NAN;
		popMultiset(multi);
		return;
	}

	// T-statistic

	double t = (mean1 - mean2) / sqrt(var1 / count1 + var2 / count2);

	if (t < 0)
		t = -t;

	// Degrees of freedom

	double nu = (var1 / count1 + var2 / count2) * (var1 / count1 + var2 / count2) / ((var1 * var1) / (count1 * count1 * (count1 - 1)) + (var2 * var2) / (count2 * count2 * (count2 - 1)));

	// P-value

	wi->value = 2 * gsl_cdf_tdist_Q(t, nu);

	// Update inputs
	popMultiset(multi);
}

WiggleIterator * TTestReduction(Multiset * multi) {
	SetComparisonData * data = (SetComparisonData *) calloc(1, sizeof(SetComparisonData));
	if (multi->count != 2 || multi->multis[0]->count < 3 || multi->multis[1]->count < 3) {
		puts("The t-test function only works for two sets with enough elements to compute variance");
		exit(1);
	}
	data->multi = multi;
	return newWiggleIterator(data, &TTestReductionPop, &SetComparisonSeek, NAN);
}

////////////////////////////////////////////////////////
// F-test
////////////////////////////////////////////////////////

typedef struct ftestData_st {
	Multiset * multi;
	int * counts;
	double * means;
	int total_count;
} FTestData;

void FTestSeek(WiggleIterator * iter, const char * chrom, int start, int finish) {
	FTestData * data = (FTestData* ) iter->data;
	seekMultiset(data->multi, chrom, start, finish);
	pop(iter);
}

void FTestReductionPop(WiggleIterator * wi) {
	if (wi->done)
		return;

	FTestData * data = (FTestData *) wi->data;
	Multiset * multi = data->multi;

	if (multi->done) {
		wi->done = true;
		return;
	}

	// Go to first position where both of the sets have at least one value
	while (!multi->inplay[0] || !multi->inplay[1]) {
		popMultiset(multi);
		if (multi->done) {
			wi->done = true;
			return;
		}
	}
	wi->chrom = multi->chrom;
	wi->start = multi->start;
	wi->finish = multi->finish;

	// Compute means
	double mean = 0;
	int groups = multi->count;
	int index, index2;
	for (index = 0; index < groups; index++) {
		Multiplexer * mplx = multi->multis[index];
		data->means[index] = 0;
		for (index2 = 0; index < mplx->count; index++) {
			if (mplx->inplay[index2])
				data->means[index] += mplx->values[index2];
			else
				data->means[index] += mplx->iters[index2]->default_value;
		}
		mean += data->means[index];
		data->means[index] /= data->counts[index];
	}
	mean /= data->total_count;

	double inter = 0;
	double intra = 0;
	for (index = 0; index < groups; index++) {
		Multiplexer * mplx = multi->multis[index];
		inter += mplx->count * (data->means[index] - mean) * (data->means[index] - mean);
		for (index2 = 0; index < mplx->count; index++) {
			if (mplx->inplay[index2])
				intra += (mplx->values[index2] - data->means[index]) * (mplx->values[index2] - data->means[index]);
			else
				intra += (mplx->iters[index2]->default_value - data->means[index]) * (mplx->iters[index2]->default_value - data->means[index]);
		}
	}

	// F-statistic
	inter /= groups - 1;
	intra /= data->total_count - groups;
	double f = inter / intra;

	// P-value
	wi->value = 2 * gsl_cdf_fdist_Q(f, multi->count - 1, data->total_count - multi->count);

	// Update inputs
	popMultiset(multi);
}

WiggleIterator * FTestReduction(Multiset * multi) {
	FTestData * data = (FTestData *) calloc(1, sizeof(FTestData));
	data->multi = multi;
	data->means = calloc(multi->count, sizeof(double));
	data->counts = calloc(multi->count, sizeof(int));
	int index;
	for (index = 0; index < multi->count; index++) {
		data->counts[index] = multi->multis[index]->count;
		data->total_count += data->counts[index];
	}
	return newWiggleIterator(data, &FTestReductionPop, &FTestSeek, NAN);
}

////////////////////////////////////////////////////////
// Mann-Whitney U (Wilcoxon rank-sum test)
////////////////////////////////////////////////////////

typedef struct valueSetPair_st {
	double value;
	bool set;
} ValueSetPair;

typedef struct mwuData_st {
	Multiset * multi;
	int n1;
	int n2;
	int N;
	// Pre-allocated table for sorting
	ValueSetPair * rankingTable;
	// For normal approximation
	bool normalApproximation;
	double mu_U, sigma_U;
} MWUData;

void MWUSeek(WiggleIterator * iter, const char * chrom, int start, int finish) {
	MWUData * data = (MWUData* ) iter->data;
	seekMultiset(data->multi, chrom, start, finish);
	pop(iter);
}

static int compareValueSetPairs(const void * A, const void * B) {
	ValueSetPair * vspA = (ValueSetPair *) A;
	ValueSetPair * vspB = (ValueSetPair *) B;
	if (vspA->value < vspB->value)
		return -1;
	if (vspA->value > vspB->value)
		return 1;
	return 0;
}

void MWUReductionPop(WiggleIterator * wi) {
	if (wi->done)
		return;

	MWUData * data = (MWUData *) wi->data;
	Multiset * multi = data->multi;

	if (multi->done) {
		wi->done = true;
		return;
	}

	// Go to first position where both of the sets have at least one value
	while (!multi->inplay[0] || !multi->inplay[1]) {
		popMultiset(multi);
		if (multi->done) {
			wi->done = true;
			return;
		}
	}
	wi->chrom = multi->chrom;
	wi->start = multi->start;
	wi->finish = multi->finish;

	// Compute measurements
	int index;
	ValueSetPair * vspPtr = data->rankingTable;

	for (index = 0; index < data->n1; index++) {
		if (multi->multis[0]->inplay[index])
			vspPtr->value = multi->values[0][index];
		else
			vspPtr->value = multi->multis[0]->iters[index]->default_value;
		if (isnan(vspPtr->value)) {
			wi->value = NAN;
			popMultiset(multi);
			return;
		}
		vspPtr->set = false;
		vspPtr++;
	}

	for (index = 0; index < data->n2; index++) {
		if (multi->multis[1]->inplay[index])
			vspPtr->value = multi->values[1][index];
		else
			vspPtr->value = multi->multis[1]->iters[index]->default_value;
		if (isnan(vspPtr->value)) {
			wi->value = NAN;
			popMultiset(multi);
			return;
		}
		vspPtr->set = true;
		vspPtr++;
	}

	qsort(data->rankingTable, data->N, sizeof(ValueSetPair), compareValueSetPairs);

	// Sum of ranks of elements of set 1
	double U1 = 0;
	// Rolling count of elements of set 1 seen prior on the list
	int prev = 0;
	// Warns you when you have a tie with the previously visited Value-Set pairs
	int ties = 0;
	int previousTies = 0;
	vspPtr = data->rankingTable;
	for (index = 0; index < data->N && prev < data->n1; index++) {
		if (!vspPtr->set) {
			U1 += index - prev;
			if (ties) {
				// Look for ties on the table prior to the current position and after the last occurence of an element of set 1.
				int index2;
				for (index2 = index + 1; index2 < data->N && data->rankingTable[index2].value == vspPtr->value && data->rankingTable[index2].set; index2++)
					previousTies++;
				U1 -= previousTies / 2.0;
				U1 += (ties - previousTies) / 2.0;
				if (previousTies == ties)
					previousTies = ties = 0;
			} else {
				int index2;
				// Look for ties with next values
				for (index2 = index + 1; index2 < data->N && data->rankingTable[index2].value == vspPtr->value; index2++)
					if (data->rankingTable[index2].set)
						ties++;
				if (ties)
					U1 += ties / 2.0;
			}
			prev++;
		}
		vspPtr++;
	}

	if (data->normalApproximation) {
		if (U1 > data->mu_U)
			wi->value = 2 * erf((data->mu_U - U1) / data->sigma_U);
		else
			wi->value = 2 * erf((U1 - data->mu_U) / data->sigma_U);
	}

	// Update inputs
	popMultiset(multi);
}

WiggleIterator * MWUReduction(Multiset * multi) {
	MWUData * data = (MWUData *) calloc(1, sizeof(MWUData));
	if (multi->count != 2 || multi->multis[0]->count == 0 || multi->multis[1]->count == 0) {
		puts("The Mann-Whitney U function only works for two non-empty sets");
		exit(1);
	}
	data->multi = multi;
	data->n1 = multi->multis[0]->count;
	data->n2 = multi->multis[1]->count;
	data->N = data->n1 + data->n2;
	data->rankingTable = calloc(data->N, sizeof(ValueSetPair));
	if (true) {
		// Ideally, tables could be used for small values of n1 and n2
		data->normalApproximation = true;
		data->mu_U = data->n1 * data->n2 / 2;
		data->sigma_U = sqrt(data->n1 * data->n2 * (data->n1 + data->n2 + 1) / 12);
	}
	return newWiggleIterator(data, &MWUReductionPop, &MWUSeek, NAN);
}

max_stars_repo_name: juettemann/WiggleTools

max_issues_repo_name: juettemann/WiggleTools

max_forks_repo_name: juettemann/WiggleTools


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

// params is not used at the moment but we must allow for
//   it to meet the required ADaPS function definition
void free_interpolate(void **interp, void *params) {
    if((*interp) != NULL) {
        SID_free(SID_FARG((interp_info *)(*interp))->x);
        SID_free(SID_FARG((interp_info *)(*interp))->y);
        gsl_interp_free(((interp_info *)(*interp))->interp);
        gsl_interp_accel_free(((interp_info *)(*interp))->accel);
        SID_free(SID_FARG * interp);
    }
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#include <gsl/gsl_spline.h>
#include "param.h"
#include "root_func.h"
#include "eos_pres.h"
#include "global.h"

double
root_func (double rho, void *params)
{
  struct param *root_func_params = (struct param *) params;
  double (*pres) (double, void *);
  pres = &eos_pres;
  return (tmp_pres - pres (rho, root_func_params));
}


text:
#pragma once
#include "iparam.h"
#include "iparamlist.h"
#include "iflag.h"
#include "iarg.h"
#include "iarglist.h"
#include "icommand.h"
#include "optioninfo.h"
#include "format.h"
#include <sfun/string_utils.h>
#include <cmdlime/usageinfoformat.h>
#include <gsl/gsl>
#include <utility>
#include <vector>
#include <memory>
#include <algorithm>
#include <iomanip>

namespace cmdlime::detail{
using namespace gsl;
namespace str = sfun::string_utils;

inline std::string adjustedToLineBreak(std::string line, std::string& text)
{
    if (!text.empty() && !isspace(text.front())){
        auto trimmedLine = str::trimFront(line);
        if (std::find_if(trimmedLine.begin(), trimmedLine.end(), [](auto ch){return std::isspace(ch);}) == trimmedLine.end())
            return line;
        while(!isspace(line.back())){
            text.insert(text.begin(), 1, line.back());
            line.pop_back();
        }
    }
    return line;
}

inline std::string popLine(std::string& text, std::size_t width, bool firstLine = false)
{
    auto newLinePos = text.find('\n');
    if (newLinePos != std::string::npos && newLinePos <= width){
        auto line = text.substr(0, newLinePos);
        text.erase(text.begin(), text.begin() + static_cast<int>(newLinePos + 1));
        if (!firstLine)
            line = str::trimFront(line);
        return line;
    }

    auto line = text.substr(0, width);
    if (text.size() < width)
        text.clear();
    else
        text.erase(text.begin(), text.begin() + static_cast<int>(width));
    if (!firstLine)
        line = str::trimFront(line);
    return adjustedToLineBreak(line, text);
}

template <typename T>
inline std::vector<not_null<T*>> getParamsByOptionality(const std::vector<not_null<T*>>& params, bool isOptional)
{
    auto result = std::vector<not_null<T*>>{};
    std::copy_if(params.begin(), params.end(), std::back_inserter(result),
                 [isOptional](auto param){return param->isOptional() == isOptional;});
    return result;
}

template <typename T>
const std::string& getName(T& option)
{
    return option.info().name();
}

template <typename T>
const std::string& getType(T& option)
{
    return option.info().type();
}

template <typename T>
const std::string& getDescription(T& option)
{
    return option.info().description();
}

template <FormatType formatType>
class UsageInfoCreator{
public:
    UsageInfoCreator(std::string programName,
                     UsageInfoFormat outputSettings,
                     const std::vector<not_null<IParam*>>& params,
                     const std::vector<not_null<IParamList*>>& paramLists,
                     std::vector<not_null<IFlag*>> flags,
                     std::vector<not_null<IArg*>> args,
                     IArgList* argList,
                     std::vector<not_null<ICommand*>> commands)
    : programName_(std::move(programName))
    , params_(getParamsByOptionality(params, false))
    , optionalParams_(getParamsByOptionality(params, true))
    , paramLists_(getParamsByOptionality(paramLists, false))
    , optionalParamLists_(getParamsByOptionality(paramLists, true))
    , flags_(std::move(flags))
    , args_(std::move(args))
    , argList_(argList)
    , commands_(std::move(commands))
    , outputSettings_(outputSettings)
    , maxOptionNameSize_(maxOptionNameLength() + outputSettings.columnsSpacing)
    {
    }

    std::string createDetailed()
    {
        return minimizedUsageInfo() +
               argsInfo() +
               paramsInfo() +
               paramListsInfo() +
               optionsInfo() +
               optionalParamListsInfo() +
               flagsInfo() +
               commandsInfo();
    }

    std::string create()
    {
        return usageInfo();
    }

private:
    using OutputFormatter = typename Format<formatType>::outputFormatter;

    std::string usageInfo()
    {
        auto result = "Usage: " + programName_ + " ";
        if (!commands_.empty())
            result += "[commands] ";

        for (auto arg : args_)
            result += OutputFormatter::argUsageName(*arg) + " ";

        for (auto param : params_)
            result += OutputFormatter::paramUsageName(*param) + " ";
        for (auto paramList : paramLists_)
            result += OutputFormatter::paramListUsageName(*paramList) + " ";
        for (auto param : optionalParams_)
           result += OutputFormatter::paramUsageName(*param) + " ";
        for (auto paramList : optionalParamLists_)
           result += OutputFormatter::paramListUsageName(*paramList) + " ";

        for (auto flag : flags_)
            result += OutputFormatter::flagUsageName(*flag) + " ";

        if (argList_)
            result += OutputFormatter::argListUsageName(*argList_);

        result += "\n";
        return result;
    }

    std::string minimizedUsageInfo()
    {
        auto result = "Usage: " + programName_ + " ";

        if (!commands_.empty())
            result += "[commands] ";

        for (auto arg : args_)
            result += OutputFormatter::argUsageName(*arg) + " ";

        for (auto param : params_)
            result += OutputFormatter::paramUsageName(*param) + " ";

        for (auto paramList : paramLists_)
            result += OutputFormatter::paramListUsageName(*paramList) + " ";

        if (!optionalParams_.empty() || !optionalParamLists_.empty())
            result += "[params] ";

        if (!flags_.empty())
            result += "[flags] ";

        if (argList_)
            result += OutputFormatter::argListUsageName(*argList_);

        result += "\n";
        return result;
    }

    std::string paramsInfo()
    {
        auto result = std::string{"Parameters:\n"};
        if (params_.empty())
            return result;
        for (const auto& param : params_){
            const auto name = OutputFormatter::paramDescriptionName(*param, outputSettings_.nameIndentation) + "\n";
            result += makeConfigFieldInfo(name, getDescription(*param));
        }
        return result;
    }

    std::string paramListsInfo()
    {
        auto result = std::string{};
        if (paramLists_.empty())
            return result;
        for (const auto paramList : paramLists_){
            const auto name = OutputFormatter::paramListDescriptionName(*paramList, outputSettings_.nameIndentation) + "\n";
            auto description = getDescription(*paramList);
            if (!description.empty())
                description += "\n(multi-value)";
            else
                description += "multi-value";
            result += makeConfigFieldInfo(name, description);
        }
        return result;
    }

    std::string optionsInfo()
    {
        if (optionalParams_.empty())
            return {};
        auto result = std::string{};
        for (const auto option : optionalParams_){
            auto description = getDescription(*option);
            if (!description.empty()){
                if (!option->defaultValue().empty())
                    description += "\n(optional, default: " + option->defaultValue() + ")";
                else
                    description += "\n(optional)";
            }
            else{
                if (!option->defaultValue().empty())
                    description += "optional, default: " + option->defaultValue();
                else
                    description += "optional";
            }
            result += makeConfigFieldInfo(OutputFormatter::paramDescriptionName(*option, outputSettings_.nameIndentation) + "\n", description);
        }
        return result;
    }

    std::string optionalParamListsInfo()
    {
        if (optionalParamLists_.empty())
            return {};
        auto result = std::string{};
        for (const auto option : optionalParamLists_){
            auto description = getDescription(*option);
            if (!description.empty()){
                description += "\n(multi-value, ";
                if (!option->defaultValue().empty())
                    description += "optional, default: " + option->defaultValue() + ")";
                else
                    description += "optional)";
            }
            else{
                description += "multi-value, ";
                if (!option->defaultValue().empty())
                    description += "optional, default: " + option->defaultValue();
                else
                    description += "optional";
            }
            result += makeConfigFieldInfo(OutputFormatter::paramListDescriptionName(*option, outputSettings_.nameIndentation) + "\n", description);
        }
        return result;
    }

    std::string argsInfo()
    {
        if (args_.empty() && !argList_)
            return {};
        auto result = std::string{"Arguments:\n"};
        for (const auto arg : args_)
            result += makeConfigFieldInfo(OutputFormatter::argDescriptionName(*arg, outputSettings_.nameIndentation) + "\n", getDescription(*arg));

        if (argList_){
            auto description = getDescription(*argList_);
            if (!description.empty()){
                description += "\n(multi-value";
                if (argList_->isOptional()){
                    if (!argList_->defaultValue().empty())
                        description += ", optional, default: " + argList_->defaultValue() + ")";
                    else
                        description +=  ", optional)";
                }
                else
                    description += ")";
            }
            else{
                description += "multi-value";
                if (argList_->isOptional()){
                    if (!argList_->defaultValue().empty())
                        description += ", optional, default: " + argList_->defaultValue();
                    else
                        description +=  ", optional";
                }
            }
            result += makeConfigFieldInfo(OutputFormatter::argListDescriptionName(*argList_, outputSettings_.nameIndentation) + "\n", description);
        }
        return result;
    }

    std::string flagsInfo()
    {
        if (flags_.empty())
            return {};
        auto result = std::string{"Flags:\n"};
        for (const auto flag : flags_)
            result += makeConfigFieldInfo(OutputFormatter::flagDescriptionName(*flag, outputSettings_.nameIndentation) + "\n", getDescription(*flag));
        return result;
    }

    std::string commandsInfo()
    {
        if (commands_.empty())
            return {};
        auto result = std::string{"Commands:\n"};
        for (const auto command : commands_){
            auto nameStream = std::stringstream{};
            if (outputSettings_.nameIndentation)
                nameStream << std::setw(outputSettings_.nameIndentation) << " ";
            nameStream << getName(*command)
                       << " [options]";
            result += makeConfigFieldInfo(nameStream.str(), getDescription(*command));
        }
        return result;
    }

    int maxOptionNameLength()
    {
        auto length = 0;
        auto updateLength = [&length](std::string name)
        {
            auto firstLine = true;
            do {
               auto nameLine = popLine(name, 100, firstLine);
               length = std::max(length, static_cast<int>(nameLine.size()));
               firstLine = false;
            } while(!name.empty());
        };

        for (auto param : params_)
            updateLength(OutputFormatter::paramDescriptionName(*param, outputSettings_.nameIndentation));
        for (auto option : optionalParams_)
            updateLength(OutputFormatter::paramDescriptionName(*option, outputSettings_.nameIndentation));
        for (auto flag : flags_)
            updateLength(OutputFormatter::flagDescriptionName(*flag, outputSettings_.nameIndentation));
        for (auto arg : args_)
            updateLength(OutputFormatter::argDescriptionName(*arg, outputSettings_.nameIndentation));
        if (argList_)
            updateLength(OutputFormatter::argListDescriptionName(*argList_, outputSettings_.nameIndentation));
        return length;
    }

    std::string makeConfigFieldInfo(std::string name, std::string description)
    {
        auto maxNameWidth = std::min(outputSettings_.maxNameColumnWidth, maxOptionNameSize_);
        const auto columnSeparatorWidth = 1;
        const auto leftColumnWidth = columnSeparatorWidth + maxNameWidth;
        const auto rightColumnWidth = static_cast<std::size_t>(outputSettings_.terminalWidth - leftColumnWidth);
        const auto descriptionWidth = rightColumnWidth - 2;
        auto stream = std::stringstream{};
        auto firstLine = true;
        while (!name.empty()){
            auto nameLine = popLine(name, static_cast<std::size_t>(maxNameWidth), firstLine);
            const auto descriptionLine = popLine(description, descriptionWidth, firstLine);
            if (firstLine)
                stream << std::setw(maxNameWidth) << std::left << nameLine << " " << descriptionLine << std::endl;
            else
                stream << std::setw(maxNameWidth) << std::left << nameLine << "   " << descriptionLine << std::endl;
            firstLine = false;
        }
        while(!description.empty()){
            auto descriptionLine = popLine(description, descriptionWidth);
            stream << std::setw(leftColumnWidth) << " " << "  " << descriptionLine << std::endl;
        }
        return stream.str();
    }

private:
    std::string programName_;
    std::vector<not_null<IParam*>> params_;
    std::vector<not_null<IParam*>> optionalParams_;
    std::vector<not_null<IParamList*>> paramLists_;
    std::vector<not_null<IParamList*>> optionalParamLists_;
    std::vector<not_null<IFlag*>> flags_;
    std::vector<not_null<IArg*>> args_;
    IArgList* argList_;
    std::vector<not_null<ICommand*>> commands_;
    UsageInfoFormat outputSettings_;
    int maxOptionNameSize_;
};

}

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_multiset.h>

BC_INLINE2(gsl_multiset_get,gsl_multiset*,size_t,size_t)


text:
#pragma once

#include <d3d11.h>
#include <DirectXMath.h>
#include <map>
#include <gsl\gsl>
#include "DrawableGameComponent.h"
#include "GaussianBlur.h"
#include "FullScreenRenderTarget.h"
#include "FullScreenQuad.h"
#include <winrt\Windows.Foundation.h>

namespace Library
{
	class FullScreenQuadMaterial;

	struct BloomSettings
	{
		float BloomThreshold{ 0.45f };
		float BlurAmount{ 2.0f };
		float BloomIntensity{ 1.25f };
		float BloomSaturation{ 1.0f };
		float SceneIntensity{ 1.0f };
		float SceneSaturation{ 1.0f };
	};

	enum class BloomDrawModes
	{
		Normal = 0,
		GlowMap,
		BlurredGlowMap,
		End
	};

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

	public:
		Bloom(Game& game, const BloomSettings& bloomSettings = DefaultBloomSettings);
		Bloom(const Bloom&) = delete;
		Bloom(Bloom&&) = default;
		Bloom& operator=(const Bloom&) = delete;
		Bloom& operator=(Bloom&&) = default;
		~Bloom() = default;

		ID3D11ShaderResourceView* SceneTexture() const;
		void SetSceneTexture(winrt::com_ptr<ID3D11ShaderResourceView> sceneTexture);

		const BloomSettings& GetBloomSettings() const;
		void SetBloomSettings(const BloomSettings& bloomSettings);

		BloomDrawModes DrawMode() const;
		const std::string& DrawModeString() const;
		void SetDrawMode(BloomDrawModes drawMode);

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

		inline static const BloomSettings DefaultBloomSettings{ 0.45f, 2.0f, 1.25f, 1.0f, 1.0f, 1.0f };

	private:
		static const std::map<BloomDrawModes, std::string> DrawModeDisplayNames;

		struct PixelCBufferPerObject
		{
			float BloomThreshold;
			float BloomIntensity;
			float BloomSaturation;
			float SceneIntensity;
			float SceneSaturation;
			DirectX::XMFLOAT3 Padding;

			PixelCBufferPerObject(const BloomSettings& bloomSettings) :
				BloomThreshold(bloomSettings.BloomThreshold), BloomIntensity(bloomSettings.BloomIntensity),
				BloomSaturation(bloomSettings.BloomSaturation), SceneIntensity(bloomSettings.SceneIntensity),
				SceneSaturation(bloomSettings.SceneSaturation), Padding()
			{
			}

			PixelCBufferPerObject& operator=(const BloomSettings& bloomSettings)
			{
				BloomThreshold = bloomSettings.BloomThreshold;
				BloomIntensity = bloomSettings.BloomIntensity;
				BloomSaturation = bloomSettings.BloomSaturation;
				SceneIntensity = bloomSettings.SceneIntensity;
				SceneSaturation = bloomSettings.SceneSaturation;

				return *this;
			}
		};

		enum class BloomShaderClass
		{
			Extract = 0,
			Composite,
			NoBloom
		};

		void DrawNormal(const GameTime& gameTime);
		void DrawGlowMap(const GameTime& gameTime);
		void DrawBlurredGlowMap(const GameTime& gameTime);

		FullScreenQuad mFullScreenQuad;
		FullScreenRenderTarget mRenderTarget;
		GaussianBlur mGaussianBlur;
		winrt::com_ptr<ID3D11ShaderResourceView> mSceneTexture;
		winrt::com_ptr<ID3D11Buffer> mPixelCBufferPerObject;
		PixelCBufferPerObject mPixelCBufferPerObjectData;
		BloomDrawModes mDrawMode{ BloomDrawModes::Normal };
		std::map<BloomDrawModes, std::function<void(const GameTime& gameTime)>> mDrawFunctions;
		std::map<BloomShaderClass, winrt::com_ptr<ID3D11ClassInstance>> mShaderClassInstances;
		BloomSettings mBloomSettings;
	};
}


text:
#ifndef SIMULATION_STEPS

#define SIMULATION_STEPS

#include <gsl/gsl_matrix.h>

extern void calculateLIFTVolMat(int N, int Np, gsl_matrix **LIFT, gsl_matrix **VolMat, gsl_matrix** MassMatrix);

extern void store_mesh(char *Mesh);
extern void initialize();
extern void time_evolution(double FinalTime);
extern void boundary_conditions();
#endif

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


text:
#pragma once

#include <gsl/gsl_assert>
#include <src/util/Logging.h>

#ifndef NDEBUG
    //#define MI_EXPECTS(x) Expects(x)
    //#define MI_ENSURES(x) Ensures(x)

#define MI_EXPECTS(x) \
  do { \
    if (!(x)) { \
        LOG(ERROR) << "assertion failed: " << #x; \
        std::terminate(); \
    } \
  } while (0)

#define MI_ENSURES(x) \
  do { \
    if (!(x)) { \
        LOG(ERROR) << "assertion failed: " << #x; \
        std::terminate(); \
    } \
  } while (0)
#else
    #define MI_EXPECTS(x)
    #define MI_ENSURES(x)
#endif


text:
#include <math.h>
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_linalg.h>
#include <getopt.h>
#include <string.h>

#include "libhdr"
#include "Library_DFE_v1.8.h"

/*
  version 1.7.3:
  -added: constant population support
  -added: generate gamma distribution deviates
  -added: generate beta distribution deviates
  -added: generate gamma+beta distribution deviates

  to compile:

  gcc -O3 -o bin/Multi_DFE_gen_v1.8 Multi_DFE_gen_v1.8.c Library_DFE_v1.8.c tmatrix_routines.c genlib.c nrlib.c nrutil.c  -lm -lgsl -lgslcblas -w

  example run:

  GSL_RNG_SEED=1 ~/Multi_DFE_est/source/bin/Multi_DFE_gen_v1.8 -N1 100 -N2 100 -nalleles 20 -t 100 -f0 0.9 -neutral 1000000 -selected 1000000 -mode 9 -conpop 1 -exp_mean 0.05 -file 1.out

  #-gamma_alpha 10 -gamma_beta 0.5
  #-exp_mean 0.02 -beta_alpha 10 -beta_beta 1
  #-exp_mean 0.02 -gamma_alpha 1000 -gamma_beta 50
  #-gamma_alpha 70 -gamma_beta 0.07 -gamma2_alpha 1000 -gamma2_beta 50
  #-exp_mean 0.05

*/

/*Global variable declaration*/
int selmode=0;
int  n_sfs=0;
int nspikes=1;//default a single spike

int N1,N2=100;
int max_n2d=2000;
double saveresults[100];
int n2_step=0, conpop=0, output_egf_mode=0;

/*Function list*/
double load_FV(int n1,int n2,int t2,double s,double f0,double *mean_FV);

/******************************************************************************/
double load_FV(int n1,int n2,int t2,double s,double f0,double *mean_FV)
{

  static double egf_vec1_lower[maxnd+1], egf_vec2_lower[maxnd+1],
    egf_vec1_upper[maxnd+1], egf_vec2_upper[maxnd+1],
    egf_vec1[maxnd+1], egf_vec2[maxnd+1],egf_vec[maxnd+1];

  double *gamma_density_vec;

  int i=0,j=0,file_size_bytes=0;
  int t2_upper=0;
  int t2_real=t2;
  int n1d=2*n1;
  int n2d=2*n2;

  char *buffer_p1_t2_lower, *buffer_p2_t2_lower,
    *buffer_p1_t2_upper, *buffer_p2_t2_upper, *buffer_const_pop;
  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  if (!conpop)
    {
      get_upper_lower_int(t2_real, &t2_lower, &t2_upper, n_t2_evaluated, t2_evaluated_vec);
      if ((t2_lower==undefined_int)||(t2_upper==undefined_int))
	{
	  return undefined;
	}
      //printf("\n%d",t2_lower);
      //printf("\n%f",s);
      file_size_bytes=compute_file_size_bytes(n2);
      //Read buffers
      ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
      buffer_p1_t2_lower = (char*) malloc (file_size_bytes);
      buffer_p1_t2_upper = (char*) malloc (file_size_bytes);

      buffer_p2_t2_lower = (char*) malloc (file_size_bytes);
      buffer_p2_t2_upper = (char*) malloc (file_size_bytes);

      read_phase1_phase2_file_into_buffer(n1,1, n2, t2_lower,
					  buffer_p1_t2_lower, file_size_bytes);
      read_phase1_phase2_file_into_buffer(n1,1, n2, t2_upper,
					  buffer_p1_t2_upper, file_size_bytes);

      read_phase1_phase2_file_into_buffer(n1,2, n2, t2_lower,
					  buffer_p2_t2_lower, file_size_bytes);
      read_phase1_phase2_file_into_buffer(n1,2, n2, t2_upper,
					  buffer_p2_t2_upper, file_size_bytes);

      get_binary_egf_vec(buffer_p1_t2_lower, n2, t2_lower, s, egf_vec1_lower);
      get_binary_egf_vec(buffer_p1_t2_upper, n2, t2_upper, s, egf_vec1_upper);

      get_binary_egf_vec(buffer_p2_t2_lower, n2, t2_lower, s, egf_vec2_lower);
      get_binary_egf_vec(buffer_p2_t2_upper, n2, t2_upper, s, egf_vec2_upper);
      //End reading buffers
      ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
      //Averaging vectors w(s) and x(s)

      compute_weighted_average_egf_vec(t2_real, t2_lower, t2_upper,
				       egf_vec1_lower, egf_vec1_upper, egf_vec1, n2d);
      compute_weighted_average_egf_vec(t2_real, t2_lower, t2_upper,
				       egf_vec2_lower, egf_vec2_upper, egf_vec2, n2d);
      ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
      //dumpvector(egf_vec1,0,210,"egf_vec1");
      //printf("\n%d\t%d",t2_lower, t2_upper);
      //dumpvector(egf_vec2,0,210,"egf_vec2");
      compute_weighted_average( mean_FV, egf_vec1, egf_vec2, n1d, n2d);

      free(buffer_p1_t2_lower);
      free(buffer_p2_t2_lower);
      free(buffer_p1_t2_upper);
      free(buffer_p2_t2_upper);
    }
  else//constant population
    {
      if (2*n1 > max_n2d)
	{
	  printf("ERROR: Value of 2*n1 %d exceeds max_n2d %d\n", n1, max_n2d);
	  gabort("Program terminating", 0);
	}
      file_size_bytes = compute_file_size_bytes(n1);
      //      printf("const pop: file_size_bytes %d\n", file_size_bytes); monitorinput();
      buffer_const_pop = (char*) malloc (file_size_bytes);

      read_const_pop_file_into_buffer(n1, buffer_const_pop, file_size_bytes);
      get_const_pop_egf_vec(s, n1, mean_FV, buffer_const_pop, 1);

      free(buffer_const_pop);
    }

  return (1);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
static int verbose_flag;
main(argc,argv)
int argc; char **argv;
{
  int i=0,j=0;
  int nalleles=0;
  double f0=0;
  int t=0;
  int nspikes=0;
  int sampleS=0,sampleN=0;

  double beta_alpha=0,beta_beta=0,gamma_alpha=0, gamma_beta=0,gamma2_alpha=0,gamma2_beta=0,exp_mean=0,exp_mean2=0,
    sel_coeff1=0,sel_coeff2=0,sel_coeff3=0,prob1=0,prob2=0,prob3=0;
  char *sfs_filename;

  /*
    location parameters measured
    Mean,squared mean,harmonic mean
  */
  double MEAN_S=0,MEAN_S2=0,MEAN_H=0,fix_prob=0;

  int c;

  while (1)
    {
      static struct option long_options[] =
	{
	  /* These options set a flag. */
	  // {"verbose", no_argument,       &verbose_flag, 1},
	  // {"brief",   no_argument,       &verbose_flag, 0},
	  /* These options don't set a flag.
	     We distinguish them by their indices. */
	  {"N1",     required_argument,0, 'a'},
	  {"N2",  required_argument,0, 'b'},
	  {"nalleles",  required_argument, 0, 'c'},
	  {"t",  required_argument, 0, 'd'},
	  {"f0",    required_argument, 0, 'e'},
	  {"neutral",    required_argument, 0, 'f'},
	  {"selected",    required_argument, 0, 'g'},
	  {"mode",    required_argument, 0, 'h'},
	  {"nspikes",    required_argument, 0, 'i'},
	  {"gamma_alpha",    required_argument, 0, 'j'},
	  {"gamma_beta",    required_argument, 0, 'k'},
	  {"beta_alpha",    required_argument, 0, 'l'},
	  {"beta_beta",    required_argument, 0, 'm'},
	  {"gamma2_alpha",    required_argument, 0, 'n'},
	  {"gamma2_beta",    required_argument, 0, 'o'},
	  {"exp_mean",    required_argument, 0, 'p'},
	  {"conpop",    required_argument, 0, 'q'},
	  {"file",    required_argument, 0, 'r'},
	  {"exp_mean2",    required_argument, 0, 's'},
	  {"s1",    required_argument, 0, 't'},
	  {"s2",    required_argument, 0, 'u'},
	  {"s3",    required_argument, 0, 'v'},
	  {"p1",    required_argument, 0, 'w'},
	  {"p2",    required_argument, 0, 'x'},
	  {"p3",    required_argument, 0, 'y'},
	  {0, 0, 0, 0}
	};
      /* getopt_long stores the option index here. */
      int option_index = 0;

      c = getopt_long_only (argc, argv, "",
                            long_options, &option_index);

      /* Detect the end of the options. */
      if (c == -1)
	break;

      switch (c)
	{
	case 0:
	  /* If this option set a flag, do nothing else now. */
	  if (long_options[option_index].flag != 0)
	    break;
	  printf ("option %s", long_options[option_index].name);
	  if (optarg)
	    printf (" with arg %s", optarg);
	  printf ("\nalleles");
	  break;

	case 'a':
	  N1=atoi(optarg);
	  break;
	case 'b':
	  N2=atoi(optarg);
	  break;
	case 'c':
	  nalleles=atoi(optarg);
	  break;
	case 'd':
	  t=atof(optarg);
	  break;
	case 'e':
	  f0=atof(optarg);
	  break;
	case 'f':
	  sampleS=atoi(optarg);
	  break;
	case 'g':
	  sampleN=atoi(optarg);
	  break;
	case 'h':
	  selmode=atoi(optarg);
	  break;
	case 'i':
	  nspikes=atoi(optarg);
	  break;
	case 'j':
	  gamma_alpha=atof(optarg);
	  break;
	case 'k':
	  gamma_beta=atof(optarg);
	  break;
	case 'l':
	  beta_alpha=atof(optarg);
	  break;
	case 'm':
	  beta_beta=atof(optarg);
	  break;
	case 'n':
	  gamma2_alpha=atof(optarg);
	  break;
	case 'o':
	  gamma2_beta=atof(optarg);
	  break;
	case 'p':
	  exp_mean=atof(optarg);
	  break;
	case 'q':
	  conpop=atoi(optarg);
	  break;
	case 'r':
	  sfs_filename=optarg;
	  break;
	case 's':
	  exp_mean2=atof(optarg);
	  break;
	case 't':
	  sel_coeff1=atof(optarg);
	  break;
	case 'u':
	  sel_coeff2=atof(optarg);
	  break;
	case 'v':
	  sel_coeff3=atof(optarg);
	  break;
	case 'w':
	  prob1=atof(optarg);
	  break;
	case 'x':
	  prob2=atof(optarg);
	  break;
	case 'y':
	  prob3=atof(optarg);
	  break;
	case '?':
	  /* getopt_long already printed an error message. */
	  break;

	default:
	  abort ();
	}
    }

  /* Print any remaining command line arguments (not options). */
  if (optind < argc)
    {
      printf ("non-option ARGV-elements: ");
      while (optind < argc)
	printf ("%s ", argv[optind++]);
      putchar ('\n');
    }

  double n1d=2*N1;
  double n2d=2*N2;
  /*find n_e*/
  double prop[4],n_es=0;
  double n_e=calculate_ne(N1,N2,t);

  /*Set up Spikes*/
  double * spikes_vec = (double*) calloc (nspikes+1, sizeof(double));
  double * prob_vec = (double*) calloc (nspikes+1, sizeof(double));

  for(i = 1; i <= nspikes; i++){
    printf("Give s for spike %d:\n", i);
    scanf("%lf",&spikes_vec[i]);
  }
  spikes_vec[0]=0;

  double sum_prob=0;
  for(i = 1; i <= nspikes-1; i++){
    printf("Give probability of Spike %d:\n", i);
    scanf("%lf",&prob_vec[i]);
    sum_prob+=prob_vec[i];
  }
  prob_vec[nspikes]=1-sum_prob;
  prob_vec[0]=0;

  double **FVSX= calloc(maxnd+1, sizeof(double *));

  for(i = 1; i <= nspikes; i++){
    FVSX[i] = calloc(maxnd+1,  sizeof(double));
  }
  double * FV0 = (double*) calloc (maxnd+1, sizeof(double));

  /*Set up Tables*/

  get_data_path(data_path);
  set_up_file_name(N1, s_evaluated_vec_file_const, s_evaluated_vec_file);
  set_up_file_name(N1, s_range_file_const, s_range_file);
  get_s_evaluated_vec(s_evaluated_vec, &n_s_evaluated, &n_s_evaluated_file,
		      s_evaluated_vec_file);
  get_s_ranges();

  if (!conpop)
    {
      set_up_file_name(N1, n2_evaluated_vec_file_const, n2_evaluated_vec_file);
      set_up_file_name(N1, t2_evaluated_vec_file_const, t2_evaluated_vec_file);
      set_up_file_name(N1, phase_1_dir_const, phase_1_dir);
      set_up_file_name(N1, phase_2_dir_const, phase_2_dir);

      get_int_evaluated_vec(t2_evaluated_vec,&n_t2_evaluated, &t2_lower,
			    &t2_step,&t2_evaluated_vec_file);
      get_int_evaluated_vec(n2_evaluated_vec,&n_n2_evaluated, &n2_lower,
			    &n2_step,&n2_evaluated_vec_file);
    }
  else
    {
      set_up_file_name(N1, "", const_pop_dir);
    }

  /* calculate neutral frequency vector*/
  load_FV(N1,N2,t,0.0,f0,FV0);

  /*Sampling of sites*/
  int * discrete0 = (int*) calloc (nalleles+2, sizeof(int));
  int * discrete1 = (int*) calloc (nalleles+2, sizeof(int));

  const gsl_rng_type * T;
  double selcoeff=0,uniform=0;

  gsl_rng_env_setup();
  T = gsl_rng_taus;
  gsl_rng  *rgen =gsl_rng_alloc(T);

  printf("%f,%f\n",gamma_alpha,gamma_beta);
  i=1;
  while (i<=sampleN)
    {// loop for sampling starts here
      i++;
      uniform=gsl_rng_uniform (rgen);
      double * FVS = (double*) calloc (maxnd+1, sizeof(double));
      selcoeff=0;
      switch(selmode)
	{
	case 0:case 1:
	  break;
	case 2:
	  selcoeff=gsl_ran_gamma (rgen, gamma_beta, 1/gamma_alpha);
	  break;
	case 3:
	  selcoeff=gsl_ran_beta (rgen, beta_alpha, beta_beta);
	  break;
	case 4:
	  if (uniform<=0.2) {selcoeff=gsl_ran_exponential(rgen,exp_mean);}
	  if (uniform>0.2) {selcoeff=gsl_ran_beta (rgen, beta_alpha, beta_beta);}
	  break;
	case 5:
	  if (uniform<=0.5) {selcoeff=gsl_ran_gamma (rgen, gamma_beta, 1/gamma_alpha);}
	  if (uniform>0.5) {selcoeff=gsl_ran_gamma (rgen, gamma2_beta, 1/gamma2_alpha);}
	  break;
	case 6:
	  if (uniform<=0.5) {selcoeff=gsl_ran_exponential(rgen, exp_mean);}
	  if (uniform>0.5) {selcoeff=gsl_ran_exponential (rgen, exp_mean2);}
	  break;
	case 7:
	  if (uniform<=0.5) {selcoeff=gsl_ran_gamma (rgen, gamma_beta, 1/gamma_alpha);}
	  if (uniform>0.5) {selcoeff=gsl_ran_exponential (rgen, exp_mean);}
	  break;
	case 8:
	  if (uniform<=0.8) {selcoeff=gsl_ran_exponential(rgen,exp_mean);}
	  if (uniform>0.8) {selcoeff=gsl_ran_beta (rgen, beta_alpha, beta_beta);}
	  break;
	case 9:
          selcoeff=gsl_ran_exponential(rgen,exp_mean);
	  break;
	case 10:
          if (uniform<=0.5) {selcoeff=0.0;}
          if (uniform>0.5) {selcoeff=0.05;}
	  break;
	case 11:
          if (uniform<=(0.3)) {selcoeff=0.0;}
          if (uniform>0.33&&uniform<=0.67) {selcoeff=0.05;}
          if (uniform>0.67) {selcoeff=0.5;}
	  break;
    case 12:
          if (uniform<=0.3) {selcoeff=0;}//gsl_ran_gaussian (rgen,1e-5);}
          if (uniform>0.3&&uniform<=0.4) {selcoeff=-0.05;}//gsl_ran_gaussian (rgen,1e-5);selcoeff+=0.005;}
          if (uniform>0.4&&uniform<=0.6) {selcoeff=-0.5;}//gsl_ran_gaussian (rgen,1e-5);selcoeff+=0.05;}
          if (uniform>0.6) {selcoeff=-5;} //gsl_ran_gaussian (rgen,1e-5);selcoeff+=0.5;}
	  break;
	case 13:
          selcoeff=uniform;
	  break;
    case 14:
	  if (uniform<=0.5) {selcoeff=gsl_ran_gamma (rgen, gamma_beta, 1/gamma_alpha);}
	  if (uniform>0.5) {selcoeff=gsl_ran_beta (rgen, beta_alpha, beta_beta);}
	  break;
    case 15:
	  if (uniform<=0.2) {selcoeff=gsl_ran_gamma (rgen, gamma_beta, 1/gamma_alpha);}
	  if (uniform>0.2) {selcoeff=gsl_ran_beta (rgen, beta_alpha, beta_beta);}
	  break;
    case 16:
	  if (uniform<=0.8) {selcoeff=gsl_ran_gamma (rgen, gamma_beta, 1/gamma_alpha);}
	  if (uniform>0.8) {selcoeff=gsl_ran_beta (rgen, beta_alpha, beta_beta);}
	  break;
 	case 17:
          if (uniform<=(0.3)) {selcoeff=0.0;}
          if (uniform>0.33&&uniform<=0.67) {selcoeff=0.03;}
          if (uniform>0.67) {selcoeff=0.06;}
	  break;
	case 18:
          if (uniform<=0.2) {selcoeff=sel_coeff1;}
          if (uniform>0.2&&uniform<=0.8) {selcoeff=sel_coeff2;}
          if (uniform>0.8) {selcoeff=sel_coeff3;}
	  break;
	case 19:
          if (uniform<=0.3) {selcoeff=0;}
          if (uniform>0.3) {selcoeff=0.05;}
	  break;
	case 20:
          if (uniform<=0.2) {selcoeff=0;}
          if (uniform>0.2&&uniform<=0.4) {selcoeff=0.05;}
          if (uniform>0.4) {selcoeff=0.5;}
	  break;
	case 21:
          if (uniform<=0.2) {selcoeff=0;}
          if (uniform>0.2&&uniform<=0.4) {selcoeff=0.02;}
          if (uniform>0.4&&uniform<=0.6) {selcoeff=0.05;}
          if (uniform>0.6) {selcoeff=0.1;}
          break;
    case 22:
	  selcoeff=exp_mean;
	  break;
    case 23:
	  selcoeff=gsl_ran_lognormal (rgen, gamma_alpha, gamma_beta);
	  break;
	case 24:
          if (uniform<=0.2) {selcoeff=sel_coeff1;}
          if (uniform>0.2) {selcoeff=sel_coeff2;}
	  break;
    case 25:
          if (uniform<=0.5) {selcoeff=sel_coeff1;}
          if (uniform>0.5) {selcoeff=sel_coeff2;}
	  break;
	case 26:
          if (uniform<=prob1) {selcoeff=sel_coeff1;}
          if (uniform>prob1) {selcoeff=sel_coeff2;}
	  break;
	case 27:
          if (uniform<=prob1) {selcoeff=sel_coeff1;}
          if (uniform>=prob1) {selcoeff=sel_coeff1;}
          if (uniform>prob1) {selcoeff=sel_coeff2;}
	  break;
	}//end switch

      //monitorinput();
      if(selcoeff>0){selcoeff=-selcoeff;}
      if (selcoeff<-100) {selcoeff=-100;}//free(FVS);continue;}//free(FVS);continue;
      //printf("%f\n",selcoeff);
      MEAN_S+=selcoeff;
      MEAN_S2+=-pow(selcoeff,2);
      MEAN_H+=1/selcoeff;

      if (selcoeff== 0)
        {
	  fix_prob += 0.5/n_e;
	}
      else
	{
	  fix_prob += kimura_fixation_prob(selcoeff, n_e);
	}

      /*calculate relative proportions*/
      n_es=-n_e*selcoeff;
      if (n_es<=0.1){prop[0]++;}
      if (n_es<=1.0){prop[1]++;}
      if (n_es<=10.0){prop[2]++;}
      if (n_es<100.0){prop[3]++;}

      load_FV(N1,N2,t,selcoeff,f0,FVS);

      egf_scaling(N2,f0,FV0,FVS);
      gsl_ran_discrete_t *r= gsl_ran_discrete_preproc (n2d,FVS);
      double s1=gsl_ran_discrete (rgen, r);
      double prob=(double)(s1)/n2d;

      int success=gsl_ran_binomial(rgen,prob,nalleles);
      if (success==nalleles){discrete1[0]++;}else{
	discrete1[success]++;
      }
      free(FVS);
      gsl_ran_discrete_free(r);

    }//end sampling

  MEAN_S/=sampleN;
  MEAN_S2/=sampleN;
  MEAN_H/=sampleN;
  MEAN_H=1/MEAN_H;//Reciprocal of the mean of the reciprocals
  fix_prob/=sampleN;

  fix_prob *=2*n_e;

  prop[3]=(prop[3]-prop[2])/sampleN;
  prop[2]=(prop[2]-prop[1])/sampleN;
  prop[1]=(prop[1]-prop[0])/sampleN;
  prop[0]=prop[0]/sampleN;

  egf_scaling(N2,f0,FV0,FV0);
  binomial_sampling(N2,nalleles,sampleS,FV0,discrete0);//neutral

  /*Print Output*/

  printf("%f\t%f\t%f\n",MEAN_S,MEAN_S2,MEAN_H);
  //output_sfs_to_file_thanasis_format(nalleles,discrete0,discrete1,sfs_filename);
  output_sfs_to_file_peter_format2(nalleles,sampleS,sampleN,discrete1,discrete0,sfs_filename);

  FILE *file1= fopen(strcat(sfs_filename,".mean"), "w" );
  fprintf(file1,"%0.16E %0.16E %0.16E ",MEAN_S,MEAN_S2,MEAN_H);

  for(i = 0; i <= 3; i++){
    fprintf(file1,"%0.16E ",prop[i]);
  }
  fprintf(file1,"%0.16E ",fix_prob);
  fclose(file1);
  ////////////////////////////////////////////////////////////////////////////////

  free(FV0);
  gsl_rng_free (rgen);

  free(discrete0);
  free(discrete1);

  for(i = 1; i <= nspikes; i++){
    free(FVSX[i]);
  }
  free(FVSX);

  return 0;
}


text:
/////////////////////////////////////////////////////////////////////////////////
//
//  Levenberg - Marquardt non-linear minimization algorithm
//  Modified and simplified by Ted Laurence to use for MLE of Poisson-distributed data; Used only for
//		double precision, without constraints
//  Copyright (C) 2004  Manolis Lourakis (lourakis at ics forth gr)
//  Institute of Computer Science, Foundation for Research & Technology - Hellas
//  Heraklion, Crete, Greece.
//
//  Modifications Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the Lawrence Livermore National Laboratory. Written by Ted Laurence (laurence2@llnl.gov)
//  LLNL-CODE-424602 All rights reserved.
//  This file is part of dlevmar_mle_der
//
//  Please also read Our Notice and GNU General Public License.
//
//  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.
//
/////////////////////////////////////////////////////////////////////////////////
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include <math.h>

#include "levmar_mle.h"
#include "lm_mle_compiler.h"

#define EPSILON       1E-12
#define ONE_THIRD     0.3333333334 /* 1.0/3.0 */
#define LM_REAL_MAX FLT_MAX
#define LM_REAL_MIN -FLT_MAX

double compute_chisq_measure(double *e, double *x, double *hx, int n, int fitType)
{
	int i;
	double chisq=0.0;

	switch (fitType)
	{
	case LM_CHISQ_MLE:
		for (i=0; i<n; i++)
		{
			if ( hx[i] > 0)
			{
				if (x[i]==0)
					chisq += 2 * hx[i];
				else
					chisq += 2 * (hx[i]-x[i]-x[i]*log(hx[i]/x[i]));
				e[i]=x[i]/hx[i]-1.0;
			}
			else
				e[i]=0.0;
		}
		break;
	case LM_CHISQ_NEYMAN:
		for (i=0; i<n; i++)
		{
			if (x[i]==0)
			{
				chisq += (x[i]-hx[i])*(x[i]-hx[i]);
				e[i]=x[i]-hx[i];
			}
			else
			{
				chisq += (x[i]-hx[i])*(x[i]-hx[i])/x[i];
				e[i]=(x[i]-hx[i])/x[i];
			}
		}
		break;
	case LM_CHISQ_EQUAL_WT:
		for (i=0; i<n; i++)
		{
			e[i]=x[i]-hx[i];
			chisq += e[i]*e[i];
		}
		break;
	}
	return chisq;
}

/*
 * This function seeks the parameter vector p that best describes the measurements vector x.
 * More precisely, given a vector function  func : R^m --> R^n with n>=m,
 * it finds p s.t. func(p) ~= x, i.e. the squared second order (i.e. L2) norm of
 * e=x-func(p) is minimized.
 *
 * This function requires an analytic Jacobian.
 *
 * Returns the number of iterations (>=0) if successful, LM_ERROR if failed
 *
 * For more details, see K. Madsen, H.B. Nielsen and O. Tingleff's lecture notes on
 * non-linear least squares at http://www.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
 */

int dlevmar_mle_der(
  void (*func)(double *p, double *hx, int m, int n, void *adata), /* functional relation describing measurements. A p \in R^m yields a \hat{x} \in  R^n */
  void (*jacf)(double *p, double *j, int m, int n, void *adata),  /* function to evaluate the Jacobian \part x / \part p */
  double *p,         /* I/O: initial parameter estimates. On output has the estimated solution */
  double *x,         /* I: measurement vector. NULL implies a zero vector */
  int m,              /* I: parameter vector dimension (i.e. #unknowns) */
  int n,              /* I: measurement vector dimension */
  int itmax,          /* I: maximum number of iterations */
  double opts[5],    /* I: minim. options [\mu, \epsilon1, \epsilon2, \epsilon3,\epsilon4]. Respectively the scale factor for initial \mu,
                       * stopping thresholds for ||J^T e||_inf, ||Dp||_2, chisq, and delta_chisq. Set to NULL for defaults to be used
                       */
  double info[LM_INFO_SZ],
					           /* O: information regarding the minimization. Set to NULL if don't care
                      * info[0]= ||e||_2 at initial p.
                      * info[1-4]=[ ||e||_2, ||J^T e||_inf,  ||Dp||_2, mu/max[J^T J]_ii ], all computed at estimated p.
                      * info[5]= # iterations,
                      * info[6]=reason for terminating: 1 - stopped by small gradient J^T e
                      *                                 2 - stopped by small Dp
                      *                                 3 - stopped by itmax
                      *                                 4 - singular matrix. Restart from current p with increased mu
                      *                                 5 - no further error reduction is possible. Restart with increased mu
                      *                                 6 - stopped by small ||e||_2
                      *                                 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
					  *									8 - stopped by small change in chisq on successful iteration (dF<eps4)
                      * info[7]= # function evaluations
                      * info[8]= # Jacobian evaluations
                      * info[9]= # linear systems solved, i.e. # attempts for reducing error
                      */
  double *work,     /* working memory at least LM_DER_WORKSZ() reals large, allocated if NULL */
  double *covar,    /* O: Covariance matrix corresponding to LS solution; mxm. Set to NULL if not needed. */
  void *adata,		/* pointer to possibly additional data, passed uninterpreted to func & jacf.
                      * Set to NULL if not needed
                      */
  int fitType)      /* Which type of fit to perform
					*	LM_CHISQ_MLE is MLE for Poisson distribution
					*	LM_CHISQ_NEYMAN	is least squares where sigma^2 is set to max(x,1)
					*	LM_CHISQ_EQUAL_WT is least squares where sigma^2 is set to 1 */
{
register int i, j, k, l;
int worksz, freework=0, issolved;
/* temp work arrays */
double *e,          /* nx1 */
	   *e_test,     /* nx1 */
       *hx,         /* \hat{x}_i, nx1 */
       *jacTe,      /* J^T e_i mx1 */
       *jac,        /* nxm */
       *jacTjac,    /* mxm */
	   *LUdecomp,   /* mxm */
       *Dp,         /* mx1 */
   *diag_jacTjac,   /* diagonal of J^T J, mx1 */
       *pDp;        /* p + Dp, mx1 */

register double mu,  /* damping constant */
                tmp; /* mainly used in matrix & vector multiplications */
double p_chisq, jacTe_inf, pDp_chisq; /* ||e(p)||_2, ||J^T e||_inf, ||e(p+Dp)||_2 */
double p_L2, Dp_L2=LM_REAL_MAX, dF, dL;
double tau, eps1, eps2, eps2_sq, eps3, eps4;
double init_p_chisq;
int nu=2, nu2, stop=0, nfev, njev=0, nlss=0;
const int nm=n*m;

gsl_matrix_view LUdecomp_view;
gsl_vector_view jacTe_view, Dp_view;
gsl_permutation * perm = gsl_permutation_alloc (m);
int signum;

  mu=jacTe_inf=0.0; /* -Wall */

  if(n<m){
    fprintf(stderr, "(): cannot solve a problem with fewer measurements [%d] than unknowns [%d]\n", n, m);
    return LM_ERROR;
  }

  if(!jacf){
    fprintf(stderr, "No function specified for computing the Jacobian");
    return LM_ERROR;
  }

  if(opts){
	  tau=opts[0];
	  eps1=opts[1];
	  eps2=opts[2];
	  eps2_sq=opts[2]*opts[2];
    eps3=opts[3];
	eps4=opts[4];
  }
  else{ // use default values
	  tau=LM_INIT_MU;
	  eps1=LM_STOP_THRESH;
	  eps2=LM_STOP_THRESH;
	  eps2_sq=LM_STOP_THRESH*LM_STOP_THRESH;
	  eps3=LM_STOP_THRESH;
	  eps4=LM_STOP_THRESH;
  }

  if(!work){
    worksz=LM_MLE_WORKSZ(m, n); //2*n+4*m + n*m + m*m;
    work=(double *)malloc(worksz*sizeof(double)); /* allocate a big chunk in one step */
    if(!work){
      fprintf(stderr, "(): memory allocation request failed\n");
      return LM_ERROR;
    }
    freework=1;
  }

  /* set up work arrays */
  e=work;
  e_test=e + n;
  hx=e_test + n;
  jacTe=hx + n;
  jac=jacTe + m;
  jacTjac=jac + nm;
  LUdecomp=jacTjac + m*m;
  Dp=LUdecomp + m*m;
  diag_jacTjac=Dp + m;
  pDp=diag_jacTjac + m;

  LUdecomp_view=gsl_matrix_view_array(LUdecomp,m,m);
  jacTe_view=gsl_vector_view_array(jacTe,m);
  Dp_view=gsl_vector_view_array(Dp,m);

  /* compute e=x - f(p) and its chisq measure*/
  (*func)(p, hx, m, n, adata); nfev=1;
  /* ### e=x-hx, p_eL2=||e|| */
  p_chisq=compute_chisq_measure(e, x, hx, n, fitType);

  init_p_chisq=p_chisq;
  if(!LM_FINITE(p_chisq)) stop=7;

  for(k=0; k<itmax && !stop; ++k){
	register int l, im;
	register double alpha, wt, *jaclm;

    /* Note that p and e have been updated at a previous iteration */

    if(p_chisq<=eps3){ /* error is small */
      stop=6;
      break;
    }

    /* Compute the Jacobian J at p,  J^T J,  J^T e,  ||J^T e||_inf and ||p||^2.
     * Since J^T J is symmetric, its computation can be sped up by computing
     * only its upper triangular part and copying it to the lower part
     */

    (*jacf)(p, jac, m, n, adata); ++njev;

    /* J^T J, J^T e */
	 for(i=m*m; i-->0; )
		jacTjac[i]=0.0;
	 for(i=m; i-->0; )
		jacTe[i]=0.0;

	 for(l=n; l-->0; ){
		jaclm=jac+l*m;
		switch (fitType)
		{
		case LM_CHISQ_MLE:
			if (hx[l]>0)
				wt=x[l]/(hx[l]*hx[l]);
			else
				wt=0.0;
			break;
		case LM_CHISQ_NEYMAN:
			if (x[l]==0)
				wt=1.0;
			else
				wt=1.0/x[l];
			break;
		case LM_CHISQ_EQUAL_WT:
			wt=1.0;
			break;
		}
		for(i=m; i-->0; ){
		  im=i*m;
		  alpha=jaclm[i]; //jac[l*m+i];
		  for(j=i+1; j-->0; ) /* j<=i computes lower triangular part only */
			jacTjac[im+j]+=jaclm[j]*alpha*wt; //jac[l*m+j]

		  /* J^T e */
		  jacTe[i]+=alpha*e[l];
		}
	 }

	 for(i=m; i-->0; ) /* copy to upper part */
		for(j=i+1; j<m; ++j)
		  jacTjac[i*m+j]=jacTjac[j*m+i];

	  /* Compute ||J^T e||_inf and ||p||^2 */
    for(i=0, p_L2=jacTe_inf=0.0; i<m; ++i){
      if(jacTe_inf < (tmp=FABS(jacTe[i]))) jacTe_inf=tmp;

      diag_jacTjac[i]=jacTjac[i*m+i]; /* save diagonal entries so that augmentation can be later canceled */
      p_L2+=p[i]*p[i];
    }

    /* check for convergence */
    if((jacTe_inf <= eps1)){
      Dp_L2=0.0; /* no increment for p in this case */
      stop=1;
      break;
    }

   /* compute initial damping factor */
    if(k==0){
      for(i=0, tmp=LM_REAL_MIN; i<m; ++i)
        if(diag_jacTjac[i]>tmp) tmp=diag_jacTjac[i]; /* find max diagonal element */
      mu=tau*tmp;
    }

    /* determine increment using adaptive damping */
    while(1){
      /* augment normal equations */
      for(i=0; i<m; ++i)
        jacTjac[i*m+i]+=mu;

      /* solve augmented equations */
      /* use the LU included with GSL*/
	  for (i=0; i<m*m; i++)
		LUdecomp[i]=jacTjac[i];

	  gsl_linalg_LU_decomp(&LUdecomp_view.matrix, perm, &signum);
	  issolved = !gsl_linalg_LU_solve (&LUdecomp_view.matrix, perm, &jacTe_view.vector, &Dp_view.vector);

	  if(issolved){
        /* compute p's new estimate and ||Dp||^2 */
        for(i=0, Dp_L2=0.0; i<m; ++i){
          pDp[i]=p[i] + (tmp=Dp[i]);
          Dp_L2+=tmp*tmp;
        }
        //Dp_L2=sqrt(Dp_L2);

        if(Dp_L2<=eps2_sq*p_L2){ /* relative change in p is small, stop */
          stop=2;
          break;
        }

//       if(Dp_L2>=(p_L2+eps2)/EPSILON*EPSILON){ /* almost singular */
//         stop=4;
//         break;
//      }

        (*func)(pDp, hx, m, n, adata); ++nfev; /* evaluate function at p + Dp */
        /* compute ||e(pDp)||_2 */
        /* ### hx=x-hx, pDp_chisq=||hx|| */
        pDp_chisq= compute_chisq_measure(e_test, x, hx, n, fitType);

		if(!LM_FINITE(pDp_chisq)){ /* chisq is not finite, most probably due to a user error.
                                  * This check makes sure that the inner loop does not run indefinitely.
                                  * Thanks to Steve Danauskas for reporting such cases
                                  */
          stop=7;
          break;
        }

        for(i=0, dL=0.0; i<m; ++i)
          dL+=Dp[i]*(mu*Dp[i]+jacTe[i]);

        dF=p_chisq-pDp_chisq;

        if(dL>0.0 && dF>0.0){ /* reduction in error, increment is accepted */
          tmp=(2.0*dF/dL-1.0);
          tmp=1.0-tmp*tmp*tmp;
          mu=mu*( (tmp>=ONE_THIRD)? tmp : ONE_THIRD );
          nu=2;

          for(i=0 ; i<m; ++i) /* update p's estimate */
            p[i]=pDp[i];

          for(i=0; i<n; ++i) /* update e and ||e||_2 */
            e[i]=e_test[i];
          p_chisq=pDp_chisq;
		  if (dF<eps4)
			  stop = 8;
          break;
        }
	  }

      /* if this point is reached, either the linear system could not be solved or
       * the error did not reduce; in any case, the increment must be rejected
       */

      mu*=nu;
      nu2=nu<<1; // 2*nu;
      if(nu2<=nu){ /* nu has wrapped around (overflown). Thanks to Frank Jordan for spotting this case */
        stop=5;
        break;
	  }
      nu=nu2;

      for(i=0; i<m; ++i) /* restore diagonal J^T J entries */
        jacTjac[i*m+i]=diag_jacTjac[i];
	} /* inner loop */
  }

  if(k>=itmax) stop=3;

  for(i=0; i<m; ++i) /* restore diagonal J^T J entries */
    jacTjac[i*m+i]=diag_jacTjac[i];

  if(info){
    info[0]=init_p_chisq;
    info[1]=p_chisq;
    info[2]=jacTe_inf;
    info[3]=Dp_L2;
    for(i=0, tmp=LM_REAL_MIN; i<m; ++i)
      if(tmp<jacTjac[i*m+i]) tmp=jacTjac[i*m+i];
    info[4]=mu/tmp;
    info[5]=(double)k;
    info[6]=(double)stop;
    info[7]=(double)nfev;
    info[8]=(double)njev;
    info[9]=(double)nlss;
  }

  if(freework) free(work);

  gsl_permutation_free(perm);

  return (stop!=4 && stop!=7)?  k : LM_ERROR;
}

max_stars_repo_name: david-hoffman/dphtools

max_issues_repo_name: david-hoffman/dphutils

max_forks_repo_name: david-hoffman/dphtools