text:
#pragma once

#include <string_view>
#include <gsl/span>

namespace kab_advent {
	auto day(gsl::span<std::string_view const> args) -> int;
}

ext: h


text:
#ifndef VEC_WRAPPER_H
#define VEC_WRAPPER_H

#include <atomic>

#include <gsl/gsl_vector.h>

class Vector {
public:
    Vector(); // null

    Vector(gsl_vector *wrap);
    Vector(const size_t length);
    Vector(const Vector& copy);
    virtual ~Vector();

    gsl_vector *get() const;
    void set(gsl_vector *newVector);
    void reset(gsl_vector *newVector);

    void reset(const size_t length);

    explicit operator bool() const;

    Vector &operator =(const Vector& copy);
    void operator *=(const Vector& other);
    double operator [](const size_t i) const;
    double& operator [](const size_t i);

    void coords();
    void uncoords();
    void basis(size_t i);

private:
    gsl_vector *m_vec;

    std::atomic<unsigned int> m_delCount;
};

#endif // VEC_WRAPPER_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


text:
#ifndef SNPINFER_DATASTRUCTURE
#define SNPINFER_DATASTRUCTURE

#define NUMPRECISION 1e-6

#include <iostream>
#include <fstream>
#include <sstream>
#include <string.h>
#include <vector>
#include <algorithm>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <iterator>
#include <omp.h>
#include <sys/stat.h>
using namespace std;

#define GSL_DLL
#include <gsl/gsl_sys.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_fit.h>

typedef struct SNPINFO {
        string chr;
	int posst, posed;//, label;
        //int err, info, st, ed, shift;
	string snpid;
	//double qual;
	//vector<char> alleles;
} SNPINFO;

#define MINUSINFINITE -1000000000
#define PI 3.14159265359

typedef struct MYDATA {
	int indN, totalN, L, ploidity; //indN = N*plodity, totalN = N, but if replicates > 1, then totalN > N
	int *indIndex; //index of individuals (true individuals, not haploids, length+1 for total copies, =totalN) (necessary when there are replicates per individual)
	float *data; //length = indindex[last]*ploidity*L;
	int *asz;
	vector<vector<int> > fixState; //fix states, ~ pop, indN of them
	vector<vector<bool> > fixAllele; //fix data without imputation, ~ data, totalN of them; NOTE: false is fix, true is to impute!
	vector<vector<int> > breaks; //must contain start and end position as double breaks

	//additional info
	vector<SNPINFO> snpinfo;
	vector<string> indinfo;
	vector<string> fyinfo, fxinfo;
} MYDATA;

typedef struct PARAUNIT {
	vector<float> allele;
	int lapseN, popN;
} PARAUNIT;

typedef struct TENSORPARA {
	//recombination parameters
	double priorW, probAlpha, defmut, *rPrior, *popr; //priorW is for recomb rate
							//probAlpha is for states proportion
	bool *rec, *lapse;
	float *recombN, trueindN;
	float *indWeight, *idw;

	//parameter for states
	int maxK, maxHapK; //maxHapK is used if number of states is fixed;//12/18/15: maxK is current number of cluster? and maxHapK is the maximum allowed if >0?
	int *pop;
	float *nh;
	vector<vector<PARAUNIT> > param; //LxKxZ, where K=maxK, Z=asz
	vector<vector<int> > flips;
	double *Vh, *Q, A;

	//control parameters
	int burnin, mcmc;
	bool shrinkclass, singlecall, imputedata;
	bool addsample; //graudually add samples, used for haplotype inference
	bool samplefull; //true for haplotype inference
	char changeAlpha; //add (>0) or reduce (<0) probAlpha at the begining, add used for haplotype inference, reduce used for stratification mapping
	int splitstop; //# of pops, beyond which no more splitting
	int refweight; //weight of reference data with known states
	char maximum; //sample or maximize, 3bits abc: a: maximum recomb, b: maximum state, c: maximum allele
	bool samplemaximum; //indicating if first sample and then take maximum, used for pop inference
	bool logscale; //indicating if probability is in local scale, used for popinfer
	double recrate; //recombination rate to be used for set up recomb priors
	int recombcut; //above which position specific recomb is counted
	double heteroVh; //0 if Vh is global, > 0 if locally Vh changes. Only implemented for ploidity = 1, and not implemented in greedyswitch function
} TENSORPARA;

typedef struct TENSORPARAP0 {
	float *recombN0;
	float *nh0; //*nh0pp;
	double *Q0;
	unsigned short *param0;//, *param0pp;
	int clustersz0, popn0, unitsz0, totalN0;
} TENSORPARA0;

inline float fast_log2(float val)
{
   int * const    exp_ptr = reinterpret_cast <int *> (&val);
   int            x = *exp_ptr;
   const int      log_2 = ((x >> 23) & 255) - 128;
   x &= ~(255 << 23);
   x += 127 << 23;
   *exp_ptr = x;

   val = ((-1.0f/3) * val + 2) * val - 2.0f/3;   // (1)

   return (val + log_2);
}

inline float fast_log(const float &val)
{
   return (fast_log2 (val) * 0.69314718f);
}

//void call_IDEAS(int, char *argv[]);

#endif


text:
//im_clam.c --
// A. Kern 3/20/15
//
// Composite Likelihood estimation of canonical IM models via AFS

//going to use NLopt for optimization

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <nlopt.h>
#include "AFS.h"
#include "adkGSL.h"
#include "cs.h"
#include "time.h"
#include <unistd.h>
#include "adkCSparse.h"
#include "AFS_ctmc.h"
#include <slepcmfn.h>
#include "AFS_ctmc_petsc.h"
#include "im_clam.h"

void import2DSFSData(const char *fileName, gsl_matrix *obsData);

PetscInt n1, n2;

afsStateSpace *stateSpace, *reducedStateSpace;

clam_lik_params *currentParams, *nextParams;

gsl_matrix *transMat;
double *res;

double lowerBounds[5] = {0.01,0.01,0.0,0.0,0.001};
double upperBounds[5] = {10.0,10.0,20.0,20.0,10.0};
PetscBool vbse = PETSC_FALSE;

static char help[] = "im_clam\n\
	Example: mpiexec -n <np> ./im_clam -s <stateSpace file> -m <mats file> -d <data file> \n\n\toptions:\n\
	\t-exp expected value mode (requires -x flag too)\n\
	\t-GIM uncertainty estimation mode (requires -x flag too)\n\
	\t-mo multiple optimizations from different start points\n\
	\t-global multi-level optimization (MLSL algo.)\n\
	\t-x <theta_2, theta_A, mig12, mig21, t_div> parameter starting values\n\
	\t-obs (prints out observed AFS as well as that expected from MLE params)\n\
	\t-u mutation rate per base pair per generation (only used to unscale parameters; default 1e-8)\n\
	\t-g generation time (gens/year; default 20)\n\
	\t-seqLen sequence length scanned for polymorphisms (used to unscale parameter)\n\
	\t-put upper bound for optimization of thetas\n\
	\t-pum upper bound for optimization of migration rates\n\
	\t-pudt upper bound for optimization of divergence time\n\
	\t-r randomSeed\n\
	\t-v verbose output\n";

int main(int argc, char **argv){
	PetscInt i,j, N,runMode;
	clock_t time1, time2;
	PetscInt nnz;
	PetscInt seed;
	char  filename[PETSC_MAX_PATH_LEN], filename2[PETSC_MAX_PATH_LEN], filename3[PETSC_MAX_PATH_LEN] ;
	double lik = 0.0;
	double mle[5] = {0.1,2,1,1,5};
 	const gsl_rng_type * T;
  	gsl_rng * r;
	PetscErrorCode ierr;
	PetscMPIInt    rank,size;
	PetscBool      flg,obsFlag;
//	Vec tmpVec;
	PetscInt dim=5;
	double snpNumber, pi_est, p, N0;
	double u=1e-8;
	double genPerYear=20;
	double put = 10.0;
	double pum = 20.0;
	double pudt=10.0;
	double propSnp;
	PetscInt seqLen;
	gsl_matrix *fi,*gi;
	FILE *infile;
	PetscInt testInt=0;
	PetscScalar one = 1.0;
	cs *ident;

	///////////////////
	/////	PETSC / Slepc library version of this code
	SlepcInitialize(&argc,&argv,(char*)0,help);
	ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
	ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
	if(rank==0){

		printf("\n.___   _____             .__\n|   | /     \\       ____ |  | _____    _____\n|   |/  \\ /  \\    _/ ___\\|  | \\__  \\  /     \\\n|   /    Y    \\   \\  \\___|  |__/ __ \\|  Y Y  \\\n|___\\____|__  /____\\___  >____(____  /__|_|  /\n            \\/_____/   \\/          \\/      \\/\n\n\n");

		printf("im_clam -- Isolation with Migration Composite Likelihood Analysis using Markov chains\n");
		//printf("A.D. Kern 2015\n///////////////////\n");
		printf("\n\n");
	}
	if(argc<2){
		printf("%s",help);
		exit(666);
	}
	time1=clock();

	ierr = PetscOptionsGetString(NULL,"-s",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetString(NULL,"-m",filename2,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetString(NULL,"-d",filename3,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetInt(NULL,"-r",&seed,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetReal(NULL,"-u",&u,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetReal(NULL,"-g",&genPerYear,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetInt(NULL,"-seqLen",&seqLen,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetReal(NULL,"-put",&put,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetReal(NULL,"-pum",&pum,&flg);CHKERRQ(ierr);
	ierr = PetscOptionsGetReal(NULL,"-pudt",&pudt,&flg);CHKERRQ(ierr);
	//set upper bounds
	upperBounds[0]=put;upperBounds[1]=put;upperBounds[2]=pum;upperBounds[3]=pum;upperBounds[4]=pudt;

	if(!flg) seed=time(NULL);

	//printf("rank:%d %d\n",rank,seed);
	MPI_Bcast(&seed,1,MPI_INT,0,PETSC_COMM_WORLD);
	//printf("rank:%d %d\n",rank,seed);

	runMode = 1;
	obsFlag=PETSC_FALSE;

	//setup RNG for starting point for optimization
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc (T);
	gsl_rng_set(r,seed);

	for(i=0;i<dim;i++){
		mle[i] = gsl_ran_flat(r,lowerBounds[i], upperBounds[i]);
	}

	PetscOptionsHasName(NULL,"-v",&flg); if(flg) vbse=PETSC_TRUE;
	PetscOptionsHasName(NULL,"-obs",&flg); if(flg) obsFlag=PETSC_TRUE;

	PetscOptionsHasName(NULL,"-exp",&flg); if(flg) runMode=2;
	ierr = PetscOptionsGetRealArray(NULL,"-x",mle,&dim,&flg);CHKERRQ(ierr);

	//double check we are set for exp mode
	if(runMode == 2 && !flg){
		printf("for exp runmode need to set parameter values using -x flag\n");
		exit(1);
	}

	PetscOptionsHasName(NULL,"-mo",&flg); if(flg) runMode=3;
	PetscOptionsHasName(NULL,"-global",&flg); if(flg) runMode=4;
	PetscOptionsHasName(NULL,"-GIM",&flg); if(flg) runMode=5;

	//import state space
	stateSpace = afsStateSpaceImportFromFile(filename);
	N = stateSpace->nstates;

	//quick peak at the mats file to set nnz
	//open file
	infile = fopen(filename2, "r");
	if (infile == NULL){
		fprintf(stderr,"Error opening mats file! ARRRRR!!!!\n");
		exit(1);
	}
	testInt=fscanf(infile,"nnz: %d", &nnz);
	fclose(infile);
	//got nnz, time for allocs

	//setup currentParams struct; alloc all arrays
	currentParams = malloc(sizeof(clam_lik_params));
	currentParams->rng = r;
	currentParams->n1 = n1 = stateSpace->states[0]->popMats[0]->size1 - 1;
	currentParams->n2 = n2 = stateSpace->states[0]->popMats[1]->size2 - 1;
	currentParams->stateSpace = stateSpace;
	currentParams->rank = rank;
//	currentParams->snpNumber = snpNumber;
	currentParams->expAFS = gsl_matrix_alloc(n1+1,n2+1);
	currentParams->expAFS2 = gsl_matrix_alloc(n1+1,n2+1);
	currentParams->map = malloc(sizeof(PetscInt)*N);
	currentParams->reverseMap = malloc(sizeof(PetscInt)*N);
	currentParams->rates = gsl_vector_alloc(stateSpace->nstates);
	currentParams->stateVec = gsl_vector_alloc(stateSpace->nstates);
	currentParams->resVec = gsl_vector_alloc(stateSpace->nstates);
	gsl_vector_set_zero(currentParams->resVec);
//	currentParams->paramVector = gsl_vector_alloc(5);
//	currentParams->paramCIUpper = gsl_vector_alloc(5);
//	currentParams->paramCILower = gsl_vector_alloc(5);
//	currentParams->mlParams = gsl_vector_alloc(5);

	currentParams->top = malloc(sizeof(double) * nnz);
	currentParams->topA = malloc(sizeof(double) * nnz);

	currentParams->move = malloc(sizeof(PetscInt) * nnz);
	currentParams->moveA = malloc(sizeof(PetscInt) * nnz);
//	printf("allocated moveType...\n");
	for(i=0;i< (nnz);i++){
		currentParams->move[i]=0;
		currentParams->top[i]=0;
		currentParams->moveA[i]=0;
		currentParams->topA[i]=0;
	}
	currentParams->dim1 = malloc(sizeof(PetscInt) * nnz);
	currentParams->dim2 = malloc(sizeof(PetscInt) * nnz);
	currentParams->dim1A = malloc(sizeof(PetscInt) * N * 10);
	currentParams->dim2A = malloc(sizeof(PetscInt) * N * 10);
	currentParams->b = malloc(N*sizeof(double));
	currentParams->expoArray = malloc(N*sizeof(double));
	currentParams->st = malloc(N*sizeof(double));
	currentParams->paramVector = gsl_vector_alloc(5);
	gsl_vector_set_all(currentParams->paramVector,1.0);
	currentParams->fEvals=0;

	currentParams->triplet = cs_spalloc(N, N, nnz + N , 1, 1); //alloc sparse mat with extra space for nonzero identity mats
	ident = cs_spalloc(N,N,N,1,1);
	for(i=0;i<N;i++) cs_entry(ident,i,i,1);
	currentParams->eye = cs_compress(ident);
	cs_spfree(ident);

	//set up some petsc matrices//////////////////////////////////////////////////////////////
	//
	//

	ierr = MatCreate(PETSC_COMM_WORLD,&currentParams->C);CHKERRQ(ierr);
//	MatSetType(C,MATMPIAIJ);
	ierr = MatSetSizes(currentParams->C, PETSC_DECIDE, PETSC_DECIDE,N,N);CHKERRQ(ierr);
	ierr = MatSetFromOptions(currentParams->C);CHKERRQ(ierr);
//	MatSetOption(currentParams->C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
	ierr = MatSetUp(currentParams->C);CHKERRQ(ierr);

	ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE, PETSC_DECIDE, N, N, NULL, &currentParams->denseMat1);CHKERRQ(ierr);
	ierr = MatSetFromOptions(currentParams->denseMat1);CHKERRQ(ierr);
	MatAssemblyBegin(currentParams->denseMat1,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(currentParams->denseMat1,MAT_FINAL_ASSEMBLY);

	///////// Set up MFN ////////////////////////////////////////////////////////////////////////////////////////
	ierr = MFNCreate(PETSC_COMM_WORLD,&currentParams->mfn);CHKERRQ(ierr);
	ierr = MFNSetFunction(currentParams->mfn,SLEPC_FUNCTION_EXP);CHKERRQ(ierr);
	ierr = MFNSetTolerances(currentParams->mfn,1e-07,PETSC_DEFAULT);CHKERRQ(ierr);

	//mapping stateSpace to reducedSpace
	currentParams->reducedStateSpace = afsStateSpaceNew();
	afsStateSpaceMapAndReducePopn(currentParams->stateSpace, currentParams->map, currentParams->reducedStateSpace, currentParams->reverseMap);
	currentParams->Na=currentParams->reducedStateSpace->nstates;

	VecCreate(PETSC_COMM_WORLD,&currentParams->ancStateVec);
	VecSetSizes(currentParams->ancStateVec,PETSC_DECIDE,currentParams->Na);
	VecSetFromOptions(currentParams->ancStateVec);
	VecDuplicate(currentParams->ancStateVec,&currentParams->ancResVec);

	//set up more matrices
	ierr = MatCreate(PETSC_COMM_WORLD,&currentParams->C2);CHKERRQ(ierr);
	MatSetType(currentParams->C2,MATMPIAIJ);
	ierr = MatSetSizes(currentParams->C2, PETSC_DECIDE, PETSC_DECIDE,currentParams->Na,currentParams->Na);CHKERRQ(ierr);
	ierr = MatSetFromOptions(currentParams->C2);CHKERRQ(ierr);
	ierr = MatSetUp(currentParams->C2);CHKERRQ(ierr);

	ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE, PETSC_DECIDE, currentParams->Na, currentParams->Na, NULL, &currentParams->denseMat2);CHKERRQ(ierr);

	ident = cs_spalloc(currentParams->Na,currentParams->Na,currentParams->Na,1,1);
	for(i=0;i<currentParams->Na;i++) cs_entry(ident,i,i,1);
	currentParams->eyeAnc = cs_compress(ident);
	cs_spfree(ident);

	//DTMC with PETSC  /////////////
//	ierr = MatCreate(PETSC_COMM_WORLD,&currentParams->D);CHKERRQ(ierr);
//	ierr = MatSetSizes(currentParams->D, PETSC_DECIDE, PETSC_DECIDE,N,N);CHKERRQ(ierr);
//	ierr = MatSetFromOptions(currentParams->D);CHKERRQ(ierr);
//	ierr = MatSetUp(currentParams->D);CHKERRQ(ierr);

//	ierr = MatCreate(PETSC_COMM_WORLD,&currentParams->D_copy);CHKERRQ(ierr);
//	ierr = MatSetSizes(currentParams->D_copy, PETSC_DECIDE, PETSC_DECIDE,N,N);CHKERRQ(ierr);
//	ierr = MatSetFromOptions(currentParams->D_copy);CHKERRQ(ierr);
//	ierr = MatSetUp(currentParams->D_copy);CHKERRQ(ierr);

	//sparse identity mat
//	ierr = MatCreate(PETSC_COMM_WORLD,&currentParams->ident);CHKERRQ(ierr);
//	ierr = MatSetSizes(currentParams->ident, PETSC_DECIDE, PETSC_DECIDE,N,N);CHKERRQ(ierr);
//	ierr = MatSetFromOptions(currentParams->ident);CHKERRQ(ierr);
//	ierr = MatSetUp(currentParams->ident);CHKERRQ(ierr);

  //	MatAssemblyBegin(currentParams->ident,MAT_FINAL_ASSEMBLY);
  //	MatAssemblyEnd(currentParams->ident,MAT_FINAL_ASSEMBLY);
//	MatShift(currentParams->ident,one);

	//dense identity mat
//	ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE, PETSC_DECIDE, N, N, NULL, &currentParams->denseIdent);CHKERRQ(ierr);
//	ierr = MatSetFromOptions(currentParams->denseIdent);CHKERRQ(ierr);

//	MatAssemblyBegin(currentParams->denseIdent,MAT_FINAL_ASSEMBLY);
  //	MatAssemblyEnd(currentParams->denseIdent,MAT_FINAL_ASSEMBLY);
//	MatShift(currentParams->denseIdent,one);

//	VecCreate(PETSC_COMM_WORLD,&currentParams->xInv);
//	VecSetSizes(currentParams->xInv,PETSC_DECIDE,N);
//	VecSetFromOptions(currentParams->xInv);
//	VecDuplicate(currentParams->xInv,&currentParams->bInv);

	///////////////////////////////////////

	//////////////////////////////

	// setup done! ///////////////
	////////////////////////////////////

	//import mats
	mcMatsImportFromFile(filename2, &nnz, currentParams->top, currentParams->move, currentParams->dim1, currentParams->dim2);
	currentParams->nnz = nnz;

	//init snpNumber to avoid compiler warning
	snpNumber=0;
	if(runMode != 2){
		//alloc data matrix and import data; estimate pi and N1
		currentParams->obsData = gsl_matrix_alloc(n1+1,n2+1);
		import2DSFSData(filename3, currentParams->obsData);
		snpNumber = gsl_matrix_sum(currentParams->obsData);
		propSnp = (float) snpNumber / (float) seqLen;
	//	pi_est = 0.0;
	//	for(i=1;i<n1;i++){
	//		for(j=0;j<n2+1;j++){
	//			p =  (float) i / n1;
	//			pi_est += (2 * p * (1.0 - p)) * gsl_matrix_get(currentParams->obsData,i,j);
	//		}
	//	}
	//	pi_est *= n1 / (n1 - 1) / snpNumber;
	//	N0=pi_est / u / 4;
	}

	////////////////////////////
	time2=clock();
	if(rank==0)printf("setup time:%f secs\n\n",(double) (time2-time1)/CLOCKS_PER_SEC);

	time1=clock();
	//expected AFS
	//printf("runMode = %d\n",runMode);
	switch(runMode){

		case 1:
		if(rank==0){
			printf("\nParameter estimation run mode\n\n");
			printf("now optimizing....\n\n");
			printf("initial parameter guess:\n");
			for(i=0;i<5;i++)printf("%f ",mle[i]);
			printf("\n\n");
		}
		maximizeLikNLOpt(&lik, currentParams, mle);
		fi = getGodambeInfoMatrix(mle, lik, currentParams);
		//get N0 estimate
		N0 = propSnp / currentParams->meanTreeLength / 4.0 / u;
		if(rank == 0){
			printf("for scaling:\nu: %e gen: %lf N0:%lf meanTreeLength:%lf seqLen:%d\n",u,genPerYear,N0,currentParams->meanTreeLength,seqLen);
			printf("Composite Likelihood estimates of IM params (scaled by 1/theta_pop1):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<5;i++)printf("%f\t",(float)mle[i]);
			printf("\n\nComposite Likelihood estimates of IM params (unscaled):\n");

			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<2;i++)printf("%f\t",(float)mle[i]*N0);
			for(i=2;i<4;i++)printf("%f\t",(float)mle[i]/N0/4);
			printf("%f\t",(float)mle[i] * N0 * 4.0/ (float)genPerYear);
			printf("\n\nUncertainty estimates of IM params (scaled by 1/theta_pop1):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<5;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(fi,i,i)));
			printf("\n\nUncertainty estimates of IM params (unscaled):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<2;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(fi,i,i))*N0);
			for(i=2;i<4;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(fi,i,i))/N0/4);
			printf("%f\t",(float)sqrt(gsl_matrix_get(fi,i,i)) * N0 * 4.0/ (float)genPerYear);
			printf("\n");

			printf("\n\nlikelihood: %lf\n",-lik);
			currentParams->nnz = nnz;
			printf("\nExpected AFS:\n");
			gsl_matrix_prettyPrint(currentParams->expAFS);
		}

		break;
		case 2:
		if(rank == 0)printf("expected value run mode\n");
		for(i=0;i<5;i++){
			gsl_vector_set(currentParams->paramVector,i,mle[i]);
		}
		calcLogAFS_IM(currentParams);
		currentParams->nnz = nnz;
		if(rank == 0){
			printf("Expected AFS:\n");
			gsl_matrix_prettyPrint(currentParams->expAFS);
			printf("parameter values used:\n");
			for(i=0;i<dim;i++)printf("%f\t",gsl_vector_get(currentParams->paramVector,i));
			printf("\n\n");
		}
		time2=clock();
		//lik = calcLikNLOpt(5,mle,NULL,currentParams);
		//

		//printf("Likelihood: %g\n",lik);
		//if(rank==0)printf("w/ CSPARSE time:%f secs\n Liklihood Func. Evals: %d\n",(double) (time2-time1)/CLOCKS_PER_SEC,currentParams->fEvals);
		time1=clock();
		break;
		case 3:
		if(rank==0){
			printf("\nParameter estimation run mode with multiple optimizations\n\n");
			printf("now optimizing....\n\n");
			printf("initial parameter guess:\n");
			for(i=0;i<5;i++)printf("%f ",mle[i]);
			printf("\n\n");
		}
		for(j = 0; j < 3; j++){

			for(i=0;i<5;i++)mle[i] = gsl_ran_flat(r,lowerBounds[i], upperBounds[i]);
			printf("optimization %d initial parameter guess:\n",j);
			for(i=0;i<5;i++)printf("%f ",mle[i]);
			maximizeLikNLOpt(&lik, currentParams, mle);
			if(rank == 0){
				printf("Composite Likelihood estimates of IM params (scaled by 1/theta_pop1):\n");
				printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
				for(i=0;i<5;i++)printf("%f\t",(float)mle[i]);
				printf("\n\nlikelihood: %lf\n",-lik);
				currentParams->nnz = nnz;
				printf("\nExpected AFS:\n");
				gsl_matrix_prettyPrint(currentParams->expAFS);
			}
		}

		break;
		case 4:
		if(rank==0){
			printf("\nParameter estimation run mode\n\n");
			printf("now optimizing....\n\n");
			printf("initial parameter guess:\n");
			for(i=0;i<5;i++)printf("%f ",mle[i]);
			printf("\n\n");
		}
		maximizeLikNLOpt_MLSL(&lik, currentParams, mle);
		fi = getGodambeInfoMatrix(mle, lik, currentParams);
		if(rank == 0){
			N0 = propSnp / currentParams->meanTreeLength / 4.0 / u;
			printf("Composite Likelihood estimates of IM params (scaled by 1/theta_pop1):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<5;i++)printf("%f\t",(float)mle[i]);
			printf("\n\nComposite Likelihood estimates of IM params (unscaled):\n");

			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<2;i++)printf("%f\t",(float)mle[i]*N0);
			for(i=2;i<4;i++)printf("%f\t",(float)mle[i]/N0/4);
			printf("%f\t",(float)mle[i] * N0 * 4.0/ (float)genPerYear);
			printf("\n\nUncertainty estimates of IM params (scaled by 1/theta_pop1):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<5;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(fi,i,i)));
			printf("\n\nUncertainty estimates of IM params (unscaled):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<2;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(fi,i,i))*N0);
			for(i=2;i<4;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(fi,i,i))/N0/4);
			printf("%f\t",(float)sqrt(gsl_matrix_get(fi,i,i)) * N0 * 4.0/ (float)genPerYear);
			printf("\n");

			printf("\n\nlikelihood: %lf\n",-lik);
			currentParams->nnz = nnz;
			printf("\nExpected AFS:\n");
			gsl_matrix_prettyPrint(currentParams->expAFS);
		}

		break;
		case 5:
		for(i=0;i<5;i++){
			gsl_vector_set(currentParams->paramVector,i,mle[i]);
		}
		if(rank==0){
			printf("\nUncertainty estimation (via Godambe Information) run mode\n\n");
			printf("MLE parameters:\n");
			for(i=0;i<5;i++)printf("%f ",mle[i]);
			printf("\n");
		}
		lik = calcLikNLOpt(5,mle,NULL,currentParams);
		gi = getGodambeInfoMatrix(mle, lik, currentParams);
		N0 = propSnp / currentParams->meanTreeLength / 4.0 / u;

		// if(rank==0){
		// 	printf("FIM:\n");
		// 	gsl_matrix_prettyPrint(fi);
		// 	printf("\nGIM:\n");
		// 	gsl_matrix_prettyPrint(gi);
		// 	printf("\n");
		// }
		if(rank==0){
			printf("\nlikelihood: %lf\n",-lik);
			printf("Composite Likelihood estimates of IM params (scaled by 1/theta_pop1):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<5;i++)printf("%f\t",(float)mle[i]);
			printf("\n\nComposite Likelihood estimates of IM params (unscaled):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<4;i++)printf("%f\t",(float)mle[i]*N0);
			printf("%f\t",(float)mle[i] * N0 * 4.0/ (float)genPerYear);
			printf("\n\nUncertainty estimates of IM params (scaled by 1/theta_pop1):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<5;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(gi,i,i)));
			printf("\n\nUncertainty estimates of IM params (unscaled):\n");
			printf("theta_pop2\ttheta_anc\tmig_1->2\tmig_2->1\tt_div\n");
			for(i=0;i<4;i++)printf("%f\t",(float) sqrt(gsl_matrix_get(gi,i,i))*N0);
			printf("%f\t",(float)sqrt(gsl_matrix_get(gi,i,i)) * N0 * 4.0/ (float)genPerYear);
			printf("\n");
		}
		break;
	}

	if(obsFlag && rank==0){
		gsl_matrix_scale(currentParams->obsData,1.0/snpNumber);
		printf("%f SNPs in dataset\n",snpNumber);
		printf("Observed AFS:\n");
		gsl_matrix_prettyPrint(currentParams->obsData);
		gsl_matrix_scale(currentParams->obsData,snpNumber);
	}

	MatDestroy(&currentParams->denseMat1);
	MatDestroy(&currentParams->denseMat2);
	MatDestroy(&currentParams->C);
	MatDestroy(&currentParams->C2);
	MatDestroy(&currentParams->D);
	MatDestroy(&currentParams->D_copy);
	MatDestroy(&currentParams->denseIdent);
	MatDestroy(&currentParams->ident);
	MFNDestroy(&currentParams->mfn);
	VecDestroy(&currentParams->ancStateVec);
	VecDestroy(&currentParams->ancResVec);
	VecDestroy(&currentParams->xInv);
	VecDestroy(&currentParams->bInv);
	gsl_rng_free (r);

	afsStateSpaceFree(currentParams->reducedStateSpace);
	time2=clock();
	if(rank==0)printf("total run time:%f secs\n Liklihood Func. Evals: %d\n",(double) (time2-time1)/CLOCKS_PER_SEC,currentParams->fEvals);
	ierr = PetscFinalize();

	return(0);
}

//import2DSFSData -- imports a matrix from fileName and stuffs it in prealloc'd obsData
void import2DSFSData(const char *fileName, gsl_matrix *obsData){
	FILE *infile;

	//open file
	infile = fopen(fileName, "r");
	if (infile == NULL){
		fprintf(stderr,"Error opening data file! ARRRRR!!!!\n");
		exit(1);
	}
	gsl_matrix_fscanf(infile,obsData);
	fclose(infile);
}

ext: c

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
/*
 * Copyright (C) 2017 Jelmer Ypma. All Rights Reserved.
 * This code is published under the L-GPL.
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 * File:   nloptrAPI.h
 * Author: Jelmer Ypma
 * Date:   3 October 2017
 *
 * This file provides an API for calling internal NLopt code from C within
 * R packages. The C functions that are registered in init_nloptr.c can be
 * accessed by external R packages.
 *
 * 03/10/2017: Initial version exposing nlopt_version.
 */

#ifndef __NLOPTRAPI_H__
#define __NLOPTRAPI_H__

#include <R_ext/Rdynload.h>
#include <R.h>
#include <Rinternals.h>

#include <nlopt.h>

/*
 * C functions can be exposed using the following template:
 *
 * RET_TYPE FUNCNAME(ARGTYPE_1 ARGNAME 1, ARGTYPE_2 ARGNAME_2)
 * {
 *     static RET_TYPE(*fun)(ARGTYPE_1, ARGTYPE_2) = NULL;
 *     if (fun == NULL) fun = (RET_TYPE(*)(ARGTYPE_1, ARGTYPE_2)) R_GetCCallable("nloptr","FUNCNAME");
 *     return fun(ARGNAME_1, ARGNAME_2);
 * }
 *
 */

inline NLOPT_EXTERN(const char *) nlopt_algorithm_name(nlopt_algorithm a)
{
    static const char *(*fun)(nlopt_algorithm) = NULL;
    if (fun == NULL) fun = (const char *(*)(nlopt_algorithm)) R_GetCCallable("nloptr","nlopt_algorithm_name");
    return fun(a);
}

inline NLOPT_EXTERN(void) nlopt_srand(unsigned long seed)
{
    static void(*fun)(unsigned long) = NULL;
    if (fun == NULL) fun = (void(*)(unsigned long)) R_GetCCallable("nloptr","nlopt_srand");
    return fun(seed);
}

inline NLOPT_EXTERN(void) nlopt_srand_time(void)
{
    static void(*fun)(void) = NULL;
    if (fun == NULL) fun = (void(*)(void)) R_GetCCallable("nloptr","nlopt_srand_time");
    return fun();
}

inline NLOPT_EXTERN(void) nlopt_version(int *major, int *minor, int *bugfix)
{
    static void(*fun)(int *, int *, int *) = NULL;
    if (fun == NULL) fun = (void(*)(int *, int *, int *)) R_GetCCallable("nloptr","nlopt_version");
    return fun(major, minor, bugfix);
}

inline NLOPT_EXTERN(nlopt_opt) nlopt_create(nlopt_algorithm algorithm, unsigned n)
{
    static nlopt_opt(*fun)(nlopt_algorithm, unsigned) = NULL;
    if (fun == NULL) fun = (nlopt_opt(*)(nlopt_algorithm, unsigned)) R_GetCCallable("nloptr","nlopt_create");
    return fun(algorithm, n);
}

inline NLOPT_EXTERN(void) nlopt_destroy(nlopt_opt opt)
{
    static void(*fun)(nlopt_opt) = NULL;
    if (fun == NULL) fun = (void(*)(nlopt_opt)) R_GetCCallable("nloptr","nlopt_destroy");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_opt) nlopt_copy(const nlopt_opt opt)
{
    static nlopt_opt(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (nlopt_opt(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_copy");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_optimize(nlopt_opt opt, double *x, double *opt_f)
{
    static nlopt_result(*fun)(nlopt_opt, double *, double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double *, double *)) R_GetCCallable("nloptr","nlopt_optimize");
    return fun(opt, x, opt_f);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, void *f_data)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, void *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, void *)) R_GetCCallable("nloptr","nlopt_set_min_objective");
    return fun(opt, f, f_data);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_max_objective(nlopt_opt opt, nlopt_func f, void *f_data)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, void *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, void *)) R_GetCCallable("nloptr","nlopt_set_max_objective");
    return fun(opt, f, f_data);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_precond_min_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, nlopt_precond, void *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, nlopt_precond, void *)) R_GetCCallable("nloptr","nlopt_set_precond_min_objective");
    return fun(opt, f, pre, f_data);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_precond_max_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, nlopt_precond, void *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, nlopt_precond, void *)) R_GetCCallable("nloptr","nlopt_set_precond_max_objective");
    return fun(opt, f, pre, f_data);
}

inline NLOPT_EXTERN(nlopt_algorithm) nlopt_get_algorithm(const nlopt_opt opt)
{
    static nlopt_algorithm(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (nlopt_algorithm(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_algorithm");
    return fun(opt);
}

inline NLOPT_EXTERN(unsigned) nlopt_get_dimension(const nlopt_opt opt)
{
    static unsigned(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (unsigned(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_dimension");
    return fun(opt);
}

/* constraints: */

inline NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bounds(nlopt_opt opt, const double *lb)
{
    static nlopt_result(*fun)(nlopt_opt, const double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, const double *)) R_GetCCallable("nloptr","nlopt_set_lower_bounds");
    return fun(opt, lb);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bounds1(nlopt_opt opt, double lb)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_lower_bounds1");
    return fun(opt, lb);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_get_lower_bounds(const nlopt_opt opt, double *lb)
{
    static nlopt_result(*fun)(const nlopt_opt, double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(const nlopt_opt, double *)) R_GetCCallable("nloptr","nlopt_get_lower_bounds");
    return fun(opt, lb);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bounds(nlopt_opt opt, const double *ub)
{
    static nlopt_result(*fun)(nlopt_opt, const double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, const double *)) R_GetCCallable("nloptr","nlopt_set_upper_bounds");
    return fun(opt, ub);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bounds1(nlopt_opt opt, double ub)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_upper_bounds1");
    return fun(opt, ub);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_get_upper_bounds(const nlopt_opt opt, double *ub)
{
    static nlopt_result(*fun)(const nlopt_opt, double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(const nlopt_opt, double *)) R_GetCCallable("nloptr","nlopt_get_upper_bounds");
    return fun(opt, ub);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_remove_inequality_constraints(nlopt_opt opt)
{
    static nlopt_result(*fun)(nlopt_opt) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt)) R_GetCCallable("nloptr","nlopt_remove_inequality_constraints");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_constraint(nlopt_opt opt,
             nlopt_func fc,
             void *fc_data,
             double tol)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, void *, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, void *, double)) R_GetCCallable("nloptr","nlopt_add_inequality_constraint");
    return fun(opt, fc, fc_data, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_add_precond_inequality_constraint(
        nlopt_opt opt, nlopt_func fc, nlopt_precond pre, void *fc_data,
        double tol)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, nlopt_precond, void *, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, nlopt_precond, void *, double)) R_GetCCallable("nloptr","nlopt_add_precond_inequality_constraint");
    return fun(opt, fc, pre, fc_data, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_mconstraint(nlopt_opt opt,
             unsigned m,
             nlopt_mfunc fc,
             void *fc_data,
             const double *tol)
{
    static nlopt_result(*fun)(nlopt_opt, unsigned, nlopt_mfunc, void *, const double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, unsigned, nlopt_mfunc, void *, const double *)) R_GetCCallable("nloptr","nlopt_add_inequality_mconstraint");
    return fun(opt, m, fc, fc_data, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_remove_equality_constraints(nlopt_opt opt)
{
    static nlopt_result(*fun)(nlopt_opt) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt)) R_GetCCallable("nloptr","nlopt_remove_equality_constraints");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_add_equality_constraint(nlopt_opt opt,
             nlopt_func h,
             void *h_data,
             double tol)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, void *, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, void *, double)) R_GetCCallable("nloptr","nlopt_add_equality_constraint");
    return fun(opt, h, h_data, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_add_precond_equality_constraint(
        nlopt_opt opt, nlopt_func h, nlopt_precond pre, void *h_data,
        double tol)
{
    static nlopt_result(*fun)(nlopt_opt, nlopt_func, nlopt_precond, void *, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, nlopt_func, nlopt_precond, void *, double)) R_GetCCallable("nloptr","nlopt_add_precond_equality_constraint");
    return fun(opt, h, pre, h_data, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_add_equality_mconstraint(nlopt_opt opt,
             unsigned m,
             nlopt_mfunc h,
             void *h_data,
             const double *tol)
{
    static nlopt_result(*fun)(nlopt_opt, unsigned, nlopt_mfunc, void *, const double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, unsigned, nlopt_mfunc, void *, const double *)) R_GetCCallable("nloptr","nlopt_add_equality_mconstraint");
    return fun(opt, m, h, h_data, tol);
}

/* stopping criteria: */

inline NLOPT_EXTERN(nlopt_result) nlopt_set_stopval(nlopt_opt opt, double stopval)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_stopval");
    return fun(opt, stopval);
}

inline NLOPT_EXTERN(double) nlopt_get_stopval(const nlopt_opt opt)
{
    static double(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (double(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_stopval");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_ftol_rel(nlopt_opt opt, double tol)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_ftol_rel");
    return fun(opt, tol);
}

inline NLOPT_EXTERN(double) nlopt_get_ftol_rel(const nlopt_opt opt)
{
    static double(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (double(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_ftol_rel");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_ftol_abs(nlopt_opt opt, double tol)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_ftol_abs");
    return fun(opt, tol);
}

inline NLOPT_EXTERN(double) nlopt_get_ftol_abs(const nlopt_opt opt)
{
    static double(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (double(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_ftol_abs");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_rel(nlopt_opt opt, double tol)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_xtol_rel");
    return fun(opt, tol);
}

inline NLOPT_EXTERN(double) nlopt_get_xtol_rel(const nlopt_opt opt)
{
    static double(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (double(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_xtol_rel");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_abs1(nlopt_opt opt, double tol)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_xtol_abs1");
    return fun(opt, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_abs(nlopt_opt opt, const double *tol)
{
    static nlopt_result(*fun)(nlopt_opt, const double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, const double *)) R_GetCCallable("nloptr","nlopt_set_xtol_abs");
    return fun(opt, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_get_xtol_abs(const nlopt_opt opt, double *tol)
{
    static nlopt_result(*fun)(nlopt_opt, double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double *)) R_GetCCallable("nloptr","nlopt_get_xtol_abs");
    return fun(opt, tol);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_maxeval(nlopt_opt opt, int maxeval)
{
    static nlopt_result(*fun)(nlopt_opt, int) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, int)) R_GetCCallable("nloptr","nlopt_set_maxeval");
    return fun(opt, maxeval);
}

inline NLOPT_EXTERN(int) nlopt_get_maxeval(const nlopt_opt opt)
{
    static int(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (int(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_maxeval");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_maxtime(nlopt_opt opt, double maxtime)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_maxtime");
    return fun(opt, maxtime);
}

inline NLOPT_EXTERN(double) nlopt_get_maxtime(const nlopt_opt opt)
{
    static double(*fun)(nlopt_opt) = NULL;
    if (fun == NULL) fun = (double(*)(nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_maxtime");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_force_stop(nlopt_opt opt)
{
    static nlopt_result(*fun)(nlopt_opt) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt)) R_GetCCallable("nloptr","nlopt_force_stop");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_force_stop(nlopt_opt opt, int val)
{
    static nlopt_result(*fun)(nlopt_opt, int) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, int)) R_GetCCallable("nloptr","nlopt_set_force_stop");
    return fun(opt, val);
}

inline NLOPT_EXTERN(int) nlopt_get_force_stop(const nlopt_opt opt)
{
    static int(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (int(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_force_stop");
    return fun(opt);
}

/* more algorithm-specific parameters */

inline NLOPT_EXTERN(nlopt_result) nlopt_set_local_optimizer(nlopt_opt opt, const nlopt_opt local_opt)
{
    static nlopt_result(*fun)(nlopt_opt, const nlopt_opt) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, const nlopt_opt)) R_GetCCallable("nloptr","nlopt_set_local_optimizer");
    return fun(opt, local_opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_population(nlopt_opt opt, unsigned pop)
{
    static nlopt_result(*fun)(nlopt_opt, unsigned) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, unsigned)) R_GetCCallable("nloptr","nlopt_set_population");
    return fun(opt, pop);
}

inline NLOPT_EXTERN(unsigned) nlopt_get_population(const nlopt_opt opt)
{
    static unsigned(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (unsigned(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_population");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_vector_storage(nlopt_opt opt, unsigned dim)
{
    static nlopt_result(*fun)(nlopt_opt, unsigned) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, unsigned)) R_GetCCallable("nloptr","nlopt_set_vector_storage");
    return fun(opt, dim);
}

inline NLOPT_EXTERN(unsigned) nlopt_get_vector_storage(const nlopt_opt opt)
{
    static unsigned(*fun)(const nlopt_opt) = NULL;
    if (fun == NULL) fun = (unsigned(*)(const nlopt_opt)) R_GetCCallable("nloptr","nlopt_get_vector_storage");
    return fun(opt);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_default_initial_step(nlopt_opt opt, const double *x)
{
    static nlopt_result(*fun)(nlopt_opt, const double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, const double *)) R_GetCCallable("nloptr","nlopt_set_default_initial_step");
    return fun(opt, x);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_initial_step(nlopt_opt opt, const double *dx)
{
    static nlopt_result(*fun)(nlopt_opt, const double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, const double *)) R_GetCCallable("nloptr","nlopt_set_initial_step");
    return fun(opt, dx);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_set_initial_step1(nlopt_opt opt, double dx)
{
    static nlopt_result(*fun)(nlopt_opt, double) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(nlopt_opt, double)) R_GetCCallable("nloptr","nlopt_set_initial_step1");
    return fun(opt, dx);
}

inline NLOPT_EXTERN(nlopt_result) nlopt_get_initial_step(const nlopt_opt opt, const double *x, double *dx)
{
    static nlopt_result(*fun)(const nlopt_opt, const double *, double *) = NULL;
    if (fun == NULL) fun = (nlopt_result(*)(const nlopt_opt, const double *, double *)) R_GetCCallable("nloptr","nlopt_get_initial_step");
    return fun(opt, x, dx);
}

#endif /* __NLOPTRAPI_H__ */

max_stars_repo_name: robalexclark/SilveR-Dev

max_issues_repo_name: robalexclark/SilveR-Dev

max_forks_repo_name: robalexclark/SilveR-Dev


text:
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 Microsoft Corporation. All rights reserved.
//
// This code is licensed under the MIT License (MIT).
//
// 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.
//
///////////////////////////////////////////////////////////////////////////////

#ifndef BITVEC_DECL_H
#define BITVEC_DECL_H

#include <gsl/gsl_assert> // for Expects
#include <gsl/gsl_util>   // for narrow_cast, narrow

using gsl::narrow_cast;

#include <cstddef>  // for ptrdiff_t, size_t, nullptr_t
#include <iterator> // for reverse_iterator, distance, random_access_...
#include <memory>

#ifdef _MSC_VER
#pragma warning(push)

// turn off some warnings that are noisy about our Expects statements
#pragma warning(disable : 4127) // conditional expression is constant
#pragma warning(disable : 4702) // unreachable code

// Turn MSVC /analyze rules that generate too much noise. TODO: fix in the tool.
#pragma warning(                                                               \
    disable : 26495) // uninitalized member when constructor calls constructor
#pragma warning(                                                               \
    disable : 26446) // parser bug does not allow attributes on some templates

#endif // _MSC_VER

namespace gb {

// shared_ptr without atomic overhead
template <typename T>
using shared_ptr_unsynchronized = std::__shared_ptr<T, __gnu_cxx::_S_single>;

// [views.constants], constants
const std::ptrdiff_t dynamic_extent = -1;

template <std::ptrdiff_t Extent = dynamic_extent> class bitvec;

// implementation details
namespace details {
template <class Span, bool IsConst> class span_iterator {
	using element_type_ = typename Span::element_type;

  public:
#ifdef _MSC_VER
	// Tell Microsoft standard library that span_iterators are checked.
	using _Unchecked_type = typename Span::pointer;
#endif

	using iterator_category = std::random_access_iterator_tag;
	using value_type = std::remove_cv_t<element_type_>;
	using difference_type = typename Span::index_type;

	using pointer = shared_ptr_unsynchronized<
	    std::conditional_t<IsConst, const element_type_, element_type_>>;

	span_iterator() = default;

	span_iterator(const Span *bitvec, typename Span::index_type idx) noexcept
	    : span_(bitvec), index_(idx) {}

	friend span_iterator<Span, true>;
	template <bool B, std::enable_if_t<!B && IsConst> * = nullptr>
	span_iterator(const span_iterator<Span, B> &other) noexcept
	    : span_iterator(other.span_, other.index_) {}

	std::conditional_t<IsConst, const bitvec<1>, bitvec<1>> operator*() const {
		Expects(index_ != span_->size());
		return std::conditional_t<IsConst, const bitvec<1>, bitvec<1>>(
		    (*span_)[index_]);
	}

	pointer operator->() const {
		Expects(index_ != span_->size());
		return span_->data() + index_;
	}

	span_iterator &operator++() {
		Expects(0 <= index_ && index_ != span_->size());
		++index_;
		return *this;
	}

	span_iterator &operator--() {
		Expects(index_ != 0 && index_ <= span_->size());
		--index_;
		return *this;
	}

	pointer operator[](difference_type n) const { return *(*this + n); }

	friend bool operator==(span_iterator lhs, span_iterator rhs) noexcept {
		return lhs.span_ == rhs.span_ && lhs.index_ == rhs.index_;
	}

	friend bool operator!=(span_iterator lhs, span_iterator rhs) noexcept {
		return !(lhs == rhs);
	}

	friend bool operator<(span_iterator lhs, span_iterator rhs) noexcept {
		return lhs.index_ < rhs.index_;
	}

#ifdef _MSC_VER
	// MSVC++ iterator debugging support; allows STL algorithms in 15.8+
	// to unwrap span_iterator to a pointer type after a range check in STL
	// algorithm calls
	friend void _Verify_range(
	    span_iterator lhs,
	    span_iterator rhs) noexcept {  // test that [lhs, rhs) forms a valid
		                               // range inside an STL algorithm
		Expects(lhs.span_ == rhs.span_ // range spans have to match
		        && lhs.index_ <= rhs.index_); // range must not be transposed
	}

	void _Verify_offset(const difference_type n) const
	    noexcept { // test that the iterator *this + n is a valid range in an
		           // STL
		// algorithm call
		Expects((index_ + n) >= 0 && (index_ + n) <= span_->size());
	}

	GSL_SUPPRESS(bounds .1)               // NO-FORMAT: attribute
	pointer _Unwrapped() const noexcept { // after seeking *this to a high water
		                                  // mark, or using one of the
		// _Verify_xxx functions above, unwrap this span_iterator to a raw
		// pointer
		return span_->data() + index_;
	}

	// Tell the STL that span_iterator should not be unwrapped if it can't
	// validate in advance, even in release / optimized builds:
	static const bool _Unwrap_when_unverified = false;
	GSL_SUPPRESS(con .3) // NO-FORMAT: attribute // TODO: false positive
	void _Seek_to(const pointer p) noexcept { // adjust the position of *this to
		                                      // previously verified location p
		// after _Unwrapped
		index_ = p - span_->data();
	}
#endif

  protected:
	const Span *span_ = nullptr;
	std::ptrdiff_t index_ = 0;
};

// Used for empty base class optimization
template <std::ptrdiff_t Ext> class extent_type {
  public:
	using index_type = std::ptrdiff_t;

	static_assert(Ext >= 0, "A fixed-size bitvec must be >= 0 in size.");

	extent_type() noexcept {}

	template <index_type Other> extent_type(extent_type<Other> ext) {
		static_assert(Other == Ext || Other == dynamic_extent,
		              "Mismatch between fixed-size extent and size of "
		              "initializing data.");
		Expects(ext.size() == Ext);
	}

	extent_type(index_type size) { Expects(size == Ext); }

	index_type size() const noexcept { return Ext; }
};

template <> class extent_type<dynamic_extent> {
  public:
	using index_type = std::ptrdiff_t;

	template <index_type Other>
	explicit extent_type(extent_type<Other> ext) : size_(ext.size()) {}

	explicit extent_type(index_type size) : size_(size) { Expects(size >= 0); }

	index_type size() const noexcept { return size_; }

  private:
	index_type size_;
};

template <std::ptrdiff_t Extent, std::ptrdiff_t Offset, std::ptrdiff_t Count>
struct calculate_subspan_type {
	using type =
	    bitvec<Count != dynamic_extent
	             ? Count
	             : (Extent != dynamic_extent ? Extent - Offset : Extent)>;
};
} // namespace details

} // namespace gb

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

#endif // BITVEC_DECL_H

max_stars_repo_name: CapacitorSet/FHE-tools

max_issues_repo_name: CapacitorSet/Glovebox

max_forks_repo_name: CapacitorSet/Glovebox


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 *  \file dnest_con.c
 *  \brief run dnest for continuum analysis
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>
#include <mpi.h>

#include "brains.h"

DNestFptrSet *fptrset_con;

/*!
 *  This function run dnest sampling for continuum.
 */
int dnest_con(int argc, char **argv)
{
  int i;

  num_params = parset.n_con_recon + num_params_var;

  par_range_model = malloc( num_params * sizeof(double *));
  par_prior_gaussian = malloc(num_params * sizeof(double *));
  for(i=0; i<num_params; i++)
  {
    par_range_model[i] = malloc(2*sizeof(double));
    par_prior_gaussian[i] = malloc(2*sizeof(double));
  }
  par_prior_model = malloc( num_params * sizeof(int));

  par_fix = (int *) malloc(num_params * sizeof(int));
  par_fix_val = (double *) malloc(num_params * sizeof(double));

  fptrset_con = dnest_malloc_fptrset();

  /* setup functions used for dnest*/
  fptrset_con->from_prior = from_prior_con;
  fptrset_con->perturb = perturb_con;
  fptrset_con->print_particle = print_particle_con;
  fptrset_con->restart_action = restart_action_con;
  fptrset_con->accept_action = accept_action_con;
  fptrset_con->kill_action = kill_action_con;
  fptrset_con->read_particle = read_particle_con;

  if(parset.flag_exam_prior != 1)
  {
    fptrset_con->log_likelihoods_cal = log_likelihoods_cal_con;
    fptrset_con->log_likelihoods_cal_initial = log_likelihoods_cal_initial_con;
    fptrset_con->log_likelihoods_cal_restart = log_likelihoods_cal_restart_con;
  }
  else
  {
    fptrset_con->log_likelihoods_cal = log_likelihoods_cal_con_exam;
    fptrset_con->log_likelihoods_cal_initial = log_likelihoods_cal_con_exam;
    fptrset_con->log_likelihoods_cal_restart = log_likelihoods_cal_con_exam;
  }

  set_par_range_con();

  /* setup fixed parameters */
  for(i=0; i<num_params; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }

  /* fix systematic error of continuum */
  if(parset.flag_con_sys_err != 1)
  {
    par_fix[0] = 1;
    par_fix_val[0] = log(1.0);
  }

  print_par_names_con();

  /* if not only print parameter name */
  if(parset.flag_para_name != 1)
    logz_con = dnest(argc, argv, fptrset_con, num_params, dnest_options_file);

  dnest_free_fptrset(fptrset_con);
  return 0;
}

/*!
 * this function set the parameter range.
 */
void set_par_range_con()
{
  int i;

  /* variability parameters */
  for(i=0; i<num_params_drw; i++)
  {
    par_range_model[i][0] = var_range_model[i][0];
    par_range_model[i][1] = var_range_model[i][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* parameters for long-term trend */
  for(i=num_params_drw; i<num_params_drw + num_params_trend; i++)
  {
    par_range_model[i][0] = var_range_model[3][0];
    par_range_model[i][1] = var_range_model[3][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }

  /* parameter for trend difference */
  for(i= num_params_drw + num_params_trend; i<num_params_var; i++)
  {
    par_range_model[i][0] = var_range_model[4 + i - (num_params_drw + num_params_trend)][0];
    par_range_model[i][1] = var_range_model[4 + i - (num_params_drw + num_params_trend)][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* continuum light curve parameters */
  for(i=num_params_var; i<num_params; i++)
  {
    par_range_model[i][0] = var_range_model[4+num_params_difftrend][0];
    par_range_model[i][1] = var_range_model[4+num_params_difftrend][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }
  return;
}

/*!
 *  print names and prior ranges for parameters
 *
 */
void print_par_names_con()
{
  if(thistask!= roottask)
    return;

  int i, j;
  FILE *fp;
  char fname[BRAINS_MAX_STR_LENGTH], str_fmt[BRAINS_MAX_STR_LENGTH];

  sprintf(fname, "%s/%s", parset.file_dir, "data/para_names_con.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
    exit(0);
  }

  strcpy(str_fmt, "%4d %-15s %10.6f %10.6f %4d %4d %15.6e\n");

  printf("# Print parameter name in %s\n", fname);

  fprintf(fp, "#*************************************************\n");
  fprint_version(fp);
  fprintf(fp, "#*************************************************\n");

  fprintf(fp, "%4s %-15s %10s %10s %4s %4s %15s\n", "#", "Par", "Min", "Max", "Prior", "Fix", "Val");
  i=0;
  fprintf(fp, str_fmt, i, "sys_err_con", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "sigmad", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "taud", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_trend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<num_params_difftrend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "diff trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<parset.n_con_recon; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "time series", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  fclose(fp);
}

/*!
 * This function generate a sample from the prior.
 */
void from_prior_con(void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    if(par_prior_model[i] == GAUSSIAN )
    {
      pm[i] = dnest_randn() * par_prior_gaussian[i][1] + par_prior_gaussian[i][0];
      dnest_wrap(&pm[i], par_range_model[i][0], par_range_model[i][1]);
    }
    else
    {
      pm[i] = var_range_model[i][0] + dnest_rand()*(par_range_model[i][1] - par_range_model[i][0]);
    }
  }

  for(i=0; i<num_params_var; i++)
  {
    if(par_fix[i] == 1)
      pm[i] = par_fix_val[i];
  }
  /* all parameters need to update at the initial step */
  which_parameter_update = -1;
}

/*!
 * This function calculate log likelihood probability.
 */
double log_likelihoods_cal_con(const void *model)
{
  double logL;
  logL = prob_con_variability_semiseparable(model);
  return logL;
}

/*!
 * This function calculate log likelihood probability at the initial step.
 */
double log_likelihoods_cal_initial_con(const void *model)
{
  double logL;
  logL = prob_con_variability_initial_semiseparable(model);
  return logL;
}

/*!
 * This function calculate log likelihood probability at the restart step.
 */
double log_likelihoods_cal_restart_con(const void *model)
{
  double logL;
  logL = prob_con_variability_initial_semiseparable(model);
  return logL;
}
/*!
 * This function generate a new move of parameters.
 */
double perturb_con(void *model)
{
  double *pm = (double *)model;
  double logH = 0.0, limit1, limit2, width, rnd;
  int which, which_level;

  /* sample variability parameters more frequently */
  do
  {
    rnd = dnest_rand();
    if(rnd < 0.1)
      which = dnest_rand_int(num_params_var);
    else
      which = dnest_rand_int(parset.n_con_recon) + num_params_var;
  }while(par_fix[which] == 1);

  which_parameter_update = which;

  /* level-dependent width */
  which_level_update = dnest_get_which_level_update();
  which_level = which_level_update > (size_levels - 10)?(size_levels-10):which_level_update;
  if( which_level > 0)
  {
    limit1 = limits[(which_level-1) * num_params *2 + which *2];
    limit2 = limits[(which_level-1) * num_params *2 + which *2 + 1];
    width = limit2 - limit1;
  }
  else
  {
    width = ( par_range_model[which][1] - par_range_model[which][0] );
  }

  if(par_prior_model[which] == GAUSSIAN)
  {
    logH -= (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
    pm[which] += dnest_randh() * width;
    dnest_wrap(&pm[which], par_range_model[which][0], par_range_model[which][1]);
    logH += (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
  }
  else
  {
    pm[which] += dnest_randh() * width;
    dnest_wrap(&(pm[which]), par_range_model[which][0], par_range_model[which][1]);
  }
  return logH;
}

/*!
 * This function print the particle into the file.
 */
void print_particle_con(FILE *fp, const void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    fprintf(fp, "%e ", pm[i] );
  }
  fprintf(fp, "\n");
}

/*!
 * This function read the particle from the file.
 */
void read_particle_con(FILE *fp, void *model)
{
  int j;
  double *psample = (double *)model;

  for(j=0; j < dnest_num_params; j++)
  {
    if(fscanf(fp, "%lf", psample+j) < 1)
    {
      printf("%f\n", *psample);
      fprintf(stderr, "#Error: Cannot read file %s.\n", options.sample_file);
      exit(0);
    }
  }
  return;
}

void restart_action_con(int iflag)
{

  return;
}

void accept_action_con()
{
  int param;

  param = which_parameter_update;
  /* only update prob when variability parameters are updated. */
  if(param < num_params_var)
  {
    prob_con_particles[which_particle_update] = prob_con_particles_perturb[which_particle_update];
  }
  return;
}

void kill_action_con(int i, int i_copy)
{
  prob_con_particles[i] = prob_con_particles[i_copy];
  return;
}

/*!
 * exam for prior
 */
double log_likelihoods_cal_con_exam(const void *model)
{
  return 0.0;
}

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
/*
 * Copyright (c) 2018-2021 Aleksas Mazeliauskas, Stefan Floerchinger,
 *                    Eduardo Grossi, and Derek Teaney
 * All rights reserved.
 *
 * FastReso is distributed under MIT license;
 * see the LICENSE file that should be present in the root
 * of the source distribution, or alternately available at:
 * https://github.com/amazeliauskas/FastReso/
 */
#ifndef FASTRESO_TFastReso_AZYHYDRO_h
#define FASTRESO_TFastReso_AZYHYDRO_h

#include "TFastReso.h"
#include <fstream>
#include <memory>
#include <map>
#include <vector>
#include <string>
#include <gsl/gsl_integration.h>

// The class storing the properties of particles and their irreducible distrubtion
// function components f_i.  Printing out routines are also defined there.
#include "TParticle_AZYHYDRO.h"

// List of parameters for GSL integration procedure.
#include "qag_params.h"

// List of parameters controlling the momentum discretization of irreducible components
// See also TParticle_AZYHYDRO class definition
#include "grid_params.h"

//! TFastReso_AZYHYDRO class is the top wrapper class of the fast freeze-out implementation.
class TFastReso_AZYHYDRO: public TFastReso {
  private:
    //! Data vector of TParticle_AZYHYDRO class of all read-in particles
    //! See TParticle_AZYHYDRO class definition.
    std::vector <std::unique_ptr<TParticle_AZYHYDRO>> fParticleData;

    //! Map between the particle name and its position in fParticleData vector;
    std::map <int, int> fIndexTable;

  public:
    //! reads the inputfile of resonances and creates TParticle_AZYHYDRO entry in fParticleData vector.
    //! This routine should be edited by the user to match their input file format.

    void read_particles_data(std::string inputname, std::string comps="Feq Fshear Fbulk Ftemp Fvel", bool verbose=false) override;

    //! The function call to initialized the particle distribution functions in fParticleData. Input parameters
    //! are the freeze-out temperature, chemical potential and (optionally), speed of sound (only needed for bulk perturbations).
    void do_thermal(double Tfo, double MuB=0, double MuI3=0, double MuS=0, double MuC=0, double Cs2=0.14) override;
    void do_thermal(int pi, double Tfo, double QMu=0, double Cs2 = 0.14 );

    //! Read in the input file of resonance decays and perfom the decays on the fly. It is important that decays
    //! are mass ordered to included all feeddown contributions to lower mass resonances.
    //! The table format is the same as decays.data in THERMINATOR 2 	[arXiv:1102.0273]
    //! Father Child1 Child2 (optional Child3) BranchingRatio ClebschGordanCoeff
    //! Note that ClebschGordanCoeff term is 0 or 1 identifying if the branching ration needs
    //! to be reweighted by actual Clebsch Gordon coefficients.
    void do_decays(std::string inputname, bool verbose=false) override;
  public:
    //! routine to get the pointer to TParticle_AZYHYDRO class for a particle particle name"
    TParticle_AZYHYDRO * getParticleByPDG(int pdgid){return fParticleData[fIndexTable[pdgid]].get();};
    TParticle_AZYHYDRO * getParticle(int i){return (i<0) ? fParticleData[fParticleData.size()+i].get()  : fParticleData[i].get();};
    int getParticleIndexByPDG(int pdgid){return fIndexTable[pdgid];};

    //! The printout of the calculated components is done by printing procedures in TParticle_AZYHYDRO class
};

#endif

max_stars_repo_name: amazeliauskas/FastReso

max_issues_repo_name: amazeliauskas/FastReso

max_forks_repo_name: amazeliauskas/FastReso


text:
#ifndef TETRA_DOS_DOS_VALUES_H
#define TETRA_DOS_DOS_VALUES_H

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include "ctetra/dos.h"
#include "elHamiltonian.h"
#include "environment.h"

gsl_matrix* cubicRecipLat(double a);

double* DosValues(Environment *env, int n, double *Es, double num_dos);

#endif // TETRA_DOS_DOS_VALUES_H


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

#pragma once

#include "Map.h"

#include "MapPoint.h"
#include "MappingKeyframe.h"
#include "InitializationData.h"

#include "Tracking\KeyframeBuilder.h"

#include "Mapping\MapPointKeyframeAssociations.h"

#include "BoW\OnlineBow.h"
#include "BoW\BaseFeatureMatcher.h"

#include "BundleAdjustment\BundleAdjust.h"

#include "Data\Types.h"
#include "Data\Pose.h"
#include "Data\Data.h"
#include "Data\MapState.h"

#include "Utils\historical_queue.h"
#include "Utils\collection.h"

#include <arcana\threading\blocking_concurrent_queue.h>

#include <gsl\gsl>

#include <shared_mutex>
#include <vector>
#include <set>

// forward-declare the friend class method for accessing the root map in unit tests
namespace UnitTests
{
    class TestHelperFunctions;
}

namespace mage
{
    class ThreadSafeMap
    {
        friend class ::UnitTests::TestHelperFunctions;

    public:
        ThreadSafeMap(
            const MageSlamSettings& settings, BaseBow& bagOfWords);

        ~ThreadSafeMap();

        void InitializeMap(const InitializationData& initializationData, thread_memory memory);

        void GetConnectedMapPoints(
            const collection<Proxy<MapPoint>>& points,
            thread_memory memory,
            loop::vector<KeyframeReprojection>& repro) const;

        // returns a sampling of mappoints
        void GetSampleOfMapPoints(temp::vector<MapPointProxy>& mapPoints, uint32_t max, uint32_t offset = 0) const;

        /*
        Inserts a new keyframe and returns all the keyframes that
        are connected to this one by the covisibility graph.
        */
        Proxy<Keyframe, proxy::Image> InsertKeyframe(
            const std::shared_ptr<const KeyframeBuilder>& kb,
            thread_memory memory);

        /*
        Inserts a new keyframe and returns all the keyframes that
        are connected to this one by the covisibility graph.
        */
        Proxy<Keyframe, proxy::Image> InsertKeyframe(
            const std::shared_ptr<const KeyframeBuilder>& keyframe,
            const gsl::span<MapPointAssociations<MapPointTrackingProxy>::Association>& extraAssociation,
            thread_memory memory);

        std::unique_ptr<KeyframeProxy> GetKeyFrameProxy(const Id<Keyframe>& kId) const;

        /*
        Merge information from a set of posited keyframes into the actual keyframes.
        */
        void UpdateKeyframesFromProxies(
            gsl::span<const KeyframeProxy> proxies,
            const OrbMatcherSettings& mergeSettings,
            std::unordered_map<Id<MapPoint>, MapPointTrackingProxy>& mapPointMerges,
            thread_memory memory);

        /*
        Updates the state of the map points in the KeyframeProxy based on the current
        information in the map.
        */
        void UpdateMapPoints(KeyframeProxy& proxy) const;

        /*
        Builds a global bundle adjust dataset of everything in the map.
        NOTE: Fixes the position of fixed and distance-tethered keyframes to preserve some semblance
        of consistent scale.
        */
        void BuildGlobalBundleAdjustData(AdjustableData& data) const;

        /*
        Retrieves keyframes similar to the image provided.
        */
        void FindSimilarKeyframes(const std::shared_ptr<const AnalyzedImage>& image, std::vector<KeyframeProxy>& output) const;

        void FindNonCovisibleSimilarKeyframeClusters(const KeyframeProxy& proxy, thread_memory memory, std::vector<std::vector<KeyframeProxy>>& output) const;

        /*
        Culls map points based on requirements related to number of keyframes they are visibile in and the results found in tracking
        Takes a list of points that have failed to appear where predicted during tracking to remove as well
        */
        void CullRecentMapPoints(
            const Id<Keyframe>& ki,
            const mira::unique_vector<Id<MapPoint>>& recentFailedMapPoints,
            thread_memory memory);

        /*
        Return the collection of all the map points
        */
        void GetMapPointsAsPositions(std::vector<Position>& positions) const;

        /*
        Return the number of map points in the map
        */
        size_t GetMapPointsCount() const;

        /*
        Return the number of keyframes in the map
        */
        size_t GetKeyframesCount() const;

        /*
        Return the collection of all the keyframes
        */
        void GetKeyframeViewMatrices(std::vector<Matrix>& viewMatrices) const;

        /*
        Creates the map points passed in and returns the keyframes that
        were added to the frames connected to Ki through the covis graph
        */
        std::vector<MappingKeyframe> CreateMapPoints(
            const Id<Keyframe>& keyframeId,
            const gsl::span<const MapPointKeyframeAssociations >& toCreate,
            thread_memory memory);

        /*
        Returns all the map points seen by this keyframe and the ones
        seen by it's neighbours in the covisibility graph. It then returns
        all the keyframes that see those points and aren't in Ki's covis graph.
        */
        unsigned int GetMapPointsAndDistantKeyframes(
            const Id<Keyframe>& Ki,
            unsigned int coVisTheta,
            thread_memory memory,
            std::vector<MapPointTrackingProxy>& mapPoints,
            std::vector<Proxy<Keyframe, proxy::Pose, proxy::Intrinsics, proxy::PoseConstraints>>& keyframes,
            std::vector<MapPointAssociation>& mapPointAssociations,
            std::vector<Id<Keyframe>>& externallyTetheredKeyframes) const;

        /*
        Updates all the map point positions and keyframe poses and
        destroys all the MapPoint outliers (Not sure if we actually destroy them or not).
        */
        void AdjustPosesAndMapPoints(
            const AdjustableData& adjusted,
            gsl::span<const std::pair<Id<MapPoint>, Id<Keyframe>>> outlierAssociations,
            thread_memory memory);

        /*
        Discards the duplicate keyframes that are
        connected to Ki in the covisibility graph.

        Returns the KeyframeIds that were removed as well as the covisability which lead to the removal decision if argument provided
        */
        void CullLocalKeyframes(
            Id<Keyframe> Ki,
            thread_memory memory,
            std::vector<std::pair<const Id<Keyframe>, const std::vector<Id<Keyframe>>>>* cullingDecisions = nullptr);

        /*
        gets keyframes connected to KId in the covisibility graph (more than the minimum required
        map points in common to be entered as a link in the covisibility graph)
        */
        void GetCovisibilityConnectedKeyframes(
            const Id<Keyframe>& kId,
            thread_memory memory,
            std::vector<MappingKeyframe>& kc) const;

        /*
        gets the list of recent points the mapping thread is still carefully evaluating. track local map
        will use this to collect more information on these points for use by the mapping thread, and the
        age of each of these map points
        */
        std::map<Id<MapPoint>, unsigned int> GetRecentlyCreatedMapPoints() const;

        /*
        * Will attempt to match the set of map points into the set of keyframes if matching criterea are met
        */
        void TryConnectMapPoints(gsl::span<const MapPointAssociations<MapPointTrackingProxy>::Association> mapPointAssociations,
            const Id<Keyframe>& insertedKeyframe,
            const TrackLocalMapSettings& trackLocalMapSettings,
            const OrbMatcherSettings& orbMatcherSettings,
            float searchRadius,
            thread_memory memory);

        /*
        Delete the state of the entire map
        */
        void Clear(thread_memory memory);

        /*
        Returns all the information pertinent when saving the map
        */
        MapState GetMapData() const;

        float GetMedianTetherDistance() const;

        float GetMapScale() const;

        /*
        Returns all the keyframes in the map
        */
        std::vector<KeyframeProxy> GetAllKeyframes() const;

        /*
        Destroys this instance of the thread safe map and returns it's
        internal map for use in another context.
        */
        static std::unique_ptr<Map> Release(std::unique_ptr<ThreadSafeMap> tsm);

    private:
        void UnSafeGetCovisibilityConnectedKeyframes(
            const Id<Keyframe>& kId,
            thread_memory memory,
            std::vector<MappingKeyframe>& kc) const;

        // the int value is the number of Keyframes this map point has existed for
        std::map<Id<MapPoint>, unsigned int> UnSafeGetRecentlyCreatedMapPoints() const;

        /*
        Performs all the logic for AdjustPosesAndMapPoints, but does not take its own
        lock and so can be used by other methods as well.
        */
        void UnSafeAdjustPosesAndMapPoints(
            const AdjustableData& adjusted,
            gsl::span<const std::pair<Id<MapPoint>, Id<Keyframe>>> outlierAssociations,
            thread_memory memory);

        // vector to record map point creations for map point culling
        historical_queue<std::vector<Proxy<MapPoint>>, 3> m_pointHistory;

        const KeyframeSettings& m_keyframeSettings;
        const MappingSettings& m_mappingSettings;
        const CovisibilitySettings& m_covisibilitySettings;
        const BundleAdjustSettings& m_bundleAdjustSettings;

        mutable std::shared_mutex m_mutex;

        std::unique_ptr<Map> m_map;

        BaseBow& m_bow;

        float m_mapScale;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#include <gsl/gsl_vector.h>
#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include "verif.h"
#include "data.h"
#include "prob_functions.h"

void EM (Dataset *data)
{
	int i, j;
	const double THRESHOLD = 1E-10;
	double Q, lastQ;

	srand(time(NULL));

	/* Initialize parameters to starting values */
	for (i = 0; i < data->numLabelers * data->numWorkflows; i++) {
		data->alpha[i] = data->priorAlpha[i];
		/*data->alpha[i] = (double) rand() / RAND_MAX * 4 - 2;*/
	}
	for (j = 0; j < data->numImages * data->numWorkflows; j++) {
		data->beta[j] = data->priorBeta[j];
		/*data->beta[j] = (double) rand() / RAND_MAX * 3;*/
	}

	Q = 0;
	EStep(data);
	Q = computeQ(data);
	/*printf("Q = %f\n", Q);*/
  int iter = 0;
	do {
  printf("Iteration %d ", iter++);
		lastQ = Q;
  printf("old Q: %f", lastQ);

		/* Re-estimate P(Z|L,alpha,beta) */
		EStep(data);
		Q = computeQ(data);
		/*printf("Q = %f; L = %f\n", Q, 0.0 computeLikelihood(data));*/
  printf("\tafter E: %f", Q);

		/*outputResults(data);*/
		MStep(data);

		Q = computeQ(data);
  printf("\tafter M: %f\n", Q);
		/*printf("Q = %f; L = %f\n", Q, 0.0 computeLikelihood(data));*/
	} while (fabs((Q - lastQ)/lastQ) > THRESHOLD);
	 outputResults(data);
}

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

    Dataset data;

	if (argc < 2) {
		fprintf(stdout, "Usage: em <data>\n");
		fprintf(stdout, "where the specified data file is formatted as described in the README file.\n");
		exit(1);
	}

	//These used to be 1
	//double prioralpha[2] = {1.0, 0.5};
	int NUMBER_OF_WORKER_POOLS;

	FILE *fi = fopen("../simulation.info", "r");
    fscanf(fi, "%d\n", &NUMBER_OF_WORKER_POOLS);

	double *prioralpha = (double *) malloc(sizeof(double) * NUMBER_OF_WORKER_POOLS);

    double temp1,temp2,temp3;
    fscanf(fi,"%lf\n",&temp1);
    fscanf(fi,"%lf,%lf\n",&temp2,&temp3);
    fscanf(fi,"%lf,%lf\n",&prioralpha[0],&prioralpha[1]);
    fclose(fi);

    prioralpha[0] = 0.5;
    prioralpha[1] = 0.5;

	double priorbeta = 0.5;
	//double prioralpha = 1.0;
	//double priorbeta = 1.0;

   //This code will now break, if extra arguments are passed, since prioralpha is
    // now an array.
  /*
  if ( argc > 2 )
    prioralpha = atof(argv[2]);
  if ( argc > 3 )
    priorbeta = atof(argv[3]);*/

  //printf("BLAH\n");
	readData(argv[1], &data, prioralpha, priorbeta);
	//	printf("LOADING TRUTH\n");
	//Load_truth(&data);
  //printf("STARTING EM\n");
	EM(&data);
	/*
  Accuracy(&data);
  int size;
  int s, samples = 50000;
  for ( size = 3; size <= 11; size+=2 ) {
    double acc = 0;
    for ( s = 0; s < samples; ++s ) {
      Randomizer( size);
      acc += Posterior( &data, size );
    }
    acc /= samples;
    printf("%d %f\n", size, acc);
  }

	free(data.priorAlpha);
	free(data.priorBeta);
	free(data.labels);
	free(data.alpha);
	free(data.beta);
	free(data.probZ1);
	free(data.probZ0);
	*/
}

max_stars_repo_name: ShreyaR/WorkerPoolSelection

max_issues_repo_name: ShreyaR/WorkerPoolSelection

max_forks_repo_name: ShreyaR/WorkerPoolSelection


text:
static char help[] = "ODE system solver example using TS.  Models two\n"
"equal-mass balls (particles) moving in the plane with either no\n"
"connection (free) or connected by a spring.  Solved via a Newtonian\n"
"formulation in cartesian coordinates (x,y) and velocities\n"
"(v=dx/dt,w=dy/dt). The system has dimension 8.\n\n";

// DEBUG check BDF2 convergence in free problem using Lagrangian formulation:
//for T in 2 3 4 5 7 8 9 10 11 12; do ./loose -ts_type bdf -ts_rtol 1.0e-$T -ts_atol 1.0e-$T; done

#include <petsc.h>

typedef struct {
    PetscBool    spring;
    PetscReal    g,     // m s-2;  acceleration of gravity
                 m,     // kg;     ball mass
                 l,     // m;      spring or rod length
                 k;     // N m-1;  spring constant
} LooseCtx;

extern PetscErrorCode SetInitial(Vec, LooseCtx*);
extern PetscErrorCode FreeExact(Vec, PetscReal, Vec, LooseCtx*);
extern PetscErrorCode NewtonRHSFcn(TS, PetscReal, Vec, Vec, void*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    PetscInt       steps;
    PetscReal      t0 = 0.0, tf = 2.0, dt = 0.1, errnorm;
    Vec            u, u0, uexact;
    TS             ts;
    LooseCtx       user;
    char           fstr[20] = "free";

    ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr;

    user.spring = PETSC_FALSE;
    user.g = 9.81;
    user.m = 58.0e-3; // 58 g for a tennis ball
    user.l = 0.5;     // spring length
    user.k = 20.0;    // spring constant (ignored by free)
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "lse_", "options for loose",
                             ""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-g", "acceleration of gravity (m s-2)",
                            "loose.c", user.g, &user.g, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-k", "spring constant (N m-1)", "loose.c",
                            user.k, &user.k, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-l", "spring length (m)", "loose.c",
                            user.l, &user.l, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-m", "mass of each ball (kg)", "loose.c",
                            user.m, &user.m, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-spring","spring case (not free)",
                            "loose.c", user.spring, &user.spring,
                            NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    ierr = VecCreate(PETSC_COMM_WORLD,&u); CHKERRQ(ierr);
    ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
    ierr = TSSetApplicationContext(ts,&user); CHKERRQ(ierr);
    ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);

    // set time axis
    ierr = TSSetTime(ts,t0); CHKERRQ(ierr);
    ierr = TSSetMaxTime(ts,tf); CHKERRQ(ierr);
    ierr = TSSetTimeStep(ts,dt); CHKERRQ(ierr);
    ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);

    // set up vecs and solver type
    ierr = VecSetSizes(u,PETSC_DECIDE,8); CHKERRQ(ierr);
    ierr = TSSetEquationType(ts, TS_EQ_ODE_EXPLICIT); CHKERRQ(ierr);
    ierr = TSSetRHSFunction(ts,NULL,NewtonRHSFcn,&user); CHKERRQ(ierr);
    ierr = TSSetType(ts,TSRK); CHKERRQ(ierr);

    // set-up of u, ts is complete
    ierr = VecSetFromOptions(u); CHKERRQ(ierr);
    ierr = TSSetFromOptions(ts); CHKERRQ(ierr);

    // set initial values and solve
    ierr = SetInitial(u, &user); CHKERRQ(ierr);
    ierr = TSSolve(ts, u); CHKERRQ(ierr);
    ierr = TSGetTime(ts, &tf); CHKERRQ(ierr);

    // numerical error in free case
    if (!user.spring) {
        ierr = VecDuplicate(u, &uexact); CHKERRQ(ierr);
        // get initial condition for evaluating exact solution
        ierr = VecDuplicate(u, &u0); CHKERRQ(ierr);
        ierr = SetInitial(u0, &user); CHKERRQ(ierr);
        ierr = FreeExact(u0, tf, uexact, &user); CHKERRQ(ierr);
        ierr = VecAXPY(u, -1.0, uexact); CHKERRQ(ierr); // u <- u - uexact
        ierr = VecNorm(u, NORM_INFINITY, &errnorm); CHKERRQ(ierr);
        VecDestroy(&u0);  VecDestroy(&uexact);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
               "numerical error at tf in free problem: |u-uexact|_inf = %.5e\n",
               errnorm); CHKERRQ(ierr);
    }

    // report
    ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr);
    if (user.spring)
        strcpy(fstr, "spring");
    ierr = PetscPrintf(PETSC_COMM_WORLD,
               "%s problem solved to tf = %.3f (%d steps)\n",
               fstr, tf, steps); CHKERRQ(ierr);

    VecDestroy(&u);  TSDestroy(&ts);
    return PetscFinalize();
}

/* Variable order used in the following functions:
u[0] = q_1 (= x_1)
u[1] = q_2 (= y_1)
u[2] = q_3 (= x_2)
u[3] = q_4 (= y_2)
u[4] = v_1 (= dx_1/dt)
u[5] = v_2 (= dy_1/dt)
u[6] = v_3 (= dx_2/dt)
u[7] = v_4 (= dy_2/dt)
*/

PetscErrorCode SetInitial(Vec u, LooseCtx *user) {
    /* Set initial conditions. */
    PetscErrorCode   ierr;
    PetscReal        *au;
    ierr = VecGetArray(u,&au); CHKERRQ(ierr);
    au[0] = 0.0;              // q1 = x1
    au[1] = 1.0;              // q2 = y1
    au[2] = 0.0;              // q3 = x2
    au[3] = au[1] + user->l;  // q4 = y2;  satisfies (1)
    au[4] = 10.0;             // v1
    au[5] = 10.0;             // v2
    au[6] = 15.0;             // v3
    au[7] = au[5];            // v4;  satisfies (2)
    ierr = VecRestoreArray(u,&au); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode FreeExact(Vec u0, PetscReal tf, Vec uexact, LooseCtx *user) {
    /* Exact solution based on parabolic motion.  Problem
        m x'' = 0,       x(0) = x0,  x'(0) = v0
        m y'' = - m g,   y(0) = y0,  y'(0) = w0
    has solution
        x(t) = x0 + v0 t,              x'(t) = v0,
        y(t) = y0 + w0 t - (g/2) t^2,  y'(t) = w0 - g t    */
    PetscErrorCode ierr;
    PetscReal *au0, *auex;
    if (user->spring) {
        SETERRQ(PETSC_COMM_SELF,7,"exact solution only implemented for free\n");
    }
    ierr = VecGetArray(u0, &au0); CHKERRQ(ierr);
    ierr = VecGetArray(uexact, &auex); CHKERRQ(ierr);
    auex[0] = au0[0] + au0[4] * tf;
    auex[1] = au0[1] + au0[5] * tf - 0.5 * user->g * tf * tf;
    auex[2] = au0[2] + au0[6] * tf;
    auex[3] = au0[3] + au0[7] * tf - 0.5 * user->g * tf * tf;
    auex[4] = au0[4];
    auex[5] = au0[5] - user->g * tf;
    auex[6] = au0[6];
    auex[7] = au0[7] - user->g * tf;
    ierr = VecRestoreArray(u0, &au0); CHKERRQ(ierr);
    ierr = VecRestoreArray(uexact, &auex); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode NewtonRHSFcn(TS ts, PetscReal t, Vec u, Vec G, void *ctx) {
    PetscErrorCode   ierr;
    LooseCtx         *user = (LooseCtx*)ctx;
    const PetscReal  *au;
    PetscReal        *aG;
    PetscReal        dspring = 0.0, cspring = 0.0, Fx1, Fy1, Fx2, Fy2;

    PetscFunctionBeginUser;
    ierr = VecGetArrayRead(u,&au); CHKERRQ(ierr);
    ierr = VecGetArray(G,&aG); CHKERRQ(ierr);
    Fx1 = 0.0;
    Fy1 = - user->m * user->g;
    Fx2 = 0.0;
    Fy2 = - user->m * user->g;
    if (user->spring) {
        dspring = PetscSqrtReal(  (au[0] - au[2]) * (au[0] - au[2])
                                + (au[1] - au[3]) * (au[1] - au[3]) );
        if (dspring == 0.0) {
            SETERRQ(PETSC_COMM_SELF,1,"exact ball collision (unlikely?)\n");
        }
        cspring = user->k * (1.0 - user->l / dspring);
        Fx1 -= cspring * (au[0] - au[2]);
        Fy1 -= cspring * (au[1] - au[3]);
        Fx2 += cspring * (au[0] - au[2]);
        Fy2 += cspring * (au[1] - au[3]);
    }
    aG[0] = au[4];
    aG[1] = au[5];
    aG[2] = au[6];
    aG[3] = au[7];
    aG[4] = Fx1 / user->m;
    aG[5] = Fy1 / user->m;
    aG[6] = Fx2 / user->m;
    aG[7] = Fy2 / user->m;
    ierr = VecRestoreArrayRead(u,&au); CHKERRQ(ierr);
    ierr = VecRestoreArray(G,&aG); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

max_stars_repo_name: bueler/dae-examples

max_issues_repo_name: bueler/dae-examples

max_forks_repo_name: bueler/dae-examples