text:
/* DomSelfAdaptFWD.c

Forward-in-time simulation of an adaptive alleles with dominance and selfing, with linked neutral fragment.

This is the BATCH version - to be used on cluster machine to produce many simulation replicates at once.

Reps program: Reads in pre-calculated tables and introduces neutral mutations based on them. See README for more information.

Simulation uses routines found with the GNU Scientific Library (GSL)
(http://www.gnu.org/software/gsl/)
Since GSL is distributed under the GNU General Public License
(http://www.gnu.org/copyleft/gpl.html), you must download it
separately from this file.

*/

/* Preprocessor statements */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stddef.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <sys/stat.h>
#include <sys/types.h>

/* Function prototypes */
void Wait();
void fitness(unsigned int N, double h, double s, unsigned int *inpop, double *outfit, double *outcf, double *fitsum);
double sumT_UI(unsigned int *Tin, unsigned int nrow);
unsigned int parchoose(unsigned int N, double *culfit, double *fitsum, const gsl_rng *r);
unsigned int bitswitch(unsigned int x);
void rec_sort(double *index, unsigned int nrow);
void generation(unsigned int N, double self, double R, unsigned int *nums, unsigned int *selin, unsigned int *selout, unsigned int **neutin, unsigned int **neutout, double *culfit, double *fitsum, unsigned int npoly, double *polypos, const gsl_rng *r);
void addpoly(unsigned int N, unsigned int **neutin, double *polyloc, unsigned int *npoly, double theta, unsigned int maxsize, const gsl_rng *r);
void reassign(unsigned int **neutin, unsigned int **neutout, unsigned int *selin, unsigned int *selout, unsigned int npoly, unsigned int N);
void reassign2(unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int N);
void polyprint(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N);
void ptrim(unsigned int size, unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int *npolyT, double *avpi);
void mutsamp(unsigned int N, unsigned int samps, unsigned int *selin, double *polypos, unsigned int **neutin, unsigned int *nums2, unsigned int npoly, unsigned int ei, unsigned int rep, const gsl_rng *r);
void popread(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N, unsigned int ei, unsigned int suffix);

void Wait(){
	printf("Press Enter to Continue");
	while( getchar() != '\n' );
	printf("\n");
}

/* Fitness calculation routine */
void fitness(unsigned int N, double h, double s, unsigned int *inpop, double *outfit, double *outcf, double *fitsum){
	unsigned int i;				/* Pop counter */
	unsigned int ac = 0;		/* Allele copies in indv */

	/* Calculating new fitnesses */
    *fitsum = 0;
	for(i = 0; i < N; i++){
		ac = 0;
		ac = (*(inpop + 2*i)) + (*(inpop + 2*i + 1));
		if(ac == 0){
			*(outfit + i) = 1;
		}else if(ac == 1){
			*(outfit + i) = 1 + h*s;
		}else if(ac == 2){
			*(outfit + i) = 1 + s;
		}
		*fitsum += *(outfit + i);
		*(outcf + i) = *fitsum;
    }
}	/* End of fitness routine */

/* Summing entire table (unsigned int) */
double sumT_UI(unsigned int *Tin, unsigned int nrow){
	unsigned int i;
	double res = 0;
	for(i = 0; i < nrow; i++){
		res += (*(Tin + i));
	}
	return res;
}

/* Parent choosing */
unsigned int parchoose(unsigned int N, double *culfit, double *fitsum, const gsl_rng *r){

	double a;					/* Random number */
	unsigned int done = 0;		/* Signal to determine whether selection complete */
	unsigned int cases = 0;		/* Determines whether one should add up or decrease */
	unsigned int choose1 = 0;	/* Chromosome to be selected */
	unsigned int firstgo = 0;	/* First try at choosing progeny */

	done = 0;
	cases = 0;
	choose1 = 0;
	a = gsl_ran_flat(r,0.0,1.0);
	firstgo = (int)N*a;

	/* Choosing case */
	if (((*(culfit + firstgo)/(*(fitsum))) >= a) && firstgo == 0){
		cases = 3;
	} else if (((*(culfit + firstgo)/(*(fitsum))) >= a) && firstgo != 0){
		cases = 2;
	} else {
		cases = 1;
	}

	switch(cases)
	{
		case 1:
			while(done != 1){
				firstgo++;
				if((*(culfit + firstgo)/(*(fitsum))) >= a){
					choose1 = firstgo;
					done = 1;
				}
			}
			break;
		case 2:
			while(done != 1){
				firstgo--;
				if((*(culfit + firstgo)/(*(fitsum))) < a || firstgo == 0){
					choose1 = (firstgo + 1);
					done = 1;
				}
			}
		break;
		case 3:
			choose1 = 0;
			break;
	}

	return choose1;
}

/* Bit-switching routine */
unsigned int bitswitch(unsigned int x){
	unsigned int ret;
	switch(x)
	{
		case 0:
			ret = 1;
			break;
		case 1:
			ret = 0;
			break;
		default:	/* If none of these cases chosen, exit with error message */
			fprintf(stderr,"Error in bitswitch: input is not 0 or 1.\n");
			exit(1);
            break;
	}
	return ret;
}

/* Sorting rec events after choosing them */
void rec_sort(double *index, unsigned int nrow){
	unsigned int j, i;		/* Sorting indices */
	double temp0;		/* For swapping */

	for(j = 1; j < nrow; j++){
		i = j;
		while( (i > 0) && ( *(index + (i - 1) ) > *(index + i) )){
			/* Swapping entries */
			temp0 = *(index + (i - 1));
			*(index + (i - 1)) =  *(index + (i));
			*(index + (i)) = temp0;
			i--;
		}
	}

}

/* Reproduction routine */
void generation(unsigned int N, double self, double R, unsigned int *nums, unsigned int *selin, unsigned int *selout, unsigned int **neutin, unsigned int **neutout, double *culfit, double *fitsum, unsigned int npoly, double *polypos, const gsl_rng *r){
	unsigned int i, j, k;		/* Pop counter, neutral marker counter, rec event counter */
	unsigned int isself = 0;	/* Is this a selfing reproduction? */
	unsigned int choose1 = 0;	/* Chromosome to be selected */
	unsigned int choose2 = 0;	/* Chromosome to be selected (2nd with outcrossing) */
	unsigned int wc1 = 0;		/* Which chrom (parent 1) */
	unsigned int wc2 = 0;		/* Which chrom (parent 2) */
	unsigned int index1 = 0;	/* Index of chosen sample 1 */
	unsigned int index2 = 0;	/* Index of chosen sample 2 */
	unsigned int nself = 0;		/* Number of selfing events */
	unsigned int selfix = 0;	/* Selfing index */
	unsigned int recix1 = 0;	/* Rec index chr 1 */
	unsigned int recix2 = 0;	/* Rec index chr 2 */
	unsigned int nrec1 = 0;		/* Number rec events chr 1 */
	unsigned int nrec2 = 0;		/* Number rec events chr 2 */

	/* First deciding number of selfing events; creating vector of events + choosing */
	nself = gsl_ran_binomial(r, self, N);
	unsigned int *selfev = calloc(nself,sizeof(unsigned int));
	gsl_ran_choose(r, selfev, nself, nums, N, sizeof(unsigned int));

	for(i = 0; i < N; i++){		/* Regenerating population */

    	/* Choosing first parent, and relevant chromosomes */
		choose1 = parchoose(N, culfit, fitsum, r);
		wc1 = gsl_ran_bernoulli(r,0.5);
		wc2 = gsl_ran_bernoulli(r,0.5);
		index1 = 2*choose1 + wc1;

		/* Copying over selected sites, first checking if selfing occurs */
		isself = 0;
		if(nself != 0){
			if(*(selfev + selfix) == i){
				isself = 1;
				selfix++;
			}
		}
		if(isself == 1){
			index2 = 2*choose1 + wc2;
		}else if(isself == 0){
			choose2 = parchoose(N, culfit, fitsum, r);
			index2 = 2*choose2 + wc2;
		}
		(*(selout + 2*i)) = (*(selin + index1));
		(*(selout + 2*i + 1)) = (*(selin + index2));

		/* Now copying over neutral fragment, accounting for recombination */
		if(R != 0){
			/* Choosing number of recombination events on arms 1, 2 */
			recix1 = 0;
			recix2 = 0;
			nrec1 = gsl_ran_poisson(r,R);
			nrec2 = gsl_ran_poisson(r,R);
			double *recev1 = calloc(nrec1 + 1,sizeof(double));
			double *recev2 = calloc(nrec2 + 1,sizeof(double));
			for(k = 0; k < nrec1; k++){
				*(recev1 + k) = gsl_ran_flat(r,0,1);
			}
			for(k = 0; k < nrec2; k++){
				*(recev2 + k) = gsl_ran_flat(r,0,1);
			}
			*(recev1 + nrec1) = 1.01;
			*(recev2 + nrec2) = 1.01;
			rec_sort(recev1,nrec1 + 1);
			rec_sort(recev2,nrec2 + 1);

			for(j = 0; j < npoly; j++){
				if( (*(polypos + j)) >= *(recev1 + recix1)){
					wc1 = bitswitch(wc1);
					index1 = 2*choose1 + wc1;
					recix1++;
				}
				if( (*(polypos + j)) >= *(recev2 + recix2)){
					wc2 = bitswitch(wc2);
					if(isself == 1){
						index2 = 2*choose1 + wc2;
					}else if(isself == 0){
						index2 = 2*choose2 + wc2;
					}
					recix2++;
				}

				*((*(neutout + 2*i)) + j) = *((*(neutin + index1)) + j);
				*((*(neutout + 2*i + 1)) + j) = *((*(neutin + index2)) + j);

			}

			free(recev1);
			free(recev2);

		}else if (R == 0){
			for(j = 0; j < npoly; j++){
				*((*(neutout + 2*i)) + j) = *((*(neutin + index1)) + j);
				*((*(neutout + 2*i + 1)) + j) = *((*(neutin + index2)) + j);
			}
		}
	}

	free(selfev);

}	/* End of reproduction routine */

/* Adding neutral polymorphism routine */
void addpoly(unsigned int N, unsigned int **neutin, double *polyloc, unsigned int *npoly, double theta, unsigned int maxsize, const gsl_rng *r){
	unsigned int i, j;
	int k;
	unsigned int newpoly = 0;			/* Number of new polymorphisms added */
	unsigned int pchr = 0;				/* Chromosome of appearance */
	unsigned int pfound = 0;			/* Flag to denote if right index found */
	double ploc = 0;					/* Location of polymorphism */
/*	FILE *ofp_poly;				 Pointer for polymorphisms output */

	newpoly = gsl_ran_poisson(r,(0.5*theta));
/*	printf("Add %d new polymorphisms\n",newpoly);*/

	for(j = 0; j < newpoly; j++){
		ploc = gsl_ran_flat(r,0.0,1.0);
/*		printf("ploc is %lf\n",ploc);		*/
		pchr = gsl_rng_uniform_int(r,2*N);
/*		printf("pchr is %d\n",pchr);		*/

		/* Inserting new polymorphism */
		if(*(npoly) == 0){
/*			printf("zero entry here\n");	*/
			*(polyloc + 0) = ploc;
				for(i = 0; i < 2*N; i++){
					*((*(neutin + i)) + 0) = 0;
				}
			*((*(neutin + pchr)) + 0) = 1;
		}else if(*(npoly) > 0){
/*			printf("non-zero entry here\n");		*/
			pfound = 0;
			for(k = (*(npoly) - 1); k >= 0; k--){
/*				printf("k is %d\n",k);*/
				if(*(polyloc + k) > ploc){
/*					printf("is here 1\n");*/
					*(polyloc + k + 1) = *(polyloc + k);
					for(i = 0; i < 2*N; i++){
						*((*(neutin + i)) + k + 1) = *((*(neutin + i)) + k);
					}
				}else if(*(polyloc + k) <= ploc){
/*					printf("is here 2\n");				*/
					pfound = 1;
					*(polyloc + k + 1) = ploc;
					for(i = 0; i < 2*N; i++){
						*((*(neutin + i)) + k + 1) = 0;
					}
					*((*(neutin + pchr)) + k + 1) = 1;
					break;
				}
			}

			if(pfound == 0){		/* In this case polymorphism is inserted at the start */
				*(polyloc + 0) = ploc;
				for(i = 0; i < 2*N; i++){
					*((*(neutin + i)) + 0) = 0;
				}
				*((*(neutin + pchr)) + 0) = 1;
			}

		}

		*(npoly) += 1;
/*		printf("n poly is %d\n",*(npoly));			*/
		if(*(npoly) == maxsize){
			fprintf(stderr,"Too many neutral polymorphisms.\n");
			exit(1);
		}

	}
}

/* Reassigning new population routine */
void reassign(unsigned int **neutin, unsigned int **neutout, unsigned int *selin, unsigned int *selout, unsigned int npoly, unsigned int N){
	unsigned int i, j;		/* pop counter, poly counter */

	for(i = 0; i < 2*N; i++){
		*(selout + i) = *(selin + i);
		for(j = 0; j < npoly; j++){
			*((*(neutout + i)) + j) = *((*(neutin + i)) + j);
		}
	}

}	/* End of reassignment routine */

/* Reassigning new population routine (after trimming) */
void reassign2(unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int N){
	unsigned int i, j;		/* pop counter, poly counter */

	for(i = 0; i < 2*N; i++){
		for(j = 0; j < npoly; j++){
			*((*(neutout + i)) + j) = *((*(neutin + i)) + j);
			if(i == 0){
				*(posout + j) = *(posin + j);
			}
		}
	}

}	/* End of reassignment routine */

/* Printing polymorphism routines */
void polyprint(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N){
	unsigned int i, j;		/* pop counter, poly counter */

	for(j = 0; j < npoly; j++){
		printf("%lf ",*(posin + j));
		for(i = 0; i < 2*N; i++){
			printf("%d ",*((*(neutin + i)) + j));
		}
		printf("\n");
	}
	Wait();

}	/* End of reassignment routine */

/* Cutting out non-polymorphic sites */
void ptrim(unsigned int size, unsigned int **neutin, unsigned int **neutout, double *posin, double *posout, unsigned int npoly, unsigned int *npolyT, double *avpi){
	unsigned int i, j;
	unsigned int newp = 0;		/* New number of polymorphisms */
	unsigned int count = 0;		/* Polymorphism count */
	double cumpi = 0;			/* Cumulative pi (diversity) */

/*	printf("Trimming routine activated\n");	*/
	for(j = 0; j < npoly; j++){
		count = 0;
		*(posout + newp) = *(posin + j);
		for(i = 0; i < size; i++){
			*((*(neutout + i)) + newp) = *((*(neutin + i)) + j);
			count += *((*(neutin + i)) + j);
		}
/*		printf("For poly %d, count is %d\n",j,count);*/
		if((count > 0) && (count < size)){
/*			printf("Keeping poly %d\n",j);	*/
			newp++;
			cumpi += ((count*(size-count))/(1.0*size*size));
/*			printf("Count is %d, freq is %lf, cumulative value is %lf\n",count,count/(1.0*size),cumpi);	*/
		}
	}

	*npolyT = newp;
	*avpi = (cumpi/(1.0*newp));
}

/* Producing mutational samples */
void mutsamp(unsigned int N, unsigned int samps, unsigned int *selin, double *polypos, unsigned int **neutin, unsigned int *nums2, unsigned int npoly, unsigned int ei, unsigned int rep, const gsl_rng *r){

	unsigned int a, j, x;
	unsigned int currsamp;
	unsigned int npolyT;
	double avpi;
	char Mout[32];				 /* String to hold filename in (Mutations) */
	FILE *ofp_mut = NULL;		 /* Pointer for data output */

	for(x = 0; x < rep; x++){

		double *polyposF = calloc(npoly,sizeof(double));						/* Position of neutral mutations (final for polymorphism sampling) */
		unsigned int *selsamp = calloc(samps,sizeof(unsigned int *));			/* State of selected locus after sim stops */
		unsigned int **nsamp = calloc(samps,sizeof(unsigned int *));			/* Table of neutral markers per individual (final for sampling) */
		unsigned int **nsampT = calloc(samps,sizeof(unsigned int *));			/* Table of neutral markers per individual (final for sampling) AFTER TRIMMING */
		for(a = 0; a < samps; a++){
			nsamp[a] = calloc(npoly,sizeof(unsigned int));
			nsampT[a] = calloc(npoly,sizeof(unsigned int));
		}

		npolyT = npoly;

		/* Choosing which haplotypes to sample */
		unsigned int *thesamp = calloc(samps,sizeof(unsigned int));
		gsl_ran_choose(r, thesamp, samps, nums2, 2*N, sizeof(unsigned int));

		for(a = 0; a < samps; a++){
			currsamp = *(thesamp + a);
			*(selsamp + a) = *(selin + currsamp);
			for(j = 0; j < npoly; j++){
				*((*(nsamp + a)) + j) = *((*(neutin + currsamp)) + j);
			}
		}

		/* Check, then include code to ONLY print out polymorphic sites in table? */
		ptrim(samps, nsamp, nsampT, polypos, polyposF, npoly, &npolyT, &avpi);

		/* Printing out sample table */
		sprintf(Mout,"Mutations/Muts_%d.dat",ei*rep + x);
		ofp_mut = fopen(Mout,"w");

		/* State of selected site */
		fprintf(ofp_mut,"0 ");
		for(a = 0; a < samps; a++){
			fprintf(ofp_mut,"%d ",*(selsamp + a));
		}
		fprintf(ofp_mut,"\n");

		for(j = 0; j < npolyT; j++){
			fprintf(ofp_mut,"%lf ",*(polyposF + j));
			for(a = 0; a < samps; a++){
				fprintf(ofp_mut,"%d ",*((*(nsampT + a)) + j));
			}
			fprintf(ofp_mut,"\n");
		}
		fclose(ofp_mut);

		for(a = 0; a < samps; a++){
			free(nsampT[a]);
			free(nsamp[a]);
		}
		free(nsampT);
		free(nsamp);
		free(thesamp);
		free(selsamp);
		free(polyposF);

	}
}

/* Reading population state from file */
void popread(unsigned int **neutin, double *posin, unsigned int npoly, unsigned int N, unsigned int ei, unsigned int suffix){
	unsigned int i, j;		/* pop counter, poly counter */
	char Pin[32];				 /* String to hold filename in (Mutations) */
	FILE *ifp_poly = NULL;		 /* Pointer for data output */

	/* Printing out polymorphism table */
	sprintf(Pin,"Pop_%d.dat",ei + suffix);
	ifp_poly = fopen(Pin,"r");
	for(j = 0; j < npoly; j++){
		/* Is this the right way of reading in the first column of file? */
		for(i = 0; i < (2*N + 1); i++){
			if(i == 0){
				fscanf(ifp_poly,"%lf",&(posin[j]));
			}else if(i > 0){
				fscanf(ifp_poly,"%d",&(neutin[(i-1)][j]));
			}
		}
		fprintf(ifp_poly,"\n");
	}
	fclose(ifp_poly);

}	/* End of population reading routine */

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

	/* Declare variables here */
	unsigned int a, i, x;		/* Counters */
	unsigned int N = 0;			/* Population Size */
	unsigned int initp = 0;		/* Initial location of mutant (individual) */
	unsigned int done = 0;		/* Is simulation done or not? */
	unsigned int dcopies = 0;	/* Number of copies of derived allele */
	unsigned int afix = 0;		/* Has allele fixed? */
	unsigned int npoly = 0;		/* Number of neutral alleles */
	unsigned int npolyB = 0;	/* BASELINE number of neutral alleles at simulation initiation */
	unsigned int npolyT = 0;	/* Number of neutral alleles (after cutting out fixed sites) */
	unsigned int nowsel = 0;	/* If derived allele selected for? */
	unsigned int samps = 0;		/* Number of samples to take */
	unsigned int suffix = 0;	/* Number of times to subsample from population */
	unsigned int base = 0;		/* Number of baseline forward-in-time simulations */
	unsigned int rps = 0;		/* Samples taken per baseline simulation */
	unsigned int INITTS = 0;	/* Initial size of mutation array */
/*	unsigned int n = 0;			sprintf counter */
	double s = 0;				/* Strength of selection */
	double h = 0;				/* Dominance level */
	double self = 0;			/* Selfing rate */
	double Rin = 0;				/* Recombination rate (input) */
	double R = 0;				/* Recombination rate (at a time) */
	double x0 = 0;				/* Initial allele frequency */
	double xfix = 0;			/* Final allele frequency (i.e. stop when allele reaches this frequency) */
	double theta = 0;			/* Rate of neutral mutation */
	double fitsum = 0;			/* Summed Fitness */
	double scurr = 0;			/* Current fitness of mutant (to account for initial neutrality) */
	double avpi = 0;
	char Pout[32];				/* String to hold filename in (Mutations) */
	char Sout[32];				/* String to hold filename in (Seed) */
	FILE *ofp_sd;				/* Pointer for seed output */
	FILE *ofp_poly;				/* Pointer for polymorphisms output */

	/* GSL random number definitions */
	const gsl_rng_type * T;
	gsl_rng * r;

	/* Reading in data from command line */
	if(argc != 14){
		fprintf(stderr,"Not enough inputs (need: N s h self R x0 xfix theta basereps reps-per-sim samples suffix npoly).\n");
		exit(1);
	}

	N = atoi(argv[1]);
	if(N <= 0){
		fprintf(stderr,"Total Population size N is zero or negative, not allowed.\n");
		exit(1);
	}

	s = strtod(argv[2],NULL);
	if(s < 0){
		fprintf(stderr,"Allele strength s is negative, not allowed.\n");
		exit(1);
	}

	h = strtod(argv[3],NULL);
	if(h < 0 || h > 1){
		fprintf(stderr,"Dominance value must lie between 0 and 1.\n");
		exit(1);
	}

	self = strtod(argv[4],NULL);
	if(self < 0 || self > 1){
		fprintf(stderr,"Selfing rate must lie between 0 and 1 (inclusive).\n");
		exit(1);
	}

	Rin = strtod(argv[5],NULL);
	if(Rin < 0){
		fprintf(stderr,"Recombination rate must be positive or zero.\n");
		exit(1);
	}
	Rin /= 2.0*N;

	x0 = strtod(argv[6],NULL);
	if(x0 < (1.0/(2.0*N)) || x0 > 1.0-(1.0/(2.0*N)) ){
		fprintf(stderr,"Initial mutant frequency must lie between 1/2N and 1-1/2N.\n");
		exit(1);
	}

	xfix = strtod(argv[7],NULL);
	if(xfix < 0 || xfix > 1.0 ){
		fprintf(stderr,"Final mutant frequency must lie between 0 and 1.\n");
		exit(1);
	}
	if(xfix <= x0){
		fprintf(stderr,"Final mutant frequency must be greater than initial mutant frequency.\n");
		exit(1);
	}

	theta = strtod(argv[8],NULL);
	if(theta <= 0){
		fprintf(stderr,"Mutation rate must be a positive value.\n");
		exit(1);
	}

	base = atoi(argv[9]);
	if(base <= 0){
		fprintf(stderr,"Number of baseline simulations must be greater than zero.\n");
		exit(1);
	}

	rps = atoi(argv[10]);
	if(rps <= 0){
		fprintf(stderr,"Number of reps per simulation must be > 0.\n");
		exit(1);
	}

	samps = atoi(argv[11]);
	if(samps <= 0){
		fprintf(stderr,"Number of samples taken must be > 0.\n");
		exit(1);
	}

	suffix = atoi(argv[12]);
	if(argv[12] < 0){
		fprintf(stderr,"File index must be greater than or equal to zero.\n");
		exit(1);
	}

	npolyB = atoi(argv[13]);
	INITTS = 2*npolyB;

	/* create a generator chosen by the
    environment variable GSL_RNG_TYPE */

    mkdir("SeedsRep/", 0777);
	gsl_rng_env_setup();
	if (!getenv("GSL_RNG_SEED")) gsl_rng_default_seed = time(0);
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	sprintf(Sout,"SeedsRep/Seed_%d.dat",suffix);
	ofp_sd = fopen(Sout,"w");
	fprintf(ofp_sd,"%lu\n",gsl_rng_default_seed);
	fclose(ofp_sd);

	mkdir("Polymorphisms/", 0777);
	mkdir("Mutations/", 0777);

	/* Vector of 0 to (N-1) for sampling random numbers */
	unsigned int *nums = calloc(N,sizeof(unsigned int));
	unsigned int *nums2 = calloc(2*N,sizeof(unsigned int));
	for(a = 0; a < 2*N; a++){
		*(nums2 + a) = a;
		if(a < N){
			*(nums + a) = a;
		}
	}

	/* Executing simulation */
	for(i = 0; i < base; i++){

/*		printf("Starting run %d\n",i);*/

		afix = 0;
		R = Rin;
		while(afix != 1){

			/* Setting up selection, neutral tables per individual */
			double *fit = calloc(N,sizeof(double));									/* Fitness of each individual */
			double *cumfit = calloc(N,sizeof(double));								/* Cumulative fitness of each individual */
			double *polypos = calloc(INITTS,sizeof(double));						/* Position of neutral mutations */
			double *polyposF = calloc(INITTS,sizeof(double));						/* Position of neutral mutations (final for polymorphism sampling) */
			unsigned int *selindvP = calloc(2*N,sizeof(unsigned int));				/* Table of derived allele states per individual (parents) */
			unsigned int *selindvO = calloc(2*N,sizeof(unsigned int));				/* Table of derived allele states per individual (offspring) */
			unsigned int **neutindvP = calloc(2*N,sizeof(unsigned int *));			/* Table of neutral markers per individual (parents) */
			unsigned int **neutindvO = calloc(2*N,sizeof(unsigned int *));			/* Table of neutral markers per individual (offspring) */
			unsigned int **neutindvF = calloc(2*N,sizeof(unsigned int *));			/* Table of neutral markers per individual (final for sampling) */
			for(a = 0; a < 2*N; a++){
				neutindvP[a] = calloc(INITTS,sizeof(unsigned int));
				neutindvO[a] = calloc(INITTS,sizeof(unsigned int));
				neutindvF[a] = calloc(INITTS,sizeof(unsigned int));
			}

			/* Reassigning ancestral polymorphism state */
			npoly = npolyB;
			popread(neutindvP, polypos, npoly, N, i, suffix);
/*
			polyprint(neutindvP, polypos, npoly, N);
			exit(1);
*/
			/* Assigning initial mutant, calculating fitness */
			initp = gsl_rng_uniform_int(r,2*N);
			*(selindvP + initp) = 1;

			if(x0 == 1.0/(2.0*N)){
				scurr = s;
				nowsel = 1;
			}else if(x0 > 1.0/(2.0*N)){
				scurr = 0;
				nowsel = 0;
			}
	/*		printf("scurr is %lf\n",scurr);	*/
			fitness(N, h, scurr, selindvP, fit, cumfit, &fitsum);

			done = 0;
			while(done != 1){

				/* Generating new population */
				generation(N,self,R,nums,selindvP,selindvO,neutindvP,neutindvO,cumfit,&fitsum,npoly,polypos,r);

				/* Neutral Mutation phase */
				addpoly(N, neutindvO, polypos, &npoly, theta, INITTS, r);

				/* Reassigning matrices */
				reassign(neutindvO, neutindvP, selindvO, selindvP, npoly, N);

				/* Garbage collection of non-polymorphic sites */
				ptrim(2*N, neutindvP, neutindvF, polypos, polyposF, npoly, &npolyT,&avpi);
				npoly = npolyT;
				reassign2(neutindvF, neutindvP, polyposF, polypos, npoly, N);

				/* Checking derived allele copies and whether it has fixed or lost */
				dcopies = sumT_UI(selindvP, 2.0*N);
				/*
				printf("Dcopies are %d\n",dcopies);
				printf("copies are %d, npoly are %d\n",dcopies,npoly);
				*/
				if(dcopies == 0){
					done = 1;
/*					printf("Selected allele lost\n");*/
				}
				if(dcopies >= (int)(2.0*N*xfix)){
					done = 1;
					afix = 1;
				}

				/* Deciding if derived allele becomes selected for, then updating fitness */
				if(((dcopies/(2.0*N)) >= x0) && nowsel == 0){
					nowsel = 1;
					scurr = s;
/*					printf("scurr is %lf\n",scurr);	*/
				}
				fitness(N, h, scurr, selindvP, fit, cumfit, &fitsum);
			}

			/* Routine to print out polymorphisms if allele fixed */
			if(afix == 1){

				printf("Selected allele fixed\n");

				/* First, remove non-polymorphic sites */
				ptrim(2*N, neutindvP, neutindvF, polypos, polyposF, npoly, &npolyT,&avpi);
				npoly = npolyT;

				/* Printing out polymorphism table */
				sprintf(Pout,"Polymorphisms/Poly_%d.dat",i + suffix);
				ofp_poly = fopen(Pout,"w");

				/* State of selected site */
				fprintf(ofp_poly,"0 ");
				for(x = 0; x < 2*N; x++){
					fprintf(ofp_poly,"%d ",*(selindvP + x));
				}
				fprintf(ofp_poly,"\n");

				/* State of neutral sites */
				for(a = 0; a < npoly; a++){
					fprintf(ofp_poly,"%lf ",*(polyposF + a));
					for(x = 0; x < 2*N; x++){
						fprintf(ofp_poly,"%d ",*((*(neutindvF + x)) + a));
					}
					fprintf(ofp_poly,"\n");
				}
				fclose(ofp_poly);

				/* Then, sample from this table to produce pseudo-coalescent results */
				mutsamp(N, samps, selindvP, polyposF, neutindvF, nums2, npoly, i + suffix, rps, r);

			}

			/* End of run; freeing memory */
			for(a = 0; a < 2*N; a++){
				free(neutindvF[a]);
				free(neutindvO[a]);
				free(neutindvP[a]);
			}
			free(neutindvF);
			free(neutindvO);
			free(neutindvP);
			free(selindvO);
			free(selindvP);
			free(polyposF);
			free(polypos);
			free(cumfit);
			free(fit);
		}

	}

	free(nums2);
	free(nums);
	gsl_rng_free(r);
	return 0;

}	/* End of main program */

/* End of File */

max_stars_repo_name: MattHartfield/DomSelfAdapt

max_issues_repo_name: MattHartfield/DomSelfAdapt

max_forks_repo_name: MattHartfield/DomSelfAdapt


text:
#include <stdio.h>
#include <gsl/gsl_rng.h>

gsl_rng * r;  /* global generator */

int
main (void)
{
  const gsl_rng_type * T;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  printf ("generator type: %s\n", gsl_rng_name (r));
  printf ("seed = %lu\n", gsl_rng_default_seed);
  printf ("first value = %lu\n", gsl_rng_get (r));

  gsl_rng_free (r);
  return 0;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdlib.h>
#include <image_lib.h>
#include <gsl/gsl_fit.h>
#include "tz_darray.h"
#include "tz_stack_io.h"
#include "tz_stack_stat.h"
#include "tz_image_lib.h"
#include "tz_stack_threshold.h"
#include "tz_fimage_lib.h"
#include "tz_error.h"

#include "private/alignstack.c"

int main(int argc, const char *argv[])
{
  const char *file1 = "/Users/zhaot/Data/nathan/2008-04-18/tiletest_R1_GR1_B1_L312.tif";
  const char *file2 = "/Users/zhaot/Data/nathan/2008-04-18/tiletest_R1_GR1_B1_L334.tif";

  //const char *file1 = "/Volumes/Mac HD 2/JF8904-PRV/11/11_1.tif";
  //const char *file2 = "/Volumes/Mac HD 2/JF8904-PRV/11/11_2.tif";

  if (argc == 3) {
    file1 = argv[1];
    file2 = argv[2];
  }

  Stack *stack1 = Read_Stack((char *) file1);
  Stack *stack2 = Read_Stack((char *) file2);

#if 0
  FMatrix *filter = Mexihat_3D_F(3, NULL);
  //FMatrix_Negative(filter);

  //double sigma[] = {2.0, 2.0, 3.0};
  //FMatrix *filter = Gaussian_3D_Filter_F(sigma, NULL);

  FMatrix *out = Filter_Stack_Fast_F(stack1, filter, NULL, 0);
  Kill_Stack(stack1);

  dim_type sub[3];
  printf("%g\n", FMatrix_Max(out, sub));

  stack1 = Scale_Float_Stack(out->array, out->dim[0], out->dim[1], out->dim[2],
			      GREY);
  Kill_FMatrix(out);
  //Stack_Threshold_RC(stack1, 0, 65535);

  out = Filter_Stack_Fast_F(stack2, filter, NULL, 0);
  Kill_Stack(stack2);
  stack2 = Scale_Float_Stack(out->array, out->dim[0], out->dim[1], out->dim[2],
			      GREY);
  Kill_FMatrix(out);
  //Stack_Threshold_RC(stack2, 0, 65535);
#endif

#if 1
  int chopoff1 = estimate_chopoff(stack1);
  printf("%d\n", chopoff1);
  Stack *substack1 = Crop_Stack(stack1, 0, 0, chopoff1,
				stack1->width, stack1->height,
				stack1->depth - chopoff1, NULL);
  int chopoff2 = estimate_chopoff(stack2);
  printf("%d\n", chopoff2);
  Stack *substack2 = Crop_Stack(stack2, 0, 0, chopoff2,
				stack2->width, stack2->height,
				stack2->depth - chopoff2, NULL);
#endif

#if 0
  Stack_Threshold_Common(substack1, 0, 65535);
  Stack_Threshold_RC(substack1, 0, 65535);
  Stack_Threshold_Common(substack2, 0, 65535);
  Stack_Threshold_RC(substack2, 0, 65535);
#endif

#if 0
  Stack_Threshold_RC(substack1, 1000, 65535);
  Stack_Threshold_RC(substack2, 1000, 65535);
#endif

  tic();
  //Write_Stack("../data/test.tif", substack1);

  float unnorm_maxcorr;
  int offset[3];
  //float score = Align_Stack_F(substack1, substack2, offset, &unnorm_maxcorr);

  int intv[] = {3, 3, 3};

  float score = Align_Stack_MR_F(substack1, substack2, intv, 1, offset,
				 &unnorm_maxcorr);

  printf("%d %d %d\n", offset[0] - stack1->width + 1, offset[1]  - stack1->height + 1, offset[2] - stack1->depth + 1);

  //printf("fixed (RC_common)\n");
  printf("%lld\n", toc());

  Kill_Stack(stack1);
  Kill_Stack(stack2);

  return 0;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include "parmt_postProcess.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "compearth.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"
#include "iscl/os/os.h"

static void getBaseAndExp(const double val, double *base, int *exp);
static void setFillColor(const int i, const int iopt, char color[32]);

/*!
 * @brief Writes the global station distribution of stations, the
 *        station names, the epicenter or moment tensor, and, if
 *        desired, some indication of station polarity.
 */
int postmt_gmtHelper_writeGlobalMap(
    struct globalMapOpts_struct globalMap,
    const int nobs, const struct sacData_struct *data)
{
    const char *fcnm = "postmt_gmtHelper_writeGlobalMap\0";
    FILE *ofl;
    char *dirName, line[256], cpick[8];
    int *pol, i, ierr, l, nw, nwd, nwu;
    size_t lenos;
    const char *forwardSlash = "/";
    const int nTimeVars = 11;
    const enum sacHeader_enum timeVarNames[11]
       = {SAC_CHAR_KA,
          SAC_CHAR_KT0, SAC_CHAR_KT1, SAC_CHAR_KT2, SAC_CHAR_KT3,
          SAC_CHAR_KT4, SAC_CHAR_KT5, SAC_CHAR_KT6, SAC_CHAR_KT7,
          SAC_CHAR_KT8, SAC_CHAR_KT9};

    dirName = os_dirname(globalMap.outputScript, &ierr);
    if (!os_path_isdir(dirName))
    {
        ierr = os_makedirs(dirName);
        if (ierr != 0)
        {
            printf("%s: Failed to make output directory: %s\n", fcnm, dirName);
            return -1;
        }
    }
    memory_free8c(&dirName);

    ofl = fopen(globalMap.outputScript, "w");
    fprintf(ofl, "#!/bin/bash\n");
    fprintf(ofl, "outps=%s\n", globalMap.psFile);
    fprintf(ofl, "olat=%f\n", globalMap.evla);
    fprintf(ofl, "olon=%f\n", globalMap.evlo);
    fprintf(ofl, "J=-JH${olon}%s6i\n", forwardSlash);
    fprintf(ofl, "R=-Rg\n");
    fprintf(ofl, "gmt pscoast $J $R -B0g30 -Di -Ggray -P -K > ${outps}\n");
    fprintf(ofl, "ts=0.15i\n");

    fprintf(ofl, "# Draw great circle arcs between source and receivers\n");
    fprintf(ofl, "gmt psxy $J $R -W1p -O -K << EOF >> ${outps}\n");
    for (i=0; i<nobs; i++)
    {
        fprintf(ofl, "%8.3f %8.3f\n", data[i].header.stlo, data[i].header.stla);
        fprintf(ofl, "%8.3f %8.3f\n", globalMap.evlo, globalMap.evla);
        if (i < nobs - 1){fprintf(ofl, "\n");}
    }
    fprintf(ofl, "EOF\n");

    fprintf(ofl, "# Plot the moment tensor\n");
    postmt_gmtHelper_makePsmecaLine(globalMap.basis, globalMap.mts,
                                    globalMap.evla, globalMap.evlo,
                                    globalMap.evdp, "Optimum\0", line);
    if (globalMap.lwantMT)
    {
        fprintf(ofl, "gmt psmeca $J $R -Sm0.15i -M -N -Gblue -W0.5p,black -O -K << EOF >> ${outps}\n");
        fprintf(ofl, "%s", line);
        fprintf(ofl, "EOF\n");
    }
    else
    {
        fprintf(ofl, "# Plot the epicenter\n");
        fprintf(ofl, "gmt psxy $J $R -Sa${ts} -Gblue -Wblack -O -K << EOF >> ${outps}\n");
        fprintf(ofl, "%8.3f %8.3f\n", globalMap.evlo, globalMap.evla);
        fprintf(ofl, "EOF\n");
    }

    if (!globalMap.lwantPolarity)
    {
        fprintf(ofl, "# Plot the stations\n");
        fprintf(ofl, "gmt psxy $J $R -St${ts} -Gred -Wblack -O -K << EOF >> ${outps}\n");
        for (i=0; i<nobs; i++)
        {
            fprintf(ofl, "%8.3f %8.3f %s\n", data[i].header.stlo,
                    data[i].header.stla, data[i].header.kstnm);
        }
        fprintf(ofl, "EOF\n");
    }
    else
    {
        // Get the polarities
        nw = 0;
        nwd = 0;
        nwu = 0;
        pol = (int *) calloc((size_t) nobs, sizeof(int));
        for (i=0; i<nobs; i++)
        {
            memset(cpick, 0, 8*sizeof(char));
            for (l=0; l<nTimeVars; l++)
            {
                ierr = sacio_getCharacterHeader(timeVarNames[l],
                                                data[i].header, cpick);
                if (ierr == 0){break;}
            }
            lenos = strlen(cpick);
            if (lenos > 0)
            {
                if (cpick[lenos-1] == '+')
                {
                    nwu = nwu + 1;
                    pol[i] = 1;
                }
                else if (cpick[lenos-1] == '-')
                {
                    nwd = nwd + 1;
                    pol[i] =-1;
                }
                else
                {
                    nw = nw + 1;
                }
            }
            else
            {
                nw = nw + 1;
            }
        }
        // Write the unknowns
        if (nw > 0)
        {
            fprintf(ofl, "# Plot the indeterminant stations\n");
            fprintf(ofl, "gmt psxy $J $R -Ss${ts} -Gred -Wblack -O -K << EOF >> ${outps}\n");
            for (i=0; i<nobs; i++)
            {
                if (pol[i] == 0)
                {
                    fprintf(ofl, "%8.3f %8.3f %s\n", data[i].header.stlo,
                            data[i].header.stla, data[i].header.kstnm);
                }
            }
            fprintf(ofl, "EOF\n");
        }
        // Write the ups
        if (nwu > 0)
        {
            fprintf(ofl, "# Plot the upward stations\n");
            fprintf(ofl, "gmt psxy $J $R -St${ts} -Gred -Wblack -O -K << EOF >> ${outps}\n");
            for (i=0; i<nobs; i++)
            {
                if (pol[i] == 1)
                {
                    fprintf(ofl, "%8.3f %8.3f %s\n", data[i].header.stlo,
                            data[i].header.stla, data[i].header.kstnm);
                }
            }
            fprintf(ofl, "EOF\n");
        }
        // Write the downs
        if (nwd > 0)
        {
            fprintf(ofl, "# Plot the down stations\n");
            fprintf(ofl, "gmt psxy $J $R -Si${ts} -Gred -Wblack -O -K << EOF >> ${outps}\n");
            for (i=0; i<nobs; i++)
            {
                if (pol[i] ==-1)
                {
                    fprintf(ofl, "%8.3f %8.3f %s\n", data[i].header.stlo,
                            data[i].header.stla, data[i].header.kstnm);
                }
            }
            fprintf(ofl, "EOF\n");
        }
        free(pol);
    }

    fprintf(ofl, "# Plot the station names\n");
    fprintf(ofl, "gmt pstext $J $R -F+f8p+a0+jCT -D0%s-4p -O << EOF >> ${outps}\n",
            forwardSlash);
    for (i=0; i<nobs; i++)
    {
        fprintf(ofl, "%8.3f %8.3f %s\n", data[i].header.stlo,
                data[i].header.stla, data[i].header.kstnm);
    }
    fprintf(ofl, "EOF\n");
    fprintf(ofl, "gmt psconvert -A -Tj ${outps}\n");
    fprintf(ofl, "rm ${outps}\n");
    fclose(ofl);
    chmod(globalMap.outputScript, 0755);
    return 0;
}
//============================================================================//
int postmt_gmtHelper_writeThetaBoxes(const bool lappend, const bool lclose,
                                     const char *outputScript,
                                     const char *psFile,
                                     const int nt,
                                     const int ithetaOpt,
                                     const double *__restrict__ thetas,
                                     const double *__restrict__ thetaHist)
{
    FILE *ofl;
    char color[32], app[8], more[8], shift[8];
    double *cumTheta, *h, hAvg, thetaAvg, tmax;
    int i, ierr;
    const bool lwritePrior = true;
    tmax = array_max64f(nt, thetaHist, &ierr);
    // Compute h
    h = memory_calloc64f(nt);
    compearth_theta2h(nt, thetas, h);
    thetaAvg = -30.0;
    memset(app, 0, 8*sizeof(char));
    memset(more, 0, 8*sizeof(char));
    memset(shift, 0, 8*sizeof(char));
    if (!lappend)
    {
        ofl = fopen(outputScript, "w");
        fprintf(ofl, "#!/bin/bash\n");
        /*
        fprintf(ofl, "gmt gmtset FONT_LABEL 12p\n");
        fprintf(ofl, "gmt gmtset MAP_LABEL_OFFSET 0.1c\n");
        */
        fprintf(ofl, "parms=\"--FONT_LABEL=10p --MAP_LABEL_OFFSET=0.1c --PROJ_LENGTH_UNIT=cm\"\n");
        fprintf(ofl, "psfl=%s\n", psFile);
    }
    else
    {
        ofl = fopen(outputScript, "a");
        strcpy(app, "-O\0");
        strcpy(more, ">\0");
        strcpy(shift, "-Y4.0\0");
    }
    fprintf(ofl,
/*
            "gmt psbasemap -JX5i/1i %s -R0/90/0/%.2f -Bg10a10:\"Dips (deg)\":/a%.2f:\"Likelihood\":WSn -P %s -K >%s ${psfl}\n",
*/
            "gmt psbasemap -JX5i/1i %s -R0/90/0/%.2f -Bpxg10a10+l:\"Dips (deg)\" -Bpya%.2f+l:\"Likelihood\" -BWSnn -P %s -K ${parms} >%s ${psfl}\n",
            shift, tmax*1.1, tmax*0.2, app, more);
    setFillColor(0, ithetaOpt, color);
    if (nt > 1)
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
    }
    else
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
    }
    fprintf(ofl, "%f %f\n", 0.0, 0.0);
    fprintf(ofl, "%f %f\n", 0.0, thetaHist[0]);
    for (i=0; i<nt-1; i++)
    {
        hAvg = 0.5*(h[i] + h[i+1]);
        compearth_h2theta(1, &hAvg, &thetaAvg);
        fprintf(ofl, "%f %f\n", thetaAvg*180.0/M_PI, thetaHist[i]);
        fprintf(ofl, "%f %f\n", thetaAvg*180.0/M_PI, 0.0);
        fprintf(ofl, "EOF\n");
        setFillColor(i+1, ithetaOpt, color);
        if (i < nt - 2 || lwritePrior)
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
        }
        else
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
        }
        fprintf(ofl, "%f %f\n", thetaAvg*180.0/M_PI, 0.0);
        fprintf(ofl, "%f %f\n", thetaAvg*180.0/M_PI, thetaHist[i+1]);
    }
    fprintf(ofl, "%f %f\n", 90.0, thetaHist[nt-1]);
    fprintf(ofl, "%f %f\n", 90.0, 0.0);
    fprintf(ofl, "EOF\n");
    // Write the prior distribution
    if (lwritePrior)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,black -O -K << EOF >> ${psfl}\n");
        fprintf(ofl, "%f %f\n", 00.0, 1.0/(double) nt);
        fprintf(ofl, "%f %f\n", 90.0, 1.0/(double) nt);
        fprintf(ofl, "EOF\n");
    }
    // Write the CDF
/*
    fprintf(ofl, "gmt psbasemap -R0/90/0/1.05 -J -Bp0.2/a0.2:\"CDF\":E -O -K >> ${psfl}\n");
*/
    fprintf(ofl, "gmt psbasemap -R0/90/0/1.05 -J -Bpx0.2 -Bpya0.2+l\"CDF\" -BE -O -K ${parms} >> ${psfl}\n");
    cumTheta = array_cumsum64f(nt, thetaHist, &ierr);
    if (lclose)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O << EOF >> ${psfl}\n");
    }
    else
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O -K << EOF >> ${psfl}\n");
    }
    for (i=0; i<nt; i++)
    {
        fprintf(ofl, "%f %f\n", thetas[i]*180/M_PI, cumTheta[i]);
    }
    fprintf(ofl, "EOF\n");
    if (lclose)
    {
        //fprintf(ofl, "gmt psconvert -A -Tj ${psfl}\n");
        //fprintf(ofl, "rm ${psfl}\n");
        fprintf(ofl, "gmt psconvert -A -Tf ${psfl}\n");
        fprintf(ofl, "/bin/rm ${psfl}\n");
    }
    memory_free64f(&cumTheta);
    memory_free64f(&h);
    fclose(ofl);
    chmod(outputScript, 0755);
    return 0;
}
//============================================================================//
int postmt_gmtHelper_writeSigmaBoxes(const bool lappend, const bool lclose,
                                     const char *outputScript,
                                     const char *psFile,
                                     const int ns,
                                     const int isigmaOpt,
                                     const double *__restrict__ sigmas,
                                     const double *__restrict__ sigmaHist)
{
    FILE *ofl;
    char color[32], app[8], more[8], shift[8];
    double *cumSigma, ds, smax;
    int i, ierr;
    const bool lwritePrior = true;
    smax = array_max64f(ns, sigmaHist, &ierr);
    // Take averages
    ds = 0.0;
    memset(app, 0, 8*sizeof(char));
    memset(more, 0, 8*sizeof(char));
    memset(shift, 0, 8*sizeof(char));
    if (!lappend)
    {
        ofl = fopen(outputScript, "w");
        fprintf(ofl, "#!/bin/bash\n");
        /*
        fprintf(ofl, "gmt gmtset FONT_LABEL 12p\n");
        fprintf(ofl, "gmt gmtset MAP_LABEL_OFFSET 0.1c\n");
        */
        fprintf(ofl, "parms=\"--FONT_LABEL=10p --MAP_LABEL_OFFSET=0.1c --PROJ_LENGTH_UNIT=cm\"\n");
        fprintf(ofl, "psfl=%s\n", psFile);
    }
    else
    {
        ofl = fopen(outputScript, "a");
        strcpy(app, "-O\0");
        strcpy(more, ">\0");
        strcpy(shift, "-Y4.0\0");
    }
    fprintf(ofl,
/*
           "gmt psbasemap -JX5i/1i %s -R-90/90/0/%.2f -Bg15a15:\"Slips (deg)\":/a%.2f:\"Likelihood\":WSn -P %s -K >%s ${psfl}\n",
*/
            "gmt psbasemap -JX5i/1i %s -R-90/90/0/%.2f -Bpxg15a15+l\"Slips (deg)\" -Bpya%.2f+l\"Likelihood\" -BWSn -P %s -K ${parms} >%s ${psfl}\n",
            shift, smax*1.1, smax*0.2, app, more);
    setFillColor(0, isigmaOpt, color);
    if (ns > 1)
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
       ds = (sigmas[1] - sigmas[0])*180.0/M_PI;
    }
    else
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
    }
    fprintf(ofl, "%.2f %f\n", -90.0, 0.0);
    fprintf(ofl, "%.2f %f\n", -90.0, sigmaHist[0]);
    for (i=0; i<ns-1; i++)
    {
        fprintf(ofl, "%.2f %f\n", -90.0 + (double) (i+1)*ds, sigmaHist[i]);
        fprintf(ofl, "%.2f %f\n", -90.0 + (double) (i+1)*ds, 0.0);
        fprintf(ofl, "EOF\n");
        setFillColor(i+1, isigmaOpt, color);
        if (i < ns - 2 || lwritePrior)
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
        }
        else
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
        }
        fprintf(ofl, "%.2f %f\n", -90.0 + (double) (i+1)*ds, 0.0);
        fprintf(ofl, "%.2f %f\n", -90.0 + (double) (i+1)*ds, sigmaHist[i+1]);
    }
    fprintf(ofl, "%.2f %f\n", 90.0, sigmaHist[ns-1]);
    fprintf(ofl, "%.2f %f\n", 90.0, 0.0);
    fprintf(ofl, "EOF\n");
    // Write the prior distribution
    if (lwritePrior)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,black -O -K << EOF >> ${psfl}\n");
        fprintf(ofl, "%f %f\n",-90.0, 1.0/(double) ns);
        fprintf(ofl, "%f %f\n", 90.0, 1.0/(double) ns);
        fprintf(ofl, "EOF\n");
    }
    // Write the CDF
/*
    fprintf(ofl, "gmt psbasemap -R-90/90/0/1.05 -J -Bp0.2/a0.2:\"CDF\":E -O -K >> ${psfl}\n");
*/
    fprintf(ofl, "gmt psbasemap -R-90/90/0/1.05 -J -Bpx0.2 -Bpya0.2+l\"CDF\" -BE -O -K ${parms} >> ${psfl}\n");
    cumSigma = array_cumsum64f(ns, sigmaHist, &ierr);
    if (lclose)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O << EOF >> ${psfl}\n");
    }
    else
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O -K << EOF >> ${psfl}\n");
    }
    for (i=0; i<ns; i++)
    {
        fprintf(ofl, "%f %f\n", sigmas[i]*180/M_PI, cumSigma[i]);
    }
    fprintf(ofl, "EOF\n");
    if (lclose)
    {
        //fprintf(ofl, "gmt psconvert -A -Tj ${psfl}\n");
        //fprintf(ofl, "rm ${psfl}\n");
        fprintf(ofl, "gmt psconvert -A -Tf ${psfl}\n");
        fprintf(ofl, "/bin/rm ${psfl}\n");
    }
    memory_free64f(&cumSigma);
    fclose(ofl);
    chmod(outputScript, 0755);
    return 0;
}

//============================================================================//
int postmt_gmtHelper_writeKappaBoxes(const bool lappend, const bool lclose,
                                     const char *outputScript,
                                     const char *psFile,
                                     const int nk,
                                     const int kappaOpt,
                                     const double *__restrict__ kappas,
                                     const double *__restrict__ kappaHist)
{
    FILE *ofl;
    char color[32], app[8], more[8], shift[8];
    double *cumKappa, dk, kmax;
    int i, ierr;
    const bool lwritePrior = true;
    kmax = array_max64f(nk, kappaHist, &ierr);
    // Take averages
    dk = 0.0;
    memset(app, 0, 8*sizeof(char));
    memset(more, 0, 8*sizeof(char));
    memset(shift, 0, 8*sizeof(char));
    if (!lappend)
    {
        ofl = fopen(outputScript, "w");
        fprintf(ofl, "#!/bin/bash\n");
        /*
        fprintf(ofl, "gmt gmtset FONT_LABEL 12p\n");
        fprintf(ofl, "gmt gmtset MAP_LABEL_OFFSET 0.1c\n");
        */
        fprintf(ofl, "parms=\"--FONT_LABEL=10p --MAP_LABEL_OFFSET=0.1c --PROJ_LENGTH_UNIT=cm\"\n");
        fprintf(ofl, "psfl=%s\n", psFile);
    }
    else
    {
        ofl = fopen(outputScript, "a");
        strcpy(app, "-O\0");
        strcpy(more, ">\0");
        strcpy(shift, "-Y4.0\0");
    }
    fprintf(ofl,
/*
            "gmt psbasemap -JX5i/1i %s -R0/360/0/%.2f -Bg30a30:\"Strike (deg)\":/a%.2f:\"Likelihood\":WSn -P %s -K >%s ${psfl}\n",
*/
            "gmt psbasemap -JX5i/1i %s -R0/360/0/%.2f -Bpxg30a30+l\"Strike (deg)\" -Bpya%.2f+l\"Likelihood\" -BWSn -P %s -K ${parms} >%s ${psfl}\n",
            shift, kmax*1.1, kmax*0.2, app, more);
    setFillColor(0, kappaOpt, color);
    if (nk > 1)
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
       dk = (kappas[1] - kappas[0])*180.0/M_PI;
    }
    else
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
    }
    fprintf(ofl, "%.2f %f\n", 0.0, 0.0);
    fprintf(ofl, "%.2f %f\n", 0.0, kappaHist[0]);
    for (i=0; i<nk-1; i++)
    {
        fprintf(ofl, "%.2f %f\n", 0.0 + (double) (i+1)*dk, kappaHist[i]);
        fprintf(ofl, "%.2f %f\n", 0.0 + (double) (i+1)*dk, 0.0);
        fprintf(ofl, "EOF\n");
        setFillColor(i+1, kappaOpt, color);
        if (i < nk - 2 || lwritePrior)
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
        }
        else
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
        }
        fprintf(ofl, "%.2f %f\n", 0.0 + (double) (i+1)*dk, 0.0);
        fprintf(ofl, "%.2f %f\n", 0.0 + (double) (i+1)*dk, kappaHist[i+1]);
    }
    fprintf(ofl, "%.2f %f\n", 360.0, kappaHist[nk-1]);
    fprintf(ofl, "%.2f %f\n", 360.0, 0.0);
    fprintf(ofl, "EOF\n");
    // Write the prior distribution
    if (lwritePrior)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,black -O -K << EOF >> ${psfl}\n");
        fprintf(ofl, "%f %f\n", 0.0,   1.0/(double) nk);
        fprintf(ofl, "%f %f\n", 360.0, 1.0/(double) nk);
        fprintf(ofl, "EOF\n");
    }
    // Write the CDF
/*
    fprintf(ofl, "gmt psbasemap -R0/360/0/1.05 -J -Bp0.2/a0.2:\"CDF\":E -O -K >> ${psfl}\n");
*/
    fprintf(ofl, "gmt psbasemap -R0/360/0/1.05 -J -Bpx0.2 -Bpya0.2+l\"CDF\" -BE -O -K ${parms} >> ${psfl}\n");
    cumKappa = array_cumsum64f(nk, kappaHist, &ierr);
    if (lclose)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O << EOF >> ${psfl}\n");
    }
    else
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O -K << EOF >> ${psfl}\n");
    }
    for (i=0; i<nk; i++)
    {
        fprintf(ofl, "%f %f\n", kappas[i]*180.0/M_PI, cumKappa[i]);
    }
    fprintf(ofl, "EOF\n");
    if (lclose)
    {
        //fprintf(ofl, "gmt psconvert -A -Tj ${psfl}\n");
        //fprintf(ofl, "rm ${psfl}\n");
        fprintf(ofl, "gmt psconvert -A -Tf ${psfl}\n");
        fprintf(ofl, "/bin/rm ${psfl}\n");
    }
    memory_free64f(&cumKappa);
    fclose(ofl);
    chmod(outputScript, 0755);
    return 0;
}
//============================================================================//
int postmt_gmtHelper_writeDepthBoxes(const bool lappend, const bool lclose,
                                     const char *outputScript,
                                     const char *psFile,
                                     const int nd,
                                     const int idepOpt,
                                     const double *__restrict__ deps,
                                     const double *__restrict__ depHist)
{
    FILE *ofl;
    char color[32], app[8], more[8], shift[8];
    double *cumDep, dd, dmax, depMin, depMax;
    int i, ierr;
    const bool lwritePrior = true;
    // Compute moment magnitudes
    dmax = array_max64f(nd, depHist, &ierr);
    dd = 0.1;
    if (nd > 1){dd = deps[1] - deps[0];};
    depMin = fmax(0.0, deps[0]    - dd/2.0);
    depMax = fmax(0.0, deps[nd-1] + dd/2.0);
    memset(app, 0, 8*sizeof(char));
    memset(more, 0, 8*sizeof(char));
    memset(shift, 0, 8*sizeof(char));
    if (!lappend)
    {
        ofl = fopen(outputScript, "w");
        fprintf(ofl, "#!/bin/bash\n");
        /*
        fprintf(ofl, "gmt gmtset FONT_LABEL 12p\n");
        fprintf(ofl, "gmt gmtset MAP_LABEL_OFFSET 0.1c\n");
        */
        fprintf(ofl, "psfl=%s\n", psFile);
    }
    else
    {
        ofl = fopen(outputScript, "a");
        strcpy(app, "-O\0");
        strcpy(more, ">\0");
        strcpy(shift, "-Y4.0\0");
    }
    fprintf(ofl,
/*
            "gmt psbasemap -JX5i/1i %s -R%.1f/%.1f/0/%.2f -Bg%fa%f:\"Depths (km)\":/a%.2f:\"Likelihood\":WSn -P %s -K >%s ${psfl}\n",
*/
            "gmt psbasemap -JX5i/1i %s -R%.1f/%.1f/0/%.2f -Bpxg%fa%f+l\"Depths (km)\" -Bpya%.2f+l\"Likelihood\" -BWSn -P %s -K ${parms} >%s ${psfl}\n",
            shift, depMin, depMax, dmax*1.1, dd, (nd-1)*dd/5.0, dmax*0.2, app, more);
    setFillColor(0, idepOpt, color);
    if (nd > 1)
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
    }
    else
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
    }
    fprintf(ofl, "%.2f %f\n", depMin, 0.0);
    fprintf(ofl, "%.2f %f\n", depMin, depHist[0]);
    for (i=0; i<nd-1; i++)
    {
        fprintf(ofl, "%.2f %f\n", depMin + (double) (i+1)*dd, depHist[i]);
        fprintf(ofl, "%.2f %f\n", depMin + (double) (i+1)*dd, 0.0);
        fprintf(ofl, "EOF\n");
        setFillColor(i+1, idepOpt, color);
        if (i < nd - 2 || lwritePrior)
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
        }
        else
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
        }
        fprintf(ofl, "%.2f %f\n", depMin + (double) (i+1)*dd, 0.0);
        fprintf(ofl, "%.2f %f\n", depMin + (double) (i+1)*dd, depHist[i+1]);
    }
    fprintf(ofl, "%.2f %f\n", depMax, depHist[nd-1]);
    fprintf(ofl, "%.2f %f\n", depMax, 0.0);
    fprintf(ofl, "EOF\n");
    // Write the prior distribution
    if (lwritePrior)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,black -O -K << EOF >> ${psfl}\n");
        fprintf(ofl, "%f %f\n", depMin, 1.0/(double) nd);
        fprintf(ofl, "%f %f\n", depMax, 1.0/(double) nd);
        fprintf(ofl, "EOF\n");
    }
    // Write the CDF
/*
    fprintf(ofl, "gmt psbasemap -R%f/%f/0/1.05 -J -Bp0.2/a0.2:\"CDF\":E -O -K >> ${psfl}\n", depMin, depMax);
*/
    fprintf(ofl, "gmt psbasemap -R%f/%f/0/1.05 -J -Bpx0.2 -Bpya0.2+l\"CDF\" -BE -O -K ${parms} >> ${psfl}\n", depMin, depMax);
    cumDep = array_cumsum64f(nd, depHist, &ierr);
    if (lclose)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O << EOF >> ${psfl}\n");
    }
    else
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O -K << EOF >> ${psfl}\n");
    }
    for (i=0; i<nd; i++)
    {
        fprintf(ofl, "%f %f\n", deps[i], cumDep[i]);
    }
    fprintf(ofl, "EOF\n");
    if (lclose)
    {
        //fprintf(ofl, "gmt psconvert -A -Tj ${psfl}\n");
        //fprintf(ofl, "rm ${psfl}\n");
        fprintf(ofl, "gmt psconvert -A -Tf ${psfl}\n");
        fprintf(ofl, "/bin/rm ${psfl}\n");
    }
    memory_free64f(&cumDep);
    fclose(ofl);
    chmod(outputScript, 0755);
    return 0;
}
//============================================================================//
int postmt_gmtHelper_writeMagnitudeBoxes(const bool lappend, const bool lclose,
                                         const char *outputScript,
                                         const char *psFile,
                                         const int nm,
                                         const int magOpt,
                                         const double *__restrict__ M0s,
                                         const double *__restrict__ magHist)
{
    FILE *ofl;
    char color[32], app[8], more[8], shift[8];
    double *cumMag, *Mw, dm, mmax, mwMin, mwMax;
    int i, ierr;
    const bool lwritePrior = true;
    // Compute moment magnitudes
    Mw = memory_calloc64f(nm);
    mmax = array_max64f(nm, magHist, &ierr);
    compearth_m02mw(nm, CE_KANAMORI_1978, M0s, Mw);
    dm = 0.1;
    if (nm > 1){dm = Mw[1] - Mw[0];};
    mwMin = Mw[0] - dm/2.0;
    mwMax = Mw[nm-1] + dm/2.0;
    memset(app, 0, 8*sizeof(char));
    memset(more, 0, 8*sizeof(char));
    memset(shift, 0, 8*sizeof(char));
    if (!lappend)
    {
        ofl = fopen(outputScript, "w");
        fprintf(ofl, "#!/bin/bash\n");
        /*
        fprintf(ofl, "gmt gmtset FONT_LABEL 12p\n");
        fprintf(ofl, "gmt gmtset MAP_LABEL_OFFSET 0.1c\n");
        */
        fprintf(ofl, "parms=\"--FONT_LABEL=10p --MAP_LABEL_OFFSET=0.1c --PROJ_LENGTH_UNIT=cm\"\n");
        fprintf(ofl, "psfl=%s\n", psFile);
    }
    else
    {
        ofl = fopen(outputScript, "a");
        strcpy(app, "-O\0");
        strcpy(more, ">\0");
        strcpy(shift, "-Y4.0\0");
    }
    fprintf(ofl,
/*
            "gmt psbasemap -JX5i/1i %s -R%.2f/%.2f/0/%.2f -Bg%fa%f:\"Magnitude (Mw)\":/a%.2f:\"Likelihood\":WSn -P %s -K >%s ${psfl}\n",
*/
            "gmt psbasemap -JX5i/1i %s -R%.2f/%.2f/0/%.2f -Bpxg%fa%f+l\"Magnitude (Mw)\" -Bpya%.2f+l\"Likelihood\" -BWSn -P %s -K ${parms} >%s ${psfl}\n",
            shift, mwMin, mwMax, mmax*1.1, dm, (nm-1)*dm/5.0, mmax*0.2, app, more);
    setFillColor(0, magOpt, color);
    if (nm > 1)
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
    }
    else
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
    }
    fprintf(ofl, "%.2f %f\n", mwMin, 0.0);
    fprintf(ofl, "%.2f %f\n", mwMin, magHist[0]);
    for (i=0; i<nm-1; i++)
    {
        fprintf(ofl, "%.2f %f\n", mwMin + (double) (i+1)*dm, magHist[i]);
        fprintf(ofl, "%.2f %f\n", mwMin + (double) (i+1)*dm, 0.0);
        fprintf(ofl, "EOF\n");
        setFillColor(i+1, magOpt, color);
        if (i < nm - 2 || lwritePrior)
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
        }
        else
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
        }
        fprintf(ofl, "%.2f %f\n", mwMin + (double) (i+1)*dm, 0.0);
        fprintf(ofl, "%.2f %f\n", mwMin + (double) (i+1)*dm, magHist[i+1]);
    }
    fprintf(ofl, "%.2f %f\n", mwMax, magHist[nm-1]);
    fprintf(ofl, "%.2f %f\n", mwMax, 0.0);
    fprintf(ofl, "EOF\n");
    // Write the prior distribution
    if (lwritePrior)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,black -O -K << EOF >> ${psfl}\n");
        fprintf(ofl, "%f %f\n", mwMin, 1.0/(double) nm);
        fprintf(ofl, "%f %f\n", mwMax, 1.0/(double) nm);
        fprintf(ofl, "EOF\n");
    }
    // Write the CDF
/*
    fprintf(ofl, "gmt psbasemap -R%f/%f/0/1.05 -J -Bp0.2/a0.2:\"CDF\":E -O -K >> ${psfl}\n", mwMin, mwMax);
*/
    fprintf(ofl, "gmt psbasemap -R%f/%f/0/1.05 -J -Bpx0.2 -Bpya0.2+l\"CDF\" -BE -O -K ${parms} >> ${psfl}\n", mwMin, mwMax);
    cumMag = array_cumsum64f(nm, magHist, &ierr);
    if (lclose)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O << EOF >> ${psfl}\n");
    }
    else
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O -K << EOF >> ${psfl}\n");
    }
    for (i=0; i<nm; i++)
    {
        fprintf(ofl, "%f %f\n", Mw[i], cumMag[i]);
    }
    fprintf(ofl, "EOF\n");
    if (lclose)
    {
        //fprintf(ofl, "gmt psconvert -A -Tj ${psfl}\n");
        //fprintf(ofl, "rm ${psfl}\n");
    }
    fclose(ofl);
    memory_free64f(&cumMag);
    memory_free64f(&Mw);
    chmod(outputScript, 0755);
    return 0;
}
//============================================================================//
int postmt_gmtHelper_writeGammaBoxes(const bool lappend, const bool lclose,
                                     const char *outputScript,
                                     const char *psFile,
                                     const int ng,
                                     const int igammaOpt,
                                     const double *__restrict__ gammas,
                                     const double *__restrict__ gammaHist)
{
    FILE *ofl;
    char color[32], app[8], more[8], shift[8];
    double *cumGamma, *v, gmax, gammaAvg, vAvg;
    int i, ierr;
    const bool lwritePrior = true;
    gmax = array_max64f(ng, gammaHist, &ierr);
    // Compute v
    v = memory_calloc64f(ng);
    compearth_gamma2v(ng, gammas, v);
    // Take averages
    gammaAvg = -30.0;
    memset(app, 0, 8*sizeof(char));
    memset(shift, 0, 8*sizeof(char));
    memset(more, 0, 8*sizeof(char));
    if (!lappend)
    {
        ofl = fopen(outputScript, "w");
        fprintf(ofl, "#!/bin/bash\n");
        /*
        fprintf(ofl, "gmt gmtset FONT_LABEL 12p\n");
        fprintf(ofl, "gmt gmtset MAP_LABEL_OFFSET 0.1c\n");
        */
        fprintf(ofl, "parms=\"--FONT_LABEL=10p --MAP_LABEL_OFFSET=0.1c --PROJ_LENGTH_UNIT=cm\"\n");
        fprintf(ofl, "psfl=%s\n", psFile);
    }
    else
    {
        ofl = fopen(outputScript, "a");
        strcpy(app, "-O\0");
        strcpy(more, ">\0");
        strcpy(shift, "-Y4.0\0");
    }
    fprintf(ofl,
/*
            "gmt psbasemap -JX5i/1i %s -R-30/30/0/%.2f -Bg5a5:\"Longitude (deg)\":/a%.2f:\"Likelihood\":WSn -P %s -K >%s ${psfl}\n",
*/
            "gmt psbasemap -JX5i/1i %s -R-30/30/0/%.2f -Bpxg5a5+l\"Longitude (deg)\" -Bpya%.2f+l\"Likelihood\" -BWSn -P %s -K ${parms} >%s ${psfl}\n",
            shift, gmax*1.1, gmax*0.2, app, more);
    setFillColor(0, igammaOpt, color);
    if (ng > 1)
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
    }
    else
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
    }
    fprintf(ofl, "%f %f\n", -30.0, 0.0);
    fprintf(ofl, "%f %f\n", -30.0, gammaHist[0]);
    for (i=0; i<ng-1; i++)
    {
        vAvg = 0.5*(v[i] + v[i+1]);
        compearth_v2gamma(1, &vAvg, &gammaAvg);
        fprintf(ofl, "%f %f\n", gammaAvg*180.0/M_PI, gammaHist[i]);
        fprintf(ofl, "%f %f\n", gammaAvg*180.0/M_PI, 0.0);
        fprintf(ofl, "EOF\n");
        setFillColor(i+1, igammaOpt, color);
        if (i < ng - 2 || lwritePrior)
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
        }
        else
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
        }
        fprintf(ofl, "%f %f\n", gammaAvg*180.0/M_PI, 0.0);
        fprintf(ofl, "%f %f\n", gammaAvg*180.0/M_PI, gammaHist[i+1]);
    }
    fprintf(ofl, "%f %f\n", 30.0, gammaHist[ng-1]);
    fprintf(ofl, "%f %f\n", 30.0, 0.0);
    fprintf(ofl, "EOF\n");
    // Write the prior distribution
    if (lwritePrior)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,black -O -K << EOF >> ${psfl}\n");
        fprintf(ofl, "%f %f\n",-30.0, 1.0/(double) ng);
        fprintf(ofl, "%f %f\n", 30.0, 1.0/(double) ng);
        fprintf(ofl, "EOF\n");
    }
    // Write the CDF
/*
    fprintf(ofl, "gmt psbasemap -R-30/30/0/1.05 -J -Bp0.2/a0.2:\"CDF\":E -O -K >> ${psfl}\n");
*/
    fprintf(ofl, "gmt psbasemap -R-30/30/0/1.05 -J -Bpx0.2 -Bpya0.2+l\"CDF\" -BE -O -K ${parms} >> ${psfl}\n");
    cumGamma = array_cumsum64f(ng, gammaHist, &ierr);
    if (lclose)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O << EOF >> ${psfl}\n");
    }
    else
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O -K << EOF >> ${psfl}\n");
    }
    for (i=0; i<ng; i++)
    {
        fprintf(ofl, "%f %f\n", gammas[i]*180.0/M_PI, cumGamma[i]);
    }
    fprintf(ofl, "EOF\n");
    if (lclose)
    {
        //fprintf(ofl, "gmt psconvert -A -Tj ${psfl}\n");
        //fprintf(ofl, "rm ${psfl}\n");
        fprintf(ofl, "gmt psconvert -A -Tf ${psfl}\n");
        fprintf(ofl, "/bin/rm ${psfl}\n");
    }
    memory_free64f(&cumGamma);
    memory_free64f(&v);
    fclose(ofl);
    chmod(outputScript, 0755);
    return 0;
}
//============================================================================//
int postmt_gmtHelper_writeBetaBoxes(const bool lappend, const bool lclose,
                                    const char *outputScript,
                                    const char *psFile,
                                    const int nb,
                                    const int ibetaOpt,
                                    const double *__restrict__ betas,
                                    const double *__restrict__ betaHist)
{
    FILE *ofl;
    char color[32], app[8], more[8], shift[8];
    double *betaCum, *u, bmax, betaAvg, uAvg;
    int i, ierr;
    const bool lwritePrior = true;
    bmax = array_max64f(nb, betaHist, &ierr);
    // Compute u
    u = memory_calloc64f(nb);
    compearth_beta2u(nb, betas, u);
    // Take averages
    betaAvg = 0.0;
    memset(app, 0, 8*sizeof(char));
    memset(more, 0, 8*sizeof(char));
    memset(shift, 0, 8*sizeof(char));
    if (!lappend)
    {
        ofl = fopen(outputScript, "w");
        fprintf(ofl, "#!/bin/bash\n");
        /*
        fprintf(ofl, "gmt gmtset FONT_LABEL 12p\n");
        fprintf(ofl, "gmt gmtset MAP_LABEL_OFFSET 0.1c\n");
        */
        fprintf(ofl, "parms=\"--FONT_LABEL=10p --MAP_LABEL_OFFSET=0.1c --PROJ_LENGTH_UNIT=cm\"\n");
        fprintf(ofl, "psfl=%s\n", psFile);
    }
    else
    {
        ofl = fopen(outputScript, "a");
        strcpy(app, "-O\0");
        strcpy(more, ">\0");
        strcpy(shift, "-Y4.0\0");
    }
    fprintf(ofl,
/*
            "gmt psbasemap -JX5i/1i %s -R-90/90/0/%.3f -Bg15a15:\"Latitude (deg)\":/a%.2f:\"Likelihood\":WSn -P %s -K >%s ${psfl}\n",
*/
            "gmt psbasemap -JX5i/1i %s -R-90/90/0/%.3f -Bpxg15a15+l\"Latitude (deg)\" -Bpya%.2f+l\"Likelihood\" -BWSn -P %s -K ${parms} >%s ${psfl}\n",
            shift, bmax*1.1, bmax*0.2, app, more);
    setFillColor(0, ibetaOpt, color);
    if (nb > 1)
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
    }
    else
    {
       fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
    }
    fprintf(ofl, "%f %f\n", 90.0 - 0.0, 0.0);
    fprintf(ofl, "%f %f\n", 90.0 - 0.0, betaHist[0]);
    for (i=0; i<nb-1; i++)
    {
        uAvg = 0.5*(u[i] + u[i+1]);
        compearth_u2beta(1, &uAvg, &betaAvg);
        //compearth_u2beta(1, 20, 2, &uAvg, 1.e-7, &betaAvg);
        fprintf(ofl, "%f %f\n", 90.0 - betaAvg*180.0/M_PI, betaHist[i]);
        fprintf(ofl, "%f %f\n", 90.0 - betaAvg*180.0/M_PI, 0.0);
        fprintf(ofl, "EOF\n");
        setFillColor(i+1, ibetaOpt, color);
        if (i < nb - 2 || lwritePrior)
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O -K << EOF >> ${psfl}\n", color);
        }
        else
        {
            fprintf(ofl, "gmt psxy -R -J -Wblack %s -O << EOF >> ${psfl}\n", color);
        }
        fprintf(ofl, "%f %f\n", 90.0 - betaAvg*180.0/M_PI, 0.0);
        fprintf(ofl, "%f %f\n", 90.0 - betaAvg*180.0/M_PI, betaHist[i+1]);
    }
    fprintf(ofl, "%f %f\n", 90.0 - M_PI*180.0/M_PI, betaHist[nb-1]);
    fprintf(ofl, "%f %f\n", 90.0 - M_PI*180.0/M_PI, 0.0);
    fprintf(ofl, "EOF\n");
    // Write the prior distribution
    if (lwritePrior)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,black -O -K << EOF >> ${psfl}\n");
        fprintf(ofl, "%f %f\n",-90.0, 1.0/(double) nb);
        fprintf(ofl, "%f %f\n", 90.0, 1.0/(double) nb);
        fprintf(ofl, "EOF\n");
    }
    // Write the CDF
/*
    fprintf(ofl, "gmt psbasemap -R-90/90/0/1.05 -J -Bp0.2/a0.2:\"CDF\":E -O -K >> ${psfl}\n");
*/
    fprintf(ofl, "gmt psbasemap -R-90/90/0/1.05 -J -Bpx0.2 -Bpya0.2+l\"CDF\" -BE -O -K ${parms} >> ${psfl}\n");
    betaCum = array_cumsum64f(nb, betaHist, &ierr);
    if (lclose)
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O << EOF >> ${psfl}\n");
    }
    else
    {
        fprintf(ofl, "gmt psxy -R -J -W1,blue -O -K << EOF >> ${psfl}\n");
    }
    for (i=0; i<nb; i++)
    {
        // the cdf is written backwards because i flip the axis to make it
        // increase left to right; this is opposite of what colatitude wants
        // to do
        fprintf(ofl, "%f %f\n", 90.0 - betas[i]*180.0/M_PI, betaCum[nb-1-i]);
    }
    fprintf(ofl, "EOF\n");
    if (lclose)
    {
        //fprintf(ofl, "gmt psconvert -A -Tj ${psfl}\n");
        //fprintf(ofl, "rm ${psfl}\n");
        fprintf(ofl, "gmt psconvert -A -Tf ${psfl}\n");
        fprintf(ofl, "/bin/rm ${psfl}\n");
    }
    memory_free64f(&u);
    memory_free64f(&betaCum);
    fclose(ofl);
    chmod(outputScript, 0755);
    return 0;
}
//============================================================================//
int postmt_gmtHelper_makeRegularHistograms(
    const int nlocs, const int nm,
    const int nb, const int ng, const int nk,
    const int ns, const int nt,
    const int nmt, const double *__restrict__ phi,
    double **__restrict__ locHist,
    double **__restrict__ magHist,
    double **__restrict__ betaHist,
    double **__restrict__ gammaHist,
    double **__restrict__ kappaHist,
    double **__restrict__ sigmaHist,
    double **__restrict__ thetaHist)
{
    double *betaH, *gammaH, *kappaH, *locH, *magH, *sigmaH, *thetaH, xsum;
    int ib, ierr, ig, ik, il, im, imt, is, it;
    locH = memory_calloc64f(nlocs);
    magH = memory_calloc64f(nm);
    betaH = memory_calloc64f(nb);
    gammaH = memory_calloc64f(ng);
    kappaH = memory_calloc64f(nk);
    sigmaH = memory_calloc64f(ns);
    thetaH = memory_calloc64f(nt);
    for (il=0; il<nlocs; il++)
    {
        for (im=0; im<nm; im++)
        {
            for (ib=0; ib<nb; ib++)
            {
                for (ig=0; ig<ng; ig++)
                {
                    for (ik=0; ik<nk; ik++)
                    {
                        for (is=0; is<ns; is++)
                        {
                            for (it=0; it<nt; it++)
                            {
                                imt = il*nm*nb*ng*nk*ns*nt
                                    + im*nb*ng*nk*ns*nt
                                    + ib*ng*nk*ns*nt
                                    + ig*nk*ns*nt
                                    + ik*ns*nt
                                    + is*nt
                                    + it;
                                locH[il]   = locH[il]   + phi[imt];
                                magH[im]   = magH[im]   + phi[imt];
                                betaH[ib]  = betaH[ib]  + phi[imt];
                                gammaH[ig] = gammaH[ig] + phi[imt];
                                kappaH[ik] = kappaH[ik] + phi[imt];
                                sigmaH[is] = sigmaH[is] + phi[imt];
                                thetaH[it] = thetaH[it] + phi[imt];
                            }
                        }
                    }
                }
            }
        }
    }
    xsum = array_sum64f(nmt, phi, &ierr);
    cblas_dscal(nlocs, 1.0/xsum, locH, 1);
    cblas_dscal(nm, 1.0/xsum, magH, 1);
    cblas_dscal(nb, 1.0/xsum, betaH, 1);
    cblas_dscal(ng, 1.0/xsum, gammaH, 1);
    cblas_dscal(nk, 1.0/xsum, kappaH, 1);
    cblas_dscal(ns, 1.0/xsum, sigmaH, 1);
    cblas_dscal(nt, 1.0/xsum, thetaH, 1);
    *locHist = locH;
    *magHist = magH;
    *betaHist = betaH;
    *gammaHist = gammaH;
    *kappaHist = kappaH;
    *sigmaHist = sigmaH;
    *thetaHist = thetaH;
    return 0;
}
//============================================================================//
static void getBaseAndExp(const double val, double *base, int *exp)
{
    char cval[64], cexp[64], temp[64];
    int i, k, l;
    bool lp1;
    memset(cval, 0, 64*sizeof(char));
    memset(cexp, 0, 64*sizeof(char));
    memset(temp, 0, 64*sizeof(char));
    sprintf(temp, "%04.2e", val);
    lp1 = true;
    k = 0;
    l = 0;
    for (i=0; i<strlen(temp); i++)
    {
        if (temp[i] == 'e')
        {
            lp1 = false;
            continue;
        }
        else
        {
            if (lp1)
            {
                cval[k] = temp[i];
                k = k + 1;
            }
            else
            {
                cexp[l] = temp[i];
                l = l + 1;
            }
        }
    }
    *exp = atoi(cexp);
    *base = atof(cval);
    return;
}

static void getBaseAndExpMT(const double *mtIn, double *mtOut, int *expOut)
{
    double xfact;
    int expWork, i;
    getBaseAndExp(mtIn[0], &mtOut[0], expOut);
    for (i=1; i<6; i++)
    {
        getBaseAndExp(mtIn[i], &mtOut[i], &expWork);
        if (expWork > *expOut){*expOut = expWork;}
    }
    // Rescale
    for (i=0; i<6; i++)
    {
        getBaseAndExp(mtIn[i], &mtOut[i], &expWork);
        xfact = pow(10.0, expWork - *expOut);
        mtOut[i] = mtOut[i]*xfact;
    }
    return;
}

int postmt_gmtHelper_makePsmecaLine(const enum compearthCoordSystem_enum basis,
                                    const double *mt,
                                    const double evla, const double evlo,
                                    const double evdp, const char *evid,
                                    char line[128])
{
    const char *fcnm = "postmt_gmtHelper_makePsmecaLine\0";
    double mtUSE[6], mtGMT[6];
    int exp, ierr;
    memset(line, 0, 128*sizeof(char));
    ierr = compearth_convertMT(1, basis, CE_USE, mt, mtUSE);
    if (ierr != 0)
    {
        printf("%s: Error switching basis\n", fcnm);
        return -1;
    }
    getBaseAndExpMT(mtUSE, mtGMT, &exp);
    sprintf(line, "%f %f %f %f %f %f %f %f %f %d %s\n",
            evlo, evla, evdp,
            mtGMT[0], mtGMT[1], mtGMT[2], mtGMT[3], mtGMT[4], mtGMT[5],
            exp, evid);
    return 0;
}

static void setFillColor(const int i, const int iopt, char color[32])
{
    memset(color, 0, 32*sizeof(char));
    if (i == iopt)
    {
        strcpy(color, "-Gyellow\0");
    }
    else
    {
        strcpy(color, "-Gred\0");
    }
    return;
}
//============================================================================//
/*
    sprintf("%7.2f %5.2f %f \n"<
            evlo, evla, evdp,
            mrr, mtt, mpp,

Columns: lon lat depth mrr mtt mpp mrt mrp mtp iexp name
-176.96 -29.25 48 7.68 0.09 -7.77 1.39 4.52 -3.26 26 X Y 010176A
*/


text:
/*
        Version 1.0

Correlation service; initial version will accept text input (optionally gzipped)
where comments and any integer (ie containing no ".") in initial columns are
skipped:

[# comments]
[int int int ] fl.oat fl.oat ...

([] means optional.)

Each row is a node, each column a TR, separated by whitespace. This is the
.1D surface fmri data format created by NIH AFNI/SUMA's 3dVol2Surf, which
interpolates voxel rows data onto each node (vertex) in a MGH FreeSurfer
surface. But the rows could be any rows data; no spatial information is
used (nor is it supplied by any recognized input arguments).

Normalizes each row by default. Saves correlationss as platform single-precision
floats with .1D.gz output an option. Uses Cblas API so link with the Cblas
implementation library of choice for your platform (Atlas, Gotoblas, MKL,
vecLib, cuBLAS, OpenCL, etc). Tested with GotoBlas2 on 64-bit Linux and vecLib
on 64-bit Mac OS X

Most CBLAS implementations are multicore, so the parallelism (within-machine)
should come from that.

Ben Singer <bdsinger@princeton.edu> April 2013, revised October 2014

*/

#define _GNU_SOURCE
#include <libgen.h>
#include <assert.h>
#include <string.h>
#include <limits.h>
#include <stdlib.h>

#ifdef __APPLE__
#include <malloc/malloc.h>
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif

#include "pnicorr.h"

// ******* Memory limits
#define BYTES_PER_GB (1000000000LL)
#define BYTES_PER_MB (1000000LL)
#define MAX_FILENAME 256
#define CBLAS_ALPHA (1.0F)
#define CBLAS_BETA (0.0F)
#define PNICORR_ALIGNMENT 64
enum { PNICORR_CORRELATIONS_TAG = 1, PNICORR_NUMTAGS };

#define PNICORR_PROC_BODY_PARAMS                                               \
  (const float *restrict full_matrix, float *restrict submatrix,               \
   float *restrict correlations, const int *restrict offset_for_task,          \
   const int *restrict rows_for_task, const int nt, const int num_rows,        \
   const int num_columns, const char *restrict outfile,                        \
   const pnicorr_iotype_t outtype, const int my_task_id)
#define PNICORR_PROC_BODY_ARGS                                                 \
  (full_matrix, submatrix, correlations, offset_for_task, rows_for_task, nt,   \
   num_rows, num_columns, outfile, outtype, my_task_id)
#define PNICORR_PARAMSET_PARAMS                                                \
  (const int my_task_id, const int *num_rows, const int *num_columns,          \
   long long *ntasks, int *num_workers, float **restrict full_matrix)
#define PNICORR_PARAMSET_ARGS                                                  \
  (my_task_id, &num_rows, &num_columns, &ntasks, &num_workers, &full_matrix)
#define PNICORR_TASK_PROCESSING_PARAMS                                         \
  (const float *restrict full_matrix, float *restrict submatrix,               \
   float *restrict correlations, const int offset_for_task,                    \
   const int rows_for_task, const int num_rows, const int num_columns)
#define PNICORR_TASK_PROCESSING_ARGS(id)                                       \
  (full_matrix, submatrix, correlations, offset_for_task[id],                  \
   rows_for_task[id], num_rows, num_columns)

#define COR_TASK_PROCESSING(id)                                                \
  pnicorr_task_processing PNICORR_TASK_PROCESSING_ARGS(id)

#ifdef HAVE_MPI
/* use MPI */
#include <mpi.h>

void pnicorr_mpi_init(int *restrict my_task_id, int *restrict num_workers);
void pnicorr_mpi_finalize(void);
void pnicorr_mpi_paramset PNICORR_PARAMSET_PARAMS;
void pnicorr_mpi_proc_body PNICORR_PROC_BODY_PARAMS;

#define COR_MPI_INIT(my_task_id, num_workers)                                  \
  pnicorr_mpi_init(&my_task_id, &num_workers)

#define COR_MPI_FINALIZE pnicorr_mpi_finalize()

#define COR_MPI_PARAMSET pnicorr_mpi_paramset PNICORR_PARAMSET_ARGS

#define COR_TASK_PROC_BODY pnicorr_mpi_proc_body PNICORR_PROC_BODY_ARGS

#else
/* no MPI */

#define NOP (void)0
#define COR_MPI_INIT(my_task_id, num_workers) NOP
#define COR_MPI_FINALIZE NOP
#define COR_MPI_PARAMSET NOP
#define COR_TASK_PROC_BODY pnicorr_proc_body PNICORR_PROC_BODY_ARGS

#endif

// ******* Helper functions

// ****************** normalize_fvec ***************
static void normalize_fvec(float *restrict full_matrix, const int num_rows,
                           const int N) {
  const int incX = 1;
  float *dptr = full_matrix;
  for (int i = 0; i < num_rows; ++i) {
    float norm2 = cblas_snrm2(N, dptr, incX);
    const float alpha = 1.0 / norm2;
    cblas_sscal(N, alpha, dptr, incX);
    dptr += N;
  }
}

void pnicorr_task_processing PNICORR_TASK_PROCESSING_PARAMS;

void pnicorr_proc_body PNICORR_PROC_BODY_PARAMS;

void pnicorr_savematrix_for_task(const int i, const int *restrict rows_for_task,
                                 const int num_rows, const char *outfile,
                                 const pnicorr_iotype_t outtype,
                                 const float *restrict correlations);

// ****************** MAIN ***************
int main(const int argc, const char *argv[]) {
  float *full_matrix = NULL;
  float *correlations = NULL;
  char outfile[MAX_FILENAME];
  int num_rows = 0, num_columns = 0;
  int my_task_id = 0;
  long long ntasks = 1;
  int num_workers = 1;
  pnicorr_iotype_t outtype = pnicorr_iotype_numiotypes;

  COR_MPI_INIT(my_task_id, num_workers);

  struct opts2struct_t *ops2s = opts2struct_create();

  if (0 == my_task_id) {

    // ***** parse input args ******
    const char *ext;
    const char *comp;
    const int pre_args = 2;
    char *basec = strdup(argv[0]);
    char *bname = basename(basec);

    if (argc < pre_args) {
      fprintf(stderr, "\nusage: %s file.1D.dset[.gz] -[no]norm -mem=MB "
                      "-iotype=1D|1Dgz|mat\n",
              bname);
      exit(EXIT_FAILURE);
    }
    const char *filename = argv[1];

    opts2struct_parseopts(ops2s, argc, argv);

    int normalize = ops2s->i[norm];
    long long maxmem = ops2s->i[mem];

    if (!strncmp(ops2s->iotype, "1D", 2)) {
      outtype = pnicorr_iotype_1D;
      ext = "1D.dset";
      comp = "";
    } else if (!strncmp(ops2s->iotype, "1Dgz", 4)) {
      outtype = pnicorr_iotype_1Dgz;
      ext = "1D.dset";
      comp = ".gz";
    } else if (!strncmp(ops2s->iotype, "mat", 3)) {
      outtype = pnicorr_iotype_mat;
      ext = "mat";
      comp = "";
    } else {
      ext = "";
      comp = "";
      ERRIF(1, "unrecognized iotype %s", ops2s->iotype);
    }

    //  **********  get full_matrix from file ************
    full_matrix = pnicorr_load_1D(filename, &num_rows, &num_columns);

    // get outfile name
    pnicorr_genoutfile(bname, num_rows, ext, comp, outfile);
    free(basec);

    //  **********  Normalize -- unless asked not to via "-nonorm" ************
    if (normalize) {
      LOG("normalizing each row/rows independently ...\n");
      TIC;
      normalize_fvec(full_matrix, num_rows, num_columns);
      TOC("normalize");
    }

    // ************  do the matrix multiplication ***************
    // (results in correlation after normalization)
    long long corbytes =
        (long long)num_rows * (long long)num_rows * (long long)sizeof(float);
    // must split up the task s.t. each task uses MAX_MEM_PER_TASK bytes
    ntasks = (long long)((float)corbytes /
                             (float)((long long)maxmem * BYTES_PER_MB) +
                         0.5);
    LOG("%lld task(s), since %d^2 floats is %.3f MB and the max MB per "
        "task is %lld MB\n",
        ntasks, num_rows, (float)corbytes / (float)BYTES_PER_MB, maxmem);

    if (ntasks <= 1) {
      // get time info immediately
      setbuf(stdout, NULL);

      // *********   all in one go; enough memory to do so *********
      LOG("number of bytes needed for correlations is small enough for a "
          "single "
          "task\n");
      LOG("allocating memory...\n");
      int err =
          posix_memalign((void *)&correlations, PNICORR_ALIGNMENT, corbytes);
      ERRIF(err, "posix_memalign");
      LOG("calling cblas_ssyrk to perform AA' ...\n");

      TIC;
      cblas_ssyrk(CblasRowMajor, CblasUpper, CblasNoTrans, num_rows,
                  num_columns, CBLAS_ALPHA, full_matrix, num_columns,
                  CBLAS_BETA, correlations, num_rows);
      TOC("calculation");

      // *** save ***
      LOG("saving %d x %d task correlations to %s\n", num_rows, num_rows,
          outfile);
      TIC;
      pnicorr_savematrix(correlations, num_rows, num_rows, "w", outtype,
                         outfile);
      TOC("saving");
      LOG("done\n");
    } // ntasks > 1
  }   // my_task_id == 0

  COR_MPI_PARAMSET;

  if (ntasks > 1) {
    float *submatrix = NULL;

    // ****************  need to do it in blocks and append correlations
    // *********
    int nt = (int)ntasks;
    int num_rows_per_task = num_rows / nt;
    if (num_rows_per_task < 1) {
      num_rows_per_task = 1;
      nt = num_rows;
    }

    int *rows_for_task = calloc(num_workers, sizeof(int));
    int *offset_for_task = calloc(num_workers, sizeof(int));
    // have task 0 (master in case of MPI) do the work for leftover rows
    // when total doesn't evenly divide ntasks
    rows_for_task[0] = num_rows - (num_rows_per_task * (nt - 1));
    offset_for_task[0] = 0;
    int maxrows_for_tasks = (rows_for_task[0] > num_rows_per_task)
                                ? rows_for_task[0]
                                : num_rows_per_task;
    for (int i = 1; i < nt; ++i) {
      offset_for_task[i] = offset_for_task[i - 1] + rows_for_task[i - 1];
      rows_for_task[i] = num_rows_per_task;
    }
    int err = posix_memalign((void *)&submatrix, PNICORR_ALIGNMENT,
                             maxrows_for_tasks * num_columns * sizeof(float));
    ERRIF(err, "posix_memalign submatrix");
    err = posix_memalign((void *)&correlations, PNICORR_ALIGNMENT,
                         maxrows_for_tasks * num_rows * sizeof(float));
    ERRIF(err, "posix_memalign correlations");

    COR_TASK_PROC_BODY;

    free(rows_for_task);
    free(offset_for_task);
    if (submatrix) {
      free(submatrix);
      submatrix = NULL;
    }
  }

  if (full_matrix) {
    free(full_matrix);
    full_matrix = NULL;
  }
  if (correlations) {
    free(correlations);
    correlations = NULL;
  }

  COR_MPI_FINALIZE;

  return EXIT_SUCCESS;
}

void pnicorr_task_processing PNICORR_TASK_PROCESSING_PARAMS {

  LOG("submatrix %p full_matrix %p\n", submatrix, full_matrix);

  memcpy(submatrix, full_matrix + offset_for_task,
         rows_for_task * num_columns * sizeof(float));
  memset(correlations, 0, num_rows * rows_for_task * sizeof(float));

  TIC;
  LOG("calling cblas_sgemm using %d num_rows (rows) at a time\n",
      rows_for_task);
  cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, rows_for_task, num_rows,
              num_columns, CBLAS_ALPHA, submatrix, num_columns, full_matrix,
              num_columns, CBLAS_BETA, correlations, num_rows);
  TOC("calculating");
}

void pnicorr_savematrix_for_task(const int i, const int *restrict rows_for_task,
                                 const int num_rows, const char *outfile,
                                 const pnicorr_iotype_t outtype,
                                 const float *restrict correlations) {
  // ******** save/append *********
  LOG("task %d: %s %d x %d task correlations to %s\n", i,
      (i == 0) ? "writing" : "appending", rows_for_task[i], num_rows, outfile);
  TIC;
  pnicorr_savematrix(correlations, rows_for_task[i], num_rows,
                     (i == 0) ? "w" : "a", outtype, outfile);
  TOC("saving");
}

#ifdef HAVE_MPI

void pnicorr_mpi_init(int *my_task_id, int *num_workers) {
  int initialized;
  MPI_Initialized(&initialized);
  if (0 == initialized) {
    MPI_Init(NULL, NULL);
  }
  MPI_Comm_rank(MPI_COMM_WORLD, my_task_id);
  MPI_Comm_size(MPI_COMM_WORLD, num_workers);
}

void pnicorr_mpi_finalize(void) {
  int initialized;
  MPI_Initialized(&initialized);
  if (initialized) {
    int finalized;
    MPI_Finalized(&finalized);
    if (0 == finalized) {
      MPI_Finalize();
    }
  }
}

void pnicorr_mpi_paramset PNICORR_PARAMSET_PARAMS {

  MPI_Bcast((void *)num_rows, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast((void *)num_columns, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast((void *)ntasks, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
  MPI_Bcast((void *)num_workers, 1, MPI_INT, 0, MPI_COMM_WORLD);

  long long fullmatrixbytes = *num_rows * *num_columns * sizeof(float);
  if (my_task_id > 0) {
    int err =
        posix_memalign((void *)full_matrix, PNICORR_ALIGNMENT, fullmatrixbytes);
    ERRIF(err, "posix_memalign");
  } else {
    ALOGIF(*num_workers != *ntasks, "setting ntasks equal to num_workers: (%d "
                                    "workers, %d tasks -> %d tasks)\n",
           *num_workers, (int)*ntasks, *num_workers);
  }
  *ntasks = *num_workers;

  MPI_Bcast((void *)(*full_matrix), *num_rows * *num_columns, MPI_FLOAT, 0,
            MPI_COMM_WORLD);

  MPI_Barrier(MPI_COMM_WORLD);
}

void pnicorr_mpi_proc_body PNICORR_PROC_BODY_PARAMS {
  COR_TASK_PROCESSING(my_task_id);

  if (0 == my_task_id) {
    pnicorr_savematrix_for_task(0, rows_for_task, num_rows, outfile, outtype,
                                correlations);
    for (int i = 1; i < nt; ++i) {
      MPI_Recv(correlations, num_rows * rows_for_task[i], MPI_FLOAT, i,
               PNICORR_CORRELATIONS_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      pnicorr_savematrix_for_task(i, rows_for_task, num_rows, outfile, outtype,
                                  correlations);
    }
  } else {
    LOG("[%d] sending %d x %d correlations\n", my_task_id,
        rows_for_task[my_task_id], num_rows);
    MPI_Send(correlations, num_rows * rows_for_task[my_task_id], MPI_FLOAT, 0,
             PNICORR_CORRELATIONS_TAG, MPI_COMM_WORLD);
  }
}

#else

void pnicorr_proc_body PNICORR_PROC_BODY_PARAMS {
  int num_rows_processed = 0;
  for (int i = 0; i < nt; ++i) {

    COR_TASK_PROCESSING(i);

    pnicorr_savematrix_for_task(i, rows_for_task, num_rows, outfile, outtype,
                                correlations);
    num_rows_processed += rows_for_task[i];

    ALOG("%d%% complete\n",
         (int)((float)num_rows_processed / (float)num_rows * 100.0f));
  }
}

#endif

max_stars_repo_name: bdsinger/pnicorr

max_issues_repo_name: bdsinger/pnicorr

max_forks_repo_name: bdsinger/pnicorr


text:
#ifndef EMPHASIS_RINIT_H_INCLUDED
#define EMPHASIS_RINIT_H_INCLUDED

#include <nlopt.h>

#if (defined(_WIN32) || defined(__WIN32__)) && !defined(__LCC__)
#  define REMP_EXPORT extern "C" __declspec(dllexport)
#else
#  define REMP_EXPORT extern "C"
#endif

#ifdef __cplusplus
extern "C" {
#endif

/* some function pointers into nlopt to be filled in by R */

REMP_EXPORT nlopt_opt(*remp_create)(nlopt_algorithm, unsigned);
REMP_EXPORT void(*remp_destroy)(nlopt_opt);
REMP_EXPORT nlopt_result(*remp_optimize)(nlopt_opt, double *, double *);
REMP_EXPORT nlopt_result(*remp_set_min_objective)(nlopt_opt, nlopt_func, void *);
REMP_EXPORT nlopt_result(*remp_set_max_objective)(nlopt_opt, nlopt_func, void *);
REMP_EXPORT nlopt_result(*remp_set_lower_bounds)(nlopt_opt, const double *);
REMP_EXPORT nlopt_result(*remp_set_lower_bounds1)(nlopt_opt, double);
REMP_EXPORT nlopt_result(*remp_set_upper_bounds)(nlopt_opt, const double *);
REMP_EXPORT nlopt_result(*remp_set_upper_bounds1)(nlopt_opt, double);
REMP_EXPORT nlopt_result(*remp_set_xtol_rel)(nlopt_opt, double);
REMP_EXPORT nlopt_result(*remp_set_xtol_abs)(nlopt_opt, double);

#ifdef __cplusplus
}
#endif

#endif

ext: h

max_stars_repo_name: franciscorichter/emphasis

max_issues_repo_name: franciscorichter/emphasis

max_forks_repo_name: franciscorichter/emphasis


text:
#ifndef DISTANCE_H_MUJPGVJL
#define DISTANCE_H_MUJPGVJL

#include <algorithm>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/median.hpp>
#include <boost/accumulators/statistics/skewness.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/variance.hpp>
#include <boost/histogram.hpp>
#include <cstddef>
#include <gsl/gsl>
#include <opencv2/core/persistence.hpp>

namespace sens_loc {

/// This namespace contains useful function for statistical analysis of
/// keypoints, descriptors and other relevant quantities.
namespace analysis {

/// Wrapper-struct to hold the common set of statistical values.
struct statistic {
    using accumulator_t = boost::accumulators::accumulator_set<
        float,
        boost::accumulators::stats<boost::accumulators::tag::count,
                                   boost::accumulators::tag::min,
                                   boost::accumulators::tag::max,
                                   boost::accumulators::tag::median,
                                   boost::accumulators::tag::mean,
                                   boost::accumulators::tag::variance(
                                       boost::accumulators::lazy),
                                   boost::accumulators::tag::skewness>>;

    std::size_t        count    = 0UL;
    float              min      = 0.0F;
    float              max      = 0.0F;
    float              mean     = 0.0F;
    float              variance = 0.0F;
    float              stddev   = 0.0F;
    float              skewness = 0.0F;
    float              median   = 0.0F;
    std::vector<float> decentils{0.0F, 0.0F, 0.0F, 0.0F, 0.0F,
                                 0.0F, 0.0F, 0.0F, 0.0F};

    /// Extract the statistical values from the original data.
    /// \pre the data must be sorted!
    static statistic make(gsl::span<const float> data) {
        using namespace std;
        Expects(is_sorted(begin(data), end(data)));

        statistic s;

        if (data.empty())
            return s;

        if (data.size() >= 10L) {
            const int      number_quantils = 10;
            std::ptrdiff_t delta           = data.size() / number_quantils;

            Expects(s.decentils.size() == 9UL);
            for (int i = 1; i < number_quantils; ++i) {
                s.decentils[i - 1] = data[i * delta];
            }
            Ensures(s.decentils.size() == 9UL);
        }
        s.median = data[data.size() / 2];

        accumulator_t accu;
        std::for_each(std::begin(data), std::end(data),
                      [&accu](float e) { accu(e); });
        using namespace boost::accumulators;
        s.count    = boost::accumulators::count(accu);
        s.min      = boost::accumulators::min(accu);
        s.max      = boost::accumulators::max(accu);
        s.mean     = boost::accumulators::mean(accu);
        s.variance = boost::accumulators::variance(accu);
        s.stddev   = std::sqrt(s.variance);
        s.skewness = boost::accumulators::skewness(accu);
        return s;
    }

    void reset() noexcept {
        count     = 0UL;
        min       = 0.0F;
        max       = 0.0F;
        mean      = 0.0F;
        variance  = 0.0F;
        stddev    = 0.0F;
        skewness  = 0.0F;
        median    = 0.0F;
        decentils = {0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F};
    }
};

/// Write-Functionality for OpenCVs-Filestorage API.
void write(cv::FileStorage& fs, const std::string& name, const statistic& stat);

/// This class processes 'float' datapoints and calculates both histograms
/// and basic statistical quantities, like the 'min', 'max', 'mean' and
/// others.
/// Its main use-case is to analyze distance information gathered while
/// processing image and matching data. It can be used for any scalar
/// quantity though.
class distance {
  public:
    using axis_t = boost::histogram::axis::regular<
        // 'float' values are tracked with the histogram.
        float,
        // do no transformation before insertion
        boost::histogram::axis::transform::id,
        // a string is the metadata for the axis (=title)
        std::string,
        // no overflow/underflow by construction.
        boost::histogram::axis::option::none_t>;
    using histo_t = decltype(boost::histogram::make_histogram(axis_t{}));

    distance() = default;
    explicit distance(gsl::span<const float> distances,
                      unsigned int           bin_count = 50U,
                      std::string            title     = "distance") noexcept
        : _bin_count{bin_count}
        , _axis_title{std::move(title)} {
        analyze(distances);
    }

    /// Analyses the provided distances. This will overwrite a previous
    /// analysis and dataset! Use this method to provide the data if the
    /// class was default constructed.
    /// \pre is_sorted(distances)
    void analyze(gsl::span<const float> distances, bool histo = true) noexcept;

    /// Return the reference to the potentially created histogram in \c analyze.
    /// If the histogram is not created, it will just be default constructed.
    /// \sa analyze
    [[nodiscard]] const histo_t& histogram() const noexcept { return _histo; }

    [[nodiscard]] std::size_t count() const noexcept { return _s.count; }
    [[nodiscard]] float       min() const noexcept { return _s.min; }
    [[nodiscard]] float       max() const noexcept { return _s.max; }
    [[nodiscard]] float       median() const noexcept { return _s.median; }
    [[nodiscard]] float       mean() const noexcept { return _s.mean; }
    [[nodiscard]] float       variance() const noexcept { return _s.variance; }
    [[nodiscard]] float       stddev() const noexcept { return _s.stddev; }
    [[nodiscard]] float       skewness() const noexcept { return _s.skewness; }

    /// Read-Only access to underlying statistical data.
    [[nodiscard]] const statistic& get_statistic() const noexcept { return _s; }

  private:
    unsigned int _bin_count  = 50UL;
    std::string  _axis_title = "distance";
    histo_t      _histo;
    statistic    _s;
};
}  // namespace analysis
}  // namespace sens_loc

#endif /* end of include guard: DISTANCE_H_MUJPGVJL */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <signal.h>
#include <gsl/gsl_rng.h>

#include "allvars.h"
#include "proto.h"

#ifdef DEBUG
#include <fenv.h>
void enable_core_dumps_and_fpu_exceptions(void)
{
  struct rlimit rlim;
  extern int feenableexcept(int __excepts);

  /* enable floating point exceptions */

  /*
     feenableexcept(FE_DIVBYZERO | FE_INVALID);
   */

  /* Note: FPU exceptions appear not to work properly
   * when the Intel C-Compiler for Linux is used
   */

  /* set core-dump size to infinity */
  getrlimit(RLIMIT_CORE, &rlim);
  rlim.rlim_cur = RLIM_INFINITY;
  setrlimit(RLIMIT_CORE, &rlim);

  /* MPICH catches the signales SIGSEGV, SIGBUS, and SIGFPE....
   * The following statements reset things to the default handlers,
   * which will generate a core file.
   */
  /*
     signal(SIGSEGV, catch_fatal);
     signal(SIGBUS, catch_fatal);
     signal(SIGFPE, catch_fatal);
     signal(SIGINT, catch_fatal);
   */

  signal(SIGSEGV, SIG_DFL);
  signal(SIGBUS, SIG_DFL);
  signal(SIGFPE, SIG_DFL);
  signal(SIGINT, SIG_DFL);

  /* Establish a handler for SIGABRT signals. */
  signal(SIGABRT, catch_abort);
}

void catch_abort(int sig)
{
  MPI_Finalize();
  exit(0);
}

void catch_fatal(int sig)
{
  terminate_processes();
  MPI_Finalize();

  signal(sig, SIG_DFL);
  raise(sig);
}

void terminate_processes(void)
{
  pid_t my_pid;
  char buf[500], hostname[500], *cp;
  char commandbuf[500];
  FILE *fd;
  int i, pid;

  sprintf(buf, "%s%s", All.OutputDir, "PIDs.txt");

  my_pid = getpid();

  if((fd = fopen(buf, "r")))
    {
      for(i = 0; i < NTask; i++)
	{
	  int ret;

	  ret = fscanf(fd, "%s %d", hostname, &pid);

	  cp = hostname;
	  while(*cp)
	    {
	      if(*cp == '.')
		*cp = 0;
	      else
		cp++;
	    }

	  if(my_pid != pid)
	    {
	      sprintf(commandbuf, "ssh %s kill -ABRT %d", hostname, pid);
	      printf("--> %s\n", commandbuf);
	      fflush(stdout);
#ifndef NOCALLSOFSYSTEM
	      ret = system(commandbuf);
#endif
	    }
	}

      fclose(fd);
    }
}

void write_pid_file(void)
{
  pid_t my_pid;
  char mode[8], buf[500];
  FILE *fd;
  int i;

  my_pid = getpid();

  sprintf(buf, "%s%s", All.OutputDir, "PIDs.txt");

  if(RestartFlag == 0)
    strcpy(mode, "w");
  else
    strcpy(mode, "a");

  for(i = 0; i < NTask; i++)
    {
      if(ThisTask == i)
	{
	  if(ThisTask == 0)
	    sprintf(mode, "w");
	  else
	    sprintf(mode, "a");

	  if((fd = fopen(buf, mode)))
	    {
	      fprintf(fd, "%s %d\n", getenv("HOST"), (int) my_pid);
	      fclose(fd);
	    }
	}

      MPI_Barrier(MPI_COMM_WORLD);
    }
}
#endif

double get_random_number(unsigned int id)
{
  return RndTable[(id % RNDTABLE)];
}

void set_random_numbers(void)
{
  int i;

  for(i = 0; i < RNDTABLE; i++)
    RndTable[i] = gsl_rng_uniform(random_generator);
}

/* returns the number of cpu-ticks in seconds that
 * have elapsed. (or the wall-clock time)
 */
double second(void)
{
#ifdef WALLCLOCK
  return MPI_Wtime();
#else
  return ((double) clock()) / CLOCKS_PER_SEC;
#endif

  /* note: on AIX and presumably many other 32bit systems,
   * clock() has only a resolution of 10ms=0.01sec
   */
}

double measure_time(void)	/* strategy: call this at end of functions to account for time in this function, and before another (nontrivial) function is called */
{
  double t, dt;

  t = second();
  dt = t - WallclockTime;
  WallclockTime = t;

  return dt;
}

/* returns the time difference between two measurements
 * obtained with second(). The routine takes care of the
 * possible overflow of the tick counter on 32bit systems.
 */
double timediff(double t0, double t1)
{
  double dt;

  dt = t1 - t0;

  if(dt < 0)			/* overflow has occured (for systems with 32bit tick counter) */
    {
#ifdef WALLCLOCK
      dt = 0;
#else
      dt = t1 + pow(2, 32) / CLOCKS_PER_SEC - t0;
#endif
    }

  return dt;
}

#ifdef X86FIX

#define _FPU_SETCW(x) asm volatile ("fldcw %0": :"m" (x));
#define _FPU_GETCW(x) asm volatile ("fnstcw %0":"=m" (x));
#define _FPU_EXTENDED 0x0300
#define _FPU_DOUBLE   0x0200

void x86_fix(void)
{
  unsigned short dummy, new_cw;
  unsigned short *old_cw;

  old_cw = &dummy;

  _FPU_GETCW(*old_cw);
  new_cw = (*old_cw & ~_FPU_EXTENDED) | _FPU_DOUBLE;
  _FPU_SETCW(new_cw);
}

#endif

void sumup_large_ints(int n, int *src, long long *res)
{
  int i, j, *numlist;

  numlist = (int *) mymalloc(NTask * n * sizeof(int));
  MPI_Allgather(src, n, MPI_INT, numlist, n, MPI_INT, MPI_COMM_WORLD);

  for(j = 0; j < n; j++)
    res[j] = 0;

  for(i = 0; i < NTask; i++)
    for(j = 0; j < n; j++)
      res[j] += numlist[i * n + j];

  myfree(numlist);
}

void sumup_longs(int n, long long *src, long long *res)
{
  int i, j;
  long long *numlist;

  numlist = (long long *) mymalloc(NTask * n * sizeof(long long));
  MPI_Allgather(src, n * sizeof(long long), MPI_BYTE, numlist, n * sizeof(long long), MPI_BYTE,
		MPI_COMM_WORLD);

  for(j = 0; j < n; j++)
    res[j] = 0;

  for(i = 0; i < NTask; i++)
    for(j = 0; j < n; j++)
      res[j] += numlist[i * n + j];

  myfree(numlist);
}

size_t sizemax(size_t a, size_t b)
{
  if(a < b)
    return b;
  else
    return a;
}

void report_VmRSS(void)
{
  pid_t my_pid;
  FILE *fd;
  char buf[1024];

  my_pid = getpid();

  sprintf(buf, "/proc/%d/status", my_pid);

  if((fd = fopen(buf, "r")))
    {
      while(1)
	{
	  if(fgets(buf, 500, fd) != buf)
	    break;

	  if(strncmp(buf, "VmRSS", 5) == 0)
	    {
	      printf("ThisTask=%d: %s", ThisTask, buf);
	    }
	  if(strncmp(buf, "VmSize", 6) == 0)
	    {
	      printf("ThisTask=%d: %s", ThisTask, buf);
	    }
	}
      fclose(fd);
    }
}

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp