text:
/* gcPresAbs.c - Calculate probability of expression in an mrna
   sample given a cel file of results. Use empiracal rank within
   a gc population to estimate probability and fisher's method
   of combining probabilities to calculate probe set probability.
*/
#include "common.h"
#include "hash.h"
#include "linefile.h"
#include "obscure.h"
#include "options.h"
#include "dystring.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>

struct row1lq
/* A single row from a affymetrix 1lq file. Also used for
   storing values and probability of expression. */
{
    struct row1lq *next;    /* Next in list. */
    char *psName;           /* Name of probe set.*/
    char *seq;              /* Sequence of probe, in weird affy format. */
    int x;                  /* X coordinate from affy cel file. */
    int y;                  /* Y coordinate from affy cel file. */
    int gcCount;            /* Number of g's and c's in sequence. */
    int repCount;           /* Number of experiments we're looking at. */
    double *intensity;      /* Intensities for each experiment. */
    double *pVal;           /* Pvals for each experiment. */
};

struct psProb
/* Probability that a probe set is expressed. */
{
    struct psProb *next; /* Next in list. */
    char psName[128];    /* Name of probe set. */
    double pVal;         /* P-Value of expression. */
};

struct replicateMap
/* Map the replicates to eachother. */
{
    struct replicateMap *next;  /* Next in list. */
    char *rootName;             /* Root name of replicates. */
    int repCount;               /* Number of replicates. */
    int *repIndexes;            /* Indexes into matrix columns. */
};

struct row1lq **gcBins = NULL;     /* Bins of lists of row1lqs for different amounts of GCs. */
struct row1lq *badBin = NULL;
struct hash *probeSetHash = NULL;  /* Hash of slRefs for each probe set where slRef->val is a row1lq. */
struct hash *killHash = NULL;      /* Hash of probes to be ignored. */
struct hash *outputValsHash= NULL; /* Hash of probe sets to output values for. */
FILE *outputPsSetSeqs = NULL;      /* Output probe sets and sequences used here. */
FILE *outputIntenVals = NULL;       /* File to output intensity for individual values to. */
FILE *outputPvals = NULL;          /* File to output pvals for individual to. */
int minGcBin = 0;                  /* Minimum GC bin. */
int maxGcBin = 0;                  /* Maximum GC bin. */
struct psProb *probePVals = NULL;  /* List of psProbe structs with pVals filled in. */
int minProbeNum = 3;               /* Ignore probe sets that have less than this many probes. */
char *placeHolder = "dummy";       /* Hold a place in a hash. */
boolean combineCols = FALSE;       /* Combine columns that are the same samples? */
int repToSort = 0;                 /* What replicate to sort the gc bins by. */
int maxReplicates = 0;             /* What are the maximum number of replicates. */
int sampleCount = 0;               /* Number of samples. */
int numBadProbes = 0;              /* Number of bad probes counted. */
FILE *gCountsFile = NULL;              /* File to output 'G' counts per probe to. */
FILE *gcCountFile = NULL;             /* File to output 'G' + 'C' counts per probe to. */
FILE *cCountsFile = NULL;              /* File to output 'C' counts per probe. */
boolean doMedianComb = FALSE;          /* Do a median rather than Fisher's combination. */
boolean doCombo = FALSE; /* Do a combination of median and fisher. */

static struct optionSpec optionSpecs[] =
/* Our acceptable options to be called with. */
{
    {"help", OPTION_BOOLEAN},
    {"1lqFile", OPTION_STRING},
    {"outFile", OPTION_STRING},
    {"badProbes", OPTION_STRING},
    {"combineReplicates", OPTION_BOOLEAN},
    {"outputDuplicates", OPTION_BOOLEAN},
    {"testVals", OPTION_BOOLEAN},
    {"countPrefix", OPTION_STRING},
    {"doMedianComb", OPTION_BOOLEAN},
    {"testMedian", OPTION_BOOLEAN},
    {"doCombo", OPTION_BOOLEAN},
    {"outputVals", OPTION_STRING},
    {"outputValsList", OPTION_STRING},
    {"outputPsSetSeqs", OPTION_STRING},
    {NULL,0}
};

static char *optionDescripts[] =
/* Description of our options for usage summary. */
{
    "Display this message.",
    "1lq file from affymetrix specifying probe layout on chip",
    "File to output matrix of probabilities to.",
    "File with sequence of probes that are known to be bad and should be ignored.",
    "Combine experiments with same name prefix.",
    "When used with combineReplicates, outputs a column for each tissue that has same prefix.",
    "Skip everything else and just run inverse chisq on each value degreeFreedom pair on command line.",
    "Output G/C probe countent counts to prefixCounts.tab",
    "Use median rather than Fisher's combination",
    "Test the median code.",
    "Do a combination of medians and fischer's combination.",
    "Output the intensity and probability to file with this prefix.",
    "Only output probe sets in this list.",
    "Output probe sets and sequences in 5'->3' orientation"
};

void usage()
/** Print usage and quit. */
{
int i=0;
warn("gcPresAbs - Calculate probability of expression in an mrna\n"
     "sample given a cel file of results. Use empiracal rank within\n"
     "a gc population to estimate probability and fisher's method\n"
     "of combining probabilities to calculate probe set probability.");
for(i=0; i<ArraySize(optionSpecs) -1; i++)
    fprintf(stderr, "  -%s -- %s\n", optionSpecs[i].name, optionDescripts[i]);
errAbort("\nusage:\n   "
	 "gcPresAbs -1lqFile=chipmousea.1lq -outFile=probeSet.matrix.probs file1.cel file2.cel ... fileN.cel");
}

int psCmp(const void *va, const void *vb)
/* Compare based on probeSet name. */
{
const struct psProb *a = *((struct psProb **)va);
const struct psProb *b = *((struct psProb **)vb);
return strcmp(a->psName, b->psName);
}

int row1lqCmp(const void *va, const void *vb)
/* Compare to sort based on intensity for the
   given repliate. */
{
const struct row1lq *a = *((struct row1lq **)va);
const struct row1lq *b = *((struct row1lq **)vb);
double bI = b->intensity[repToSort];
double aI = a->intensity[repToSort];
if(aI < bI)
    return 1;
if(aI > bI)
    return -1;
return 0;
}

int countGc(char *sequence)
/* Count the numbers of G's+C's in the sequence. */
{
int gCount = 0, cCount =0;
tolowers(sequence);
gCount = countChars(sequence, 'g');
if(gCountsFile)
    fprintf(gCountsFile, "%d\n", gCount);
cCount = countChars(sequence, 'c');
if(cCountsFile)
    fprintf(cCountsFile, "%d\n", cCount);
if(gcCountFile)
    fprintf(gcCountFile, "%d\n", cCount + gCount);
return cCount + gCount;
}

void initGcBins(int minGc, int maxGc)
/* Intialize the gcBins structure for maxGc-minGc bins. */
{
int numBins = maxGc - minGc + 1;
minGcBin = minGc;
maxGcBin = maxGc;
AllocArray(gcBins, numBins);
}

int binForGc(int gcCount)
/* Return the correct bin for this gc count. */
{
int count = 0;
if(gcCount < minGcBin)
    gcCount = minGcBin;
if(gcCount > maxGcBin)
    gcCount = maxGcBin;
return gcCount - minGcBin;
}

void readBadProbes(char *fileName)
/* Read in a list that contains bad probes sequence should be in first row. */
{
struct lineFile *lf = NULL;
char *line = NULL;
killHash = newHash(15);
lf = lineFileOpen(fileName, TRUE);
while(lineFileNextReal(lf, &line))
    {
    tolowers(line);
    hashAdd(killHash, line, placeHolder);
    }
lineFileClose(&lf);
}

void read1lqFile(char *name, struct row1lq ***pRowArray, int *pRowCount)
/* Read in the 1lq popluating the gcBins and hasing by probe set name.
   Need to index rows by a few different ways:
   1) Index in file, for filling in intensity vals from cel files.
   2) Gc Bin for calculating the rank within a set of probes that have same gc content.
   3) Probe set for calculating the probability that the probe set as a whole is expressed.
*/
{
struct lineFile *lf = lineFileOpen(name, TRUE);
struct row1lq **array1lq = NULL;
struct row1lq *rowList = NULL, *row=NULL, *rowNext = NULL;
int i;
char *words[50];
char *line = NULL;
int count =0;
boolean found = FALSE;
struct slRef *refList = NULL, *ref = NULL;
while(lineFileNextReal(lf, &line))
    {
    if(stringIn("X\tY\t", line))
	{
	found = TRUE;
	break;
	}
    }
if(!found)
    errAbort("Couldn't find \"X\tY\t\" in 1lq file: %s", name);

while(lineFileChopNextTab(lf, words, sizeof(words)))
    {
    AllocVar(row);
    row->psName = cloneString(words[5]);
    row->seq = cloneString(words[2]);
    row->x = atoi(words[0]);
    row->y = atoi(words[1]);
    slAddHead(&rowList, row);
    (*pRowCount)++;
    }
slReverse(&rowList);
AllocArray(array1lq, (*pRowCount));

probeSetHash = newHash(15);

for(row = rowList; row != NULL; row = rowNext)
    {
    int bin = -1;
    struct slRef *ref = NULL, *reflist = NULL;

    rowNext = row->next;

    /* Keep track of the row by its index in the file. */
    array1lq[count++] = row;

    /* Skip things that are on the kill list. */
    tolowers(row->seq);
    if(sameString(row->seq, "!") || (killHash != NULL && hashFindVal(killHash, row->seq) != NULL))
	{
	slAddHead(&badBin, row);
	numBadProbes++;
	continue;
	}

    /* Add it to the correct gc bin. */
    row->gcCount = countGc(row->seq);
    bin = binForGc(row->gcCount);
    slAddHead(&gcBins[bin], row);

    /* Create a list of references indexed by the probe sets. */
    AllocVar(ref);
    ref->val = row;
    refList = hashFindVal(probeSetHash, row->psName);
    if(refList == NULL)
	hashAddUnique(probeSetHash, row->psName, ref);
    else
	slAddTail(&refList, ref);
    }

for(i = minGcBin; i < maxGcBin; i++)
    slReverse(&gcBins[i-minGcBin]);

//lineFileClose(&lf);
*pRowArray = array1lq;
}

void fillInMatrix(double **matrix, int celIx, int rowCount, char *celFile)
/* Read in the values of the cel file into matrix column celIx. */
{
struct lineFile *lf = lineFileOpen(celFile, TRUE);
char *words[50];
char *line = NULL;
boolean found = FALSE;
int count =0;
/* Read out the header. */
while(lineFileNextReal(lf, &line))
    {
    if(stringIn("[INTENSITY]", line))
	{
	found = TRUE;
	break;
	}
    }
if(!found)
    errAbort("Couldn't find \"[INTENSITY]\" in cel file: %s", celFile);
count =0;
lineFileNextReal(lf, &line); /* NumberCells... */
lineFileNextReal(lf, &line); /* CellHeader...*/
/* Read in the data. */
while(lineFileChopNext(lf, words, sizeof(words)))
    {
    if(*(words[0]) == '[')
	break;
    assert(count+1 <= rowCount);
    matrix[count][celIx] = atof(words[2]);
    count++;
    if(lf->bytesInBuf == 0)
	warn("Something weird is happening...");
    }
//warn("Closing line file.");
lineFileClose(&lf);
}

int xy2i(int x, int y)
/* Take in the x,y coordinates on the array and
   return the index into the array. */
{
return (712 * y) + x;
}

int doubleComp(const void *a, const void *b)
/* Return difference to sort on doubles. */
{
double x = (*(double *)a);
double y = (*(double *)b);
if(x > y)
    return 1;
else if(y > x)
    return -1;
else
    return 0;
}

/* double quantileArray(double *array, int count, double quantile) */
/* /\* give the quantile at a given position of the array. *\/ */
/* { */
/* double *results = NULL; */
/* boolean isEven = FALSE; */
/* int i = 0; */
/* double quant = 0; */

/* assert(count); */
/* assert(quantile >=0 && quantile <=1); */
/* AllocArray(results, count); */
/* if(floor(count * quantile) == ceil(count *quantile)) */
/*     isEven = TRUE; */

/* if(isEven && (count % 2 != 0)) */
/*     errAbort("Wrong even for quantile %.2f with count %d", quantile, count); */
/* if(!isEven && (count % 2 == 0)) */
/*     errAbort("Wrong odd for quantile %.2f with count %d", quantile, count); */

/* for(i = 0; i < count; i++) */
/*     results[i] = array[i]; */
/* qsort(results, count, sizeof(results[0]), doubleComp); */
/* if(isEven) */
/*     quant = (results[(int)floor(count * quantile)-1] + results[(int)floor(count*quantile)])/2; */
/* else */
/*     quant = results[(int)floor(count * quantile)]; */
/* freez(&results); */
/* return quant; */
/* } */

/* double medianArray(double *array, int count)  */
/* /\* Give the median, or a quantile of an array. *\/ */
/* { */
/* return quantileArray(array, count, quantile); */
/* } */

double medianArray(double *array, int count)
/* give the median of the array. */
{
double *results = NULL;
boolean isEven = FALSE;
int i = 0;
double median = 0;

assert(count);
AllocArray(results, count);
if(count % 2 == 0)
    isEven = TRUE;

for(i = 0; i < count; i++)
    results[i] = array[i];
qsort(results, count, sizeof(results[0]), doubleComp);
if(isEven)
    median = (results[count/2-1] + results[count/2])/2;
else
    median = results[count/2];
freez(&results);
return median;
}

void printDoubleArray(double *array, int count)
/* Print the elements of a double array. */
{
int i = 0;
for(i = 0; i<count; i++)
    printf("%.3f, ", array[i]);
printf("\n");
}

double medianRefList(struct slRef *refList)
/** Return the median value of the row list. */
{
struct slRef *ref = NULL;
struct row1lq *row = NULL;
int count = 0, currentIx = 0;
boolean isEven = FALSE;
double median = 0;
double *results = NULL;
int repIx = 0;

/* How much memory do we need? */
for(ref = refList; ref != NULL; ref = ref->next)
    {
    row = ref->val;
    count += row->repCount;
    }

assert(count > 0);
AllocArray(results, count);
/* Fill in the values. */
for(ref = refList; ref != NULL; ref = ref->next)
    {
    row = ref->val;
    for(repIx = 0; repIx < row->repCount; repIx++)
	{
	results[currentIx++] = row->pVal[repIx];
	}
    }

/* Get the median. */
median = medianArray(results, count);

freez(&results);
return 1 - median;
}

void probeSetCalcPval(void *val)
/* Loop through the slRef list calculating pVal based
   on row1lq->pval using fishers method for combining probabilities.
   Create a psProbe val and add it to the probePVals list. */
{
struct slRef *ref = NULL, *refList = NULL;
struct psProb *ps = NULL;
struct row1lq *row = NULL;
boolean debug = FALSE;
double probProduct = 0;
int probeCount = 0;
int probRepCount = 0;
int i = 0, j = 0;
refList = val;
probeCount = slCount(refList);
if(outputPsSetSeqs != NULL)
    {
    assert(refList);
    row = refList->val;
    fprintf(outputPsSetSeqs, "%s\t", row->psName);
    for(ref = refList; ref != NULL; ref = ref->next)
	{
	row = ref->val;
	reverseBytes(row->seq, strlen(row->seq));
	if(ref->next != NULL)
	    fprintf(outputPsSetSeqs, "%s,", row->seq);
	else
	    fprintf(outputPsSetSeqs, "%s", row->seq);

	reverseBytes(row->seq, strlen(row->seq));
	}
    fprintf(outputPsSetSeqs, "\n");
    }
/* Allocate some memory. */
AllocVar(ps);
row = refList->val;
ps->pVal = 0;
safef(ps->psName, sizeof(ps->psName), "%s", row->psName);

if(probeCount < minProbeNum)
    {
    ps->pVal=0;
    }
else if(doMedianComb)
    {
    ps->pVal = medianRefList(refList);
    }
else
    {
    for(ref = refList; ref != NULL; ref = ref->next)
	{
	row = ref->val;
	if(doCombo)
	    {
	    double med = medianArray(row->pVal, row->repCount);
	    if(med == 0)
		med = .00000000001;
	    probRepCount++;
	    probProduct += log(med);
	    }
	else
	    {
	    for(i = 0; i < row->repCount; i++)
		{
		probRepCount++;
		if(row->pVal[i] != 0)
		    probProduct = probProduct + log(row->pVal[i]);
		else
		    probProduct = probProduct + log(.000000001);
		}
	    }
	}
/* Fisher's method for combining probabilities. */
    ps->pVal = gsl_cdf_chisq_P(-2*probProduct,2*probRepCount);
    }
if(ps->pVal > 1 || ps->pVal < 0)
    warn("%s pVal is %.10f, wrong!");
slAddHead(&probePVals, ps);
}

doInvChiSq(int argc, char *argv[])
/* Test routine to make sure that gsl_cdf_chisq_P() is working correctly. */
{
int i = 0;
double result = 0;
for(i = 1; i < argc; i+=2)
    {
    result = gsl_cdf_chisq_P(atof(argv[i]), atoi(argv[i+1]));
    printf("%s\t%s\t%.2f\n", argv[i], argv[i+1], result);
    }
}

void printIntensities(struct row1lq *rowList, int limit)
/* Loop through and print out the row intensity for a list. */
{
struct row1lq *row = NULL;
int count = 0;
if(limit == -1)
    limit = slCount(rowList);

for(row = rowList; count < limit && row != NULL; row = row->next)
    {
    fprintf(stdout, "%.1f\n", row->intensity[repToSort]);
    count++;
    }
}

void fillInProbeSetPvals(double **matrix, int rowCount, int colCount, int repCount, int *repIndexes,
			 int sampleIx, struct row1lq **row1lqArray, double ***probeSetPValsP,
			 char ***probeSetNamesP, int *psCountP, double **probMatrix)
/* Calculate and fill in the probe set probabilities of expression given
   the data in the matrix. */
{
int i = 0, j = 0;
int rank = 0;
int binCount = 0;
struct psProb *probe = NULL;
struct row1lq *row = NULL;
int psCount = 0;

/* Fill in the values for the rows. */
for(i = 0; i < rowCount; i++)
    {
    row1lqArray[i]->repCount = repCount;
    for(j = 0; j < repCount; j++)
	{
	row1lqArray[i]->intensity[j] = matrix[i][repIndexes[j]];
	}
    }
/* Sort the gc bins by their rank and fill
   in the pVals as an empiracal rank. */
for(i = minGcBin; i <= maxGcBin; i++)
    {
    /* Sort for each replicate. */
    for(j = 0; j < repCount; j++)
	{
	rank = 1;
	repToSort = j;
	slSort(&gcBins[i-minGcBin], row1lqCmp);
	binCount = slCount(gcBins[i-minGcBin]);
	for(row = gcBins[i-minGcBin]; row != NULL; row = row->next)
	    {
	    row->pVal[j] = ((double)rank)/binCount;
	    probMatrix[xy2i(row->x, row->y)][repIndexes[j]] = row->pVal[j];
	    rank++;
	    }
	}
    }
/* Do the bad bin. */
for(row = badBin; row != NULL; row = row->next)
    {
    for(j = 0; j < repCount; j++)
	{
	row->pVal[j] = 0;
	}
    }
hashTraverseVals(probeSetHash, probeSetCalcPval);
if(outputPsSetSeqs == NULL)
    carefulClose(&outputPsSetSeqs);
slSort(&probePVals, psCmp);
psCount = slCount(probePVals);
/* If this is the first time through Allocate some memory. */
if(*probeSetPValsP == NULL)
    {
    (*psCountP) = psCount;
    AllocArray((*probeSetNamesP), psCount);
    AllocArray((*probeSetPValsP), psCount);
    for(i = 0; i < psCount; i++)
	{
	AllocArray((*probeSetPValsP)[i], sampleCount);
	}
    i = 0;
    for(probe = probePVals; probe != NULL; probe = probe->next)
	{
	(*probeSetNamesP)[i++] = cloneString(probe->psName);
	}
    }

i = 0;
for(probe = probePVals; probe != NULL; probe = probe->next)
    {
    (*probeSetPValsP)[i++][sampleIx] = probe->pVal;
    }
slFreeList(&probePVals);
}

struct replicateMap *createReplicateMap(char **expNames, int expCount)
/* Group replicates together if required or
   just make a replicate for each one if no replicates
   are required. */
{
struct replicateMap *rMapList = NULL, *rMap = NULL;
int i = 0, j = 0;
boolean combineReplicates = optionExists("combineReplicates");
boolean *used = NULL;
char expBuff[256], nameBuff[256];
char *mark = NULL;

/* If we're combining results look for
   expNames that have the same root (string before
   first ".") */
if(combineReplicates)
    {
    AllocArray(used, expCount);
    for(i = 0; i < expCount; i++)
	{
	/* If we haven't seen this experiment before
	   start a new replicateMap. */
	if(!used[i])
	    {
	    safef(nameBuff, sizeof(nameBuff), "%s", expNames[i]);
	    mark = strchr(nameBuff, '.');
	    if(mark != NULL)
		*mark = '\0';
	    AllocVar(rMap);
	    used[i] = TRUE;
	    rMap->rootName = cloneString(nameBuff);
	    rMap->repCount = 1;
	    AllocArray(rMap->repIndexes, expCount);
	    rMap->repIndexes[0] = i;

	    /* Loop through the experiments finding the ones
	       that are replicates. */
	    for(j = i+1; j < expCount; j++)
		{
		safef(nameBuff, sizeof(nameBuff), "%s", expNames[j]);
		mark = strchr(nameBuff, '.');
		if(mark != NULL)
		    *mark = '\0';
		if(sameString(nameBuff, rMap->rootName))
		    {
		    used[j] = TRUE;
		    rMap->repIndexes[rMap->repCount++] = j;
		    }
		}
	    maxReplicates = max(maxReplicates, rMap->repCount);
	    slAddHead(&rMapList, rMap);
	    }
	}
    }
else /* Make each expName a new replicate. */
    {
    for(i = 0; i < expCount; i++)
	{
	AllocVar(rMap);
	rMap->rootName = cloneString(expNames[i]);
	rMap->repCount = 1;
	AllocArray(rMap->repIndexes, rMap->repCount);
	rMap->repIndexes[0] = i;
	maxReplicates = max(maxReplicates, rMap->repCount);
	slAddHead(&rMapList, rMap);
	}
    }

slReverse(&rMapList);
/* Do a summmary of the replicates... */
/* for(rMap = rMapList; rMap != NULL; rMap = rMap->next) */
/*     { */
/*     printf("%s\t", rMap->rootName); */
/*     for(i = 0; i < rMap->repCount; i++) */
/* 	{ */
/* 	printf("%s,", expNames[rMap->repIndexes[i]]); */
/* 	} */
/*     printf("\n"); */
/*     } */
return rMapList;
}

boolean reportThisPSet(char *name)
/* Return TRUE if thes pset should be outputed for intensity and pvals. */
{
if(outputValsHash == NULL)
    return TRUE;
return hashLookup(outputValsHash, name) != NULL;
}

void outputPvalValues(struct row1lq **row1lqArray, double **probMatrix,
		      char **expNames, int rowCount, int colCount)
/* Write out the pvals of the probe sets over all the tissues. */
{
int rowIx = 0, colIx = 0;
struct row1lq *row = NULL;
assert(outputPvals);
/* Print column headers. */
fprintf(outputPvals, "%s\t", "probeSetName");
for(colIx = 0; colIx < colCount-1; colIx++)
    fprintf(outputPvals, "%s\t", expNames[colIx]);
fprintf(outputPvals
, "%s\n", expNames[colIx]);

/* Print all the data. */
for(rowIx = 0; rowIx < rowCount; rowIx++)
    {
    row = row1lqArray[rowIx];
    if(reportThisPSet(row->psName))
	{
	fprintf(outputPvals,"%d\t%s",rowIx, row->psName);
	for(colIx = 0; colIx < colCount; colIx++)
	    {
	    fprintf(outputPvals, "\t%.5f", probMatrix[rowIx][colIx]);
	    }
	fputc('\n', outputPvals);
	}
    }
}

void outputValues(double **matrix, struct row1lq **row1lqArray, char **expNames, int rowCount, int colCount)
/* Write out the intensity values for each probe. */
{
int rowIx = 0, colIx = 0;
struct row1lq *row = NULL;
assert(outputIntenVals);
fprintf(outputIntenVals, "%s\t", "probeSetName");
for(colIx = 0; colIx < colCount-1; colIx++)
    fprintf(outputIntenVals, "%s\t", expNames[colIx]);
fprintf(outputIntenVals, "%s\n", expNames[colIx]);

for(rowIx = 0; rowIx < rowCount; rowIx++)
    {
    row = row1lqArray[rowIx];
    if(reportThisPSet(row->psName))
	{
	assert(rowIx == xy2i(row->x, row->y));
 	fprintf(outputIntenVals, "%d\t%s", rowIx, row->psName);
/* 	fprintf(outputIntenVals, "%s\t%s\t%d\t%d\t%d\t%d", row->psName, row->seq, row->x, row->y, rowIx, xy2i(row->x, row->y)); */
	for(colIx = 0; colIx < colCount; colIx++)
	    fprintf(outputIntenVals, "\t%.1f", matrix[rowIx][colIx]);
	fputc('\n', outputIntenVals);
	}
    }
}

void gcPresAbs(char *outFile, char *file1lq, int celCount, char *celFiles[])
/* Output all of the probe sets and their pvals for expression as
   calculated from the intensities in the cel files. */
{
struct row1lq **rowArray = NULL, *row=NULL;
char **expNames = NULL, **expOutNames = NULL;
double **matrix = NULL, **probMatrix = NULL;
double **probeSetPVals = NULL; /* Matrix of pVals of expression with probe sets ordered altphabetically. */
char **probeSetNames = NULL;   /* List of probe set names indexing the probeSetPVals matrix. */
char *badProbeName = NULL;     /* Name of file with bad probes. */
int psCount = 0;
int i=0, j=0, k=0;
int rowCount=0;
FILE *out = NULL;
struct replicateMap *rMapList = NULL, *rMap = NULL;
boolean outputDuplicates = optionExists("outputDuplicates");
dotForUserInit(1);
initGcBins(8,17); /* Determined by where most of the data lives...
		     On the altmousea chip the data drops of significantly
		     before 8 (15,732 3% <= 8) and after 18 (9636 2% >=18) */

/* See if we have a kill list. */
badProbeName = optionVal("badProbes", NULL);
if(badProbeName != NULL)
    {
    warn("Reading the kill list");
    readBadProbes(badProbeName);
    }

/* Read in a 1lq file and put it in an array, lists and hash
   all at the same time. */
warn("Reading 1lq file: %s", file1lq);
read1lqFile(file1lq, &rowArray, &rowCount);

/* Lets initialize some memory. */
AllocArray(expNames, celCount);
AllocArray(expOutNames, celCount);
AllocArray(matrix, rowCount);
AllocArray(probMatrix, rowCount);
for(i=0; i<rowCount; i++)
    {
    AllocArray(matrix[i], celCount);
    AllocArray(probMatrix[i], celCount);
    }

/* Read in the cel files. */
for(i=0; i<celCount; i++)
    {
    char *lastSlash = NULL;
    dotForUser();
    expNames[i] = cloneString(celFiles[i]);
    lastSlash = strrchr(expNames[i], '/');
    if(lastSlash != NULL)
	expNames[i] = lastSlash+1;
    fillInMatrix(matrix, i, rowCount, celFiles[i]);
    }

if(optionExists("outputVals"))
    outputValues(matrix, rowArray, expNames, rowCount, celCount);
rMapList = createReplicateMap(expNames, celCount);
sampleCount = slCount(rMapList);
/* Allocate enough memory for the maximum number of replicates. */
for(i = 0; i < rowCount; i++)
    {
    AllocArray(rowArray[i]->intensity, maxReplicates);
    AllocArray(rowArray[i]->pVal, maxReplicates);
    }

warn("Calculating pVals");
//for(i=0; i<celCount; i++)
i = 0;
for(rMap = rMapList; rMap != NULL; rMap = rMap->next)
    {
    dotForUser();
    for(j = 0; j < rMap->repCount; j++)
	expOutNames[i+j] = cloneString(expNames[rMap->repIndexes[j]]);
    fillInProbeSetPvals(matrix, rowCount, celCount, rMap->repCount, rMap->repIndexes, i++, rowArray,
			&probeSetPVals, &probeSetNames, &psCount, probMatrix);
    }
if(optionExists("outputVals"))
    outputPvalValues(rowArray, probMatrix, expOutNames, rowCount, celCount);

warn("Done.");
warn("Found %d bad probes.", numBadProbes);
warn("Writing out the results.");
out = mustOpen(outFile, "w");

if(outputDuplicates)
    {
    /* Output each name, this only works if duplicates are
       next to each other. */
    for(rMap = rMapList; rMap->next != NULL; rMap = rMap->next)
	{
	for(i = 0; i < rMap->repCount; i++)
	    fprintf(out, "%s\t", expNames[rMap->repIndexes[i]]);
	}
    for(i = 0; i < rMap->repCount - 1; i++)
	fprintf(out, "%s\t", expNames[rMap->repIndexes[i]]);
    fprintf(out, "%s\n", expNames[rMap->repIndexes[i]]);

    for(i = 0; i < psCount; i++)
	{
	fprintf(out, "%s\t", probeSetNames[i]);
	for(j = 0; j < sampleCount; j++)
	    {
	    /* Output a column for each repetition. */
	    rMap = slElementFromIx(rMapList, j);
	    for(k = 0; k < rMap->repCount -1; k++)
		fprintf(out, "%.5g\t", probeSetPVals[i][j]);

	    /* If this is the last repetion then print the newline
	       otherwise continue with tabs. */
	    if(j == sampleCount -1)
		fprintf(out, "%.5g\n", probeSetPVals[i][j]);
	    else
		fprintf(out, "%.5g\t", probeSetPVals[i][j]);
	    }
	}
    }
else /* Only print one column per rMap. */
    {
    for(rMap = rMapList; rMap->next != NULL; rMap = rMap->next)
	fprintf(out, "%s\t", rMap->rootName);
    fprintf(out, "%s\n", rMap->rootName);

    for(i = 0; i < psCount; i++)
	{
	fprintf(out, "%s\t", probeSetNames[i]);
	for(j = 0; j < sampleCount -1; j++)
	    fprintf(out, "%.5g\t", probeSetPVals[i][j]);
	fprintf(out, "%.5g\n", probeSetPVals[i][j]);
	}
    }
carefulClose(&out);
warn("Done");
}

void initCountFiles(char *prefix)
/* open files for probe counts. */
{
struct dyString *dy = newDyString(1024);
dyStringClear(dy);
dyStringPrintf(dy, "%s.gcCounts.tab", prefix);
gcCountFile = mustOpen(dy->string, "w");

dyStringClear(dy);
dyStringPrintf(dy, "%s.gCounts.tab", prefix);
gCountsFile = mustOpen(dy->string, "w");

dyStringClear(dy);
dyStringPrintf(dy, "%s.cCounts.tab", prefix);
cCountsFile = mustOpen(dy->string, "w");
dyStringFree(&dy);
}

void closeCountFiles()
/* Close the counting file handles. */
{
carefulClose(&gcCountFile);
carefulClose(&gCountsFile);
carefulClose(&cCountsFile);
}

void setupTestRow(struct slRef **refList, double *array, int count)
/* Create a test row and add it to the list. */
{
int i = 0;
struct row1lq *row = NULL;
struct slRef *ref = NULL;
AllocVar(row);
row->repCount = count;
AllocArray(row->pVal, row->repCount);
for(i = 0; i < count; i++)
    row->pVal[i] = array[i];

AllocVar(ref);
ref->val = row;
slAddHead(refList, ref);
}

boolean sameDoubles(double a, double b)
/* return TRUE if two doubles are the same to 6 significant
   digits. */
{
char buffA[256];
char buffB[256];
safef(buffA, sizeof(buffA), "%.6f", a);
safef(buffB, sizeof(buffB), "%.6f", b);
return sameString(buffA, buffB);
}

void testMedian()
/* Quick testing of the median function. */
{
struct slRef *refList = NULL;
int count = 0, i = 0;
double first[] = {.7, .8, .9, .6, .5};
double second[] = {.8, .7, .6};
double third[] = {.6, .5};
double median = 0;
/* set up the first guy. */
setupTestRow(&refList, first, ArraySize(first));

median = medianRefList(refList);
if(!sameDoubles(median, .3))
    errAbort("Median should be .3, but it isn't.");
setupTestRow(&refList, second, ArraySize(second));

median = medianRefList(refList);
if(!sameDoubles(median, .3))
    errAbort("Median should still be .3, but it isn't.");
setupTestRow(&refList, third, ArraySize(third));
median = medianRefList(refList);
if(!sameDoubles(median, .35))
    errAbort("Median should be .7, but it isn't.");
warn("Passed all median tests.");
exit(0);
}

void initOutputValsList()
/* Read in the output vals list from the file and inititalize
   the outputValsHash. */
{
struct lineFile *lf = NULL;
char *s = NULL;
char *fileName = optionVal("outputValsList", NULL);
assert(fileName);
lf = lineFileOpen(fileName, TRUE);
outputValsHash = newHash(10);
while(lineFileNextReal(lf, &s))
    hashAdd(outputValsHash, s, NULL);
lineFileClose(&lf);
}

void initOutputVals()
/* Open two files for writing the probe level intensities
   and probabilities to. */
{
struct dyString *s = newDyString(512);
char *filePrefix = optionVal("outputVals", NULL);
assert(filePrefix);
dyStringPrintf(s, "%s.inten", filePrefix);
outputIntenVals = mustOpen(s->string, "w");
dyStringClear(s);
dyStringPrintf(s, "%s.pvals", filePrefix);
outputPvals = mustOpen(s->string, "w");
dyStringFree(&s);
}

int main(int argc, char *argv[])
/* Everybodys favorite function. */
{
char *affy1lqName = NULL;
char *outFileName = NULL;
int origCount = argc;
optionInit(&argc, argv, optionSpecs);
if(optionExists("testMedian"))
    testMedian();
if(optionExists("help") || origCount < 3)
    usage();
if(optionExists("countPrefix"))
    initCountFiles(optionVal("countPrefix", NULL));
if(optionExists("doMedianComb"))
    doMedianComb = TRUE;
if(optionExists("doCombo"))
    doCombo = TRUE;
if(optionExists("outputValsList"))
    initOutputValsList();
if(optionExists("outputVals"))
    initOutputVals();
if(optionExists("outputPsSetSeqs"))
    outputPsSetSeqs = mustOpen(optionVal("outputPsSetSeqs", NULL), "w");
affy1lqName = optionVal("1lqFile", NULL);
outFileName = optionVal("outFile", NULL);
if(optionExists("testVals"))
    doInvChiSq(argc, argv);
else if(affy1lqName == NULL || outFileName == NULL)
    errAbort("Need to specify 1lqFile and outFile, use -help for usage.");
else
    gcPresAbs(outFileName, affy1lqName, argc-1, argv+1);

if(optionExists("countPrefix"))
    closeCountFiles();
return 0;
}

max_stars_repo_name: andypohl/kent

max_issues_repo_name: andypohl/kent

max_forks_repo_name: andypohl/kent


text:
/*----------------------------------------------------------------------------*/
/*                                                                            */
/*  quad_golden.c                                                             */
/*                                                                            */
/*  Copyright (C) 2007 James Howse                                            */
/*  Copyright (C) 2009 Brian Gough                                             */
/*                                                                            */
/*  This program is free software; you can redistribute it and/or modify      */
/*  it under the terms of the GNU General Public License as published by      */
/*  the Free Software Foundation; either version 3 of the License, or (at     */
/*  your option) any later version.                                           */
/*                                                                            */
/*  This program is distributed in the hope that it will be useful, but       */
/*  WITHOUT ANY WARRANTY; without even the implied warranty of                */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU         */
/*  General Public License for more details.                                  */
/*                                                                            */
/*  You should have received a copy of the GNU General Public License         */
/*  along with this program; if not, write to the Free Software               */
/*  Foundation, Inc., 51 Franklin Street, Fifth Floor,                        */
/*  Boston, MA 02110-1301, USA.                                               */
/*                                                                            */
/*  ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::  */
/*                                                                            */
/*  This algorithm performs univariate minimization (i.e., line search).      */
/*  It requires only objective function values g(x) to compute the minimum.   */
/*  The algorithm maintains an interval of uncertainty [a,b] and a point x    */
/*  in the interval [a,b] such that a < x < b, and g(a) > g(x) and            */
/*  g(x) < g(b).  The algorithm also maintains the three points with the      */
/*  smallest objective values x, v and w such that g(x) < g(v) < g(w).  The   */
/*  algorithm terminates when max( x - a, b - x ) <  2(r |x| + t) where r     */
/*  and t are small positive reals.  At a given iteration, the algorithm      */
/*  first fits a quadratic through the three points (x, g(x)), (v, g(v))      */
/*  and (w, g(w)) and computes the location of the minimum u of the           */
/*  resulting quadratic.  If u is in the interval [a,b] then g(u) is          */
/*  computed.  If u is not in the interval [a,b], and either v < x and        */
/*  w < x, or v > x and w > x (i.e., the quadratic is extrapolating), then    */
/*  a point u' is computed using a safeguarding procedure and g(u') is        */
/*  computed.  If u is not in the interval [a,b], and the quadratic is not    */
/*  extrapolating, then a point u'' is computed using approximate golden      */
/*  section and g(u'') is computed.  After evaluating g() at the              */
/*  appropriate new point, a, b, x, v, and w are updated and the next         */
/*  iteration is performed.  The algorithm is based on work presented in      */
/*  the following references.                                                 */
/*                                                                            */
/*  Algorithms for Minimization without derivatives                           */
/*  Richard Brent                                                             */
/*  Prentice-Hall Inc., Englewood Cliffs, NJ, 1973                            */
/*                                                                            */
/*  Safeguarded Steplength Algorithms for Optimization using Descent Methods  */
/*  Philip E. Gill and Walter Murray                                          */
/*  Division of Numerical Analysis and Computing                              */
/*  National Physical Laboratory, Teddington, United Kingdom                  */
/*  NPL Report NAC 37, August 1974                                            */
/*                                                                            */
/*----------------------------------------------------------------------------*/
#include <config.h>

#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_min.h>

#include "min.h"

#define REL_ERR_VAL   1.0e-06
#define ABS_ERR_VAL   1.0e-10
#define GOLDEN_MEAN   0.3819660112501052	/* (3 - sqrt(5))/2 */
#define GOLDEN_RATIO  1.6180339887498950	/* (1 + sqrt(5))/2 */

#define DEBUG_PRINTF(x) /* do nothing */

typedef struct
{
  double step_size, stored_step, prev_stored_step;
  double x_prev_small, f_prev_small, x_small, f_small;
  unsigned int num_iter;
}
quad_golden_state_t;

static int
quad_golden_init (void *vstate, gsl_function * f, double x_minimum,
		  double f_minimum, double x_lower, double f_lower,
		  double x_upper, double f_upper)
{
  quad_golden_state_t *state = (quad_golden_state_t *) vstate;

  /* For the original behavior, the first value for x_minimum_minimum
     passed in by the user should be a golden section step but we
     don't enforce this here. */

  state->x_prev_small = x_minimum;
  state->x_small = x_minimum;

  state->f_prev_small = f_minimum;
  state->f_small = f_minimum;

  state->step_size = 0.0;
  state->stored_step = 0.0;
  state->prev_stored_step = 0.0;
  state->num_iter = 0;

  x_lower = 0 ;  /* avoid warnings about unused variables */
  x_upper = 0 ;
  f_lower = 0 ;
  f_upper = 0 ;
  f = 0;

  return GSL_SUCCESS;
}

static int
quad_golden_iterate (void *vstate, gsl_function * f, double *x_minimum,
		     double *f_minimum, double *x_lower, double *f_lower,
		     double *x_upper, double *f_upper)
{
  quad_golden_state_t *state = (quad_golden_state_t *) vstate;

  const double x_m = *x_minimum;
  const double f_m = *f_minimum;

  const double x_l = *x_lower;
  const double x_u = *x_upper;

  const double x_small = state->x_small;
  const double f_small = state->f_small;

  const double x_prev_small = state->x_prev_small;
  const double f_prev_small = state->f_prev_small;

  double stored_step = state->stored_step; /* update on exit */
  double prev_stored_step = state->prev_stored_step; /* update on exit */
  double step_size = state->step_size; /* update on exit */

  double quad_step_size = prev_stored_step;

  double x_trial;
  double x_eval, f_eval;

  double x_midpoint = 0.5 * (x_l + x_u);
  double tol = REL_ERR_VAL * fabs (x_m) + ABS_ERR_VAL; /* total error tolerance */

  if (fabs (stored_step) - tol > -2.0 * GSL_DBL_EPSILON)
    {
      /* Fit quadratic */
      double c3 = (x_m - x_small) * (f_m - f_prev_small);
      double c2 = (x_m - x_prev_small) * (f_m - f_small);
      double c1 = (x_m - x_prev_small) * c2 - (x_m - x_small) * c3;

      c2 = 2.0 * (c2 - c3);

      if (fabs (c2) > GSL_DBL_EPSILON)	/* if( c2 != 0 ) */
	{
	  if (c2 > 0.0)
	    c1 = -c1;

	  c2 = fabs (c2);

	  quad_step_size = c1 / c2;
	}
      else
	{
	  /* Handle case where c2 ~=~ 0  */
	  /* Insure that the line search will NOT take a quadratic
	     interpolation step in this iteration */
	  quad_step_size = stored_step;
	}

      prev_stored_step = stored_step;
      stored_step = step_size;
    }

  x_trial = x_m + quad_step_size;

  if (fabs (quad_step_size) < fabs (0.5 * prev_stored_step) && x_trial > x_l && x_trial < x_u)
    {
      /* Take quadratic interpolation step */
      step_size = quad_step_size;

      /* Do not evaluate function too close to x_l or x_u */
      if ((x_trial - x_l) < 2.0 * tol || (x_u - x_trial) < 2.0 * tol)
        {
          step_size = (x_midpoint >= x_m ? +1.0 : -1.0) * fabs(tol);
        }

      DEBUG_PRINTF(("quadratic step: %g\n", step_size));
    }
  else if ((x_small != x_prev_small && x_small < x_m && x_prev_small < x_m) ||
           (x_small != x_prev_small && x_small > x_m && x_prev_small > x_m))
    {
      /* Take safeguarded function comparison step */
      double outside_interval, inside_interval;

      if (x_small < x_m)
	{
	  outside_interval = x_l - x_m;
	  inside_interval = x_u - x_m;
	}
      else
	{
	  outside_interval = x_u - x_m;
	  inside_interval = x_l - x_m;
	}

      if (fabs (inside_interval) <= tol)
	{
          /* Swap inside and outside intervals */
          double tmp = outside_interval;
	  outside_interval = inside_interval;
	  inside_interval = tmp;
	}

      {
        double step = inside_interval;
        double scale_factor;

        if (fabs (outside_interval) < fabs (inside_interval))
          {
            scale_factor = 0.5 * sqrt (-outside_interval / inside_interval);
          }
        else
          {
            scale_factor = (5.0 / 11.0) * (0.1 - inside_interval / outside_interval);
          }

        state->stored_step = step;
        step_size = scale_factor * step;
      }

      DEBUG_PRINTF(("safeguard step: %g\n", step_size));
    }
  else
    {
      /* Take golden section step */
      double step;

      if (x_m < x_midpoint)
        {
          step = x_u - x_m;
        }
      else
        {
          step = x_l - x_m;
        }

      state->stored_step = step;
      step_size = GOLDEN_MEAN * step;

      DEBUG_PRINTF(("golden step: %g\n", step_size));
    }

  /* Do not evaluate function too close to x_minimum */
  if (fabs (step_size) > tol)
    {
      x_eval = x_m + step_size;
    }
  else
    {
      x_eval = x_m + (step_size >= 0 ? +1.0 : -1.0) * fabs(tol);
    }

  /* Evaluate function at the new point x_eval */
  SAFE_FUNC_CALL(f, x_eval, &f_eval);

  /* Update {x,f}_lower, {x,f}_upper, {x,f}_prev_small, {x,f}_small, and {x,f}_minimum */
  if (f_eval <= f_m)
    {
      if (x_eval < x_m)
	{
          *x_upper = x_m;
          *f_upper = f_m;
        }
      else
      	{
          *x_lower = x_m;
          *f_upper = f_m;
        }

      state->x_prev_small = x_small;
      state->f_prev_small = f_small;

      state->x_small = x_m;
      state->f_small = f_m;

      *x_minimum = x_eval;
      *f_minimum = f_eval;
    }
  else
    {
      if (x_eval < x_m)
        {
          *x_lower = x_eval;
          *f_lower = f_eval;
        }
      else
        {
          *x_upper = x_eval;
          *f_upper = f_eval;
        }

      if (f_eval <= f_small || fabs (x_small - x_m) < 2.0 * GSL_DBL_EPSILON)
	{
	  state->x_prev_small = x_small;
	  state->f_prev_small = f_small;

	  state->x_small = x_eval;
	  state->f_small = f_eval;
	}
      else if (f_eval <= f_prev_small ||
	       fabs (x_prev_small - x_m) < 2.0 * GSL_DBL_EPSILON ||
	       fabs (x_prev_small - x_small) < 2.0 * GSL_DBL_EPSILON)
	{
	  state->x_prev_small = x_eval;
	  state->f_prev_small = f_eval;
	}
    }

  /* Update stored values for next iteration */

  state->stored_step = stored_step;
  state->prev_stored_step = prev_stored_step;
  state->step_size = step_size;
  state->num_iter++;

  DEBUG_PRINTF(("[%d] Final State: %g  %g  %g\n", state->num_iter, x_l, x_m, x_u));

  return GSL_SUCCESS;
}

static const gsl_min_fminimizer_type quad_golden_type = { "quad-golden",	/* name */
  sizeof (quad_golden_state_t),
  &quad_golden_init,
  &quad_golden_iterate
};

const gsl_min_fminimizer_type *gsl_min_fminimizer_quad_golden =
  &quad_golden_type;

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#ifndef __GSL_MATRIX_USHORT_H__
#define __GSL_MATRIX_USHORT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_ushort.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned short * data;
  gsl_block_ushort * block;
  int owner;
} gsl_matrix_ushort;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_view;

typedef _gsl_matrix_ushort_view gsl_matrix_ushort_view;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_const_view;

typedef const _gsl_matrix_ushort_const_view gsl_matrix_ushort_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_ushort *
gsl_matrix_ushort_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_ushort *
gsl_matrix_ushort_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_block (gsl_block_ushort * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_matrix (gsl_matrix_ushort * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector_ushort *
gsl_vector_ushort_alloc_row_from_matrix (gsl_matrix_ushort * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector_ushort *
gsl_vector_ushort_alloc_col_from_matrix (gsl_matrix_ushort * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_ushort_free (gsl_matrix_ushort * m);

/* Views */

GSL_EXPORT
_gsl_matrix_ushort_view
gsl_matrix_ushort_submatrix (gsl_matrix_ushort * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_ushort_view
gsl_matrix_ushort_row (gsl_matrix_ushort * m, const size_t i);

GSL_EXPORT
_gsl_vector_ushort_view
gsl_matrix_ushort_column (gsl_matrix_ushort * m, const size_t j);

GSL_EXPORT
_gsl_vector_ushort_view
gsl_matrix_ushort_diagonal (gsl_matrix_ushort * m);

GSL_EXPORT
_gsl_vector_ushort_view
gsl_matrix_ushort_subdiagonal (gsl_matrix_ushort * m, const size_t k);

GSL_EXPORT
_gsl_vector_ushort_view
gsl_matrix_ushort_superdiagonal (gsl_matrix_ushort * m, const size_t k);

GSL_EXPORT
_gsl_matrix_ushort_view
gsl_matrix_ushort_view_array (unsigned short * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_ushort_view
gsl_matrix_ushort_view_array_with_tda (unsigned short * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector (gsl_vector_ushort * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector_with_tda (gsl_vector_ushort * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_submatrix (const gsl_matrix_ushort * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_row (const gsl_matrix_ushort * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_column (const gsl_matrix_ushort * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_diagonal (const gsl_matrix_ushort * m);

GSL_EXPORT
_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subdiagonal (const gsl_matrix_ushort * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_superdiagonal (const gsl_matrix_ushort * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array (const unsigned short * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array_with_tda (const unsigned short * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector (const gsl_vector_ushort * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector_with_tda (const gsl_vector_ushort * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT unsigned short   gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x);

GSL_EXPORT unsigned short * gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_EXPORT const unsigned short * gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_ushort_set_zero (gsl_matrix_ushort * m);
GSL_EXPORT void gsl_matrix_ushort_set_identity (gsl_matrix_ushort * m);
GSL_EXPORT void gsl_matrix_ushort_set_all (gsl_matrix_ushort * m, unsigned short x);

GSL_EXPORT int gsl_matrix_ushort_fread (FILE * stream, gsl_matrix_ushort * m) ;
GSL_EXPORT int gsl_matrix_ushort_fwrite (FILE * stream, const gsl_matrix_ushort * m) ;
GSL_EXPORT int gsl_matrix_ushort_fscanf (FILE * stream, gsl_matrix_ushort * m);
GSL_EXPORT int gsl_matrix_ushort_fprintf (FILE * stream, const gsl_matrix_ushort * m, const char * format);

GSL_EXPORT int gsl_matrix_ushort_memcpy(gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);
GSL_EXPORT int gsl_matrix_ushort_swap(gsl_matrix_ushort * m1, gsl_matrix_ushort * m2);

GSL_EXPORT int gsl_matrix_ushort_swap_rows(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_ushort_swap_columns(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_ushort_swap_rowcol(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_ushort_transpose (gsl_matrix_ushort * m);
GSL_EXPORT int gsl_matrix_ushort_transpose_memcpy (gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);

GSL_EXPORT unsigned short gsl_matrix_ushort_max (const gsl_matrix_ushort * m);
GSL_EXPORT unsigned short gsl_matrix_ushort_min (const gsl_matrix_ushort * m);
GSL_EXPORT void gsl_matrix_ushort_minmax (const gsl_matrix_ushort * m, unsigned short * min_out, unsigned short * max_out);

GSL_EXPORT void gsl_matrix_ushort_max_index (const gsl_matrix_ushort * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_ushort_min_index (const gsl_matrix_ushort * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_ushort_minmax_index (const gsl_matrix_ushort * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_ushort_isnull (const gsl_matrix_ushort * m);

GSL_EXPORT int gsl_matrix_ushort_add (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_EXPORT int gsl_matrix_ushort_sub (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_EXPORT int gsl_matrix_ushort_mul_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_EXPORT int gsl_matrix_ushort_div_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_EXPORT int gsl_matrix_ushort_scale (gsl_matrix_ushort * a, const double x);
GSL_EXPORT int gsl_matrix_ushort_add_constant (gsl_matrix_ushort * a, const double x);
GSL_EXPORT int gsl_matrix_ushort_add_diagonal (gsl_matrix_ushort * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_ushort_get_row(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t i);
GSL_EXPORT int gsl_matrix_ushort_get_col(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t j);
GSL_EXPORT int gsl_matrix_ushort_set_row(gsl_matrix_ushort * m, const size_t i, const gsl_vector_ushort * v);
GSL_EXPORT int gsl_matrix_ushort_set_col(gsl_matrix_ushort * m, const size_t j, const gsl_vector_ushort * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
unsigned short
gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
unsigned short *
gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (unsigned short *) (m->data + (i * m->tda + j)) ;
}

extern inline
const unsigned short *
gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const unsigned short *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_USHORT_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#define MPICH_SKIP_MPICXX 1
#define OMPI_SKIP_MPICXX 1
#include <Python.h>
#include <petsc.h>
#include "libpetsc4py/libpetsc4py.c"


text:
//
// Created by pedram pakseresht on 2/18/21.
//

static char help[] = "Simple Petsc program";

#include <petsc.h>
#include "incompressibleFlow.h"
#include "mesh.h"
#include "particleInertial.h"
#include "particles.h"
#include "particleInitializer.h"
#include "petscviewer.h"

static PetscErrorCode uniform_u(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *u, void *ctx) {
    //PetscInt d;
    //for (d = 0; d < Dim; ++d)
    //    u[d] = 1.0;
    //return 0;
    if (X[1]==1.0){
        u[0] = 1.0;
    }
    else {
        u[0] = 0.0;
    }
        u[1] =0.0;
    // u[0] = X[1] - X[1]*X[1] ;
    // u[1] = 0.0;
    return 0;
}

static PetscErrorCode uniform_u_t(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *u, void *ctx) {
    u[0] = 0.0;
    u[1] = 0.0;
    return 0;
}

static PetscErrorCode uniform_p(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *p, void *ctx) {
    p[0] = 0.0;
    return 0;
}

static PetscErrorCode uniform_T(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *T, void *ctx) {
    T[0] = 0.0;
    return 0;
}
static PetscErrorCode uniform_T_t(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *T, void *ctx) {
    T[0] = 0.0;
    return 0;
}

static PetscErrorCode SetInitialConditions(TS ts, Vec u) {

    PetscErrorCode (*initFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {uniform_u,uniform_p,uniform_T};
    // u here is the solution vector including velocity, temperature and pressure fields.
    DM dm;
    PetscReal t;
    PetscErrorCode ierr;
    PetscFunctionBegin;

    ierr = TSGetDM(ts, &dm);
    CHKERRQ(ierr);
    ierr = TSGetTime(ts, &t);

    CHKERRQ(ierr);

    ierr = DMProjectFunction(dm, 0.0, initFuncs, NULL, INSERT_ALL_VALUES, u);
    CHKERRQ(ierr);
    // get the flow to apply the completeFlowInitialization method
    ierr = IncompressibleFlow_CompleteFlowInitialization(dm, u);
    CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

static PetscErrorCode MonitorFlowAndParticleError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) {
    // PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
    // void *ctxs[3];
    DM dm;
    PetscDS ds;
    Vec v;
    PetscReal ferrors[3];
    PetscInt f;
    PetscErrorCode ierr;
    PetscInt num;

    PetscFunctionBeginUser;
    ierr = TSGetDM(ts, &dm);
    CHKERRQ(ierr);
    ierr = DMGetDS(dm, &ds);
    CHKERRQ(ierr);

    // get the particle data from the context
    ParticleData particlesData = (ParticleData)ctx;
    PetscInt particleCount;
    ierr = DMSwarmGetSize(particlesData->dm, &particleCount);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // compute the average particle location
    const PetscReal *coords;
    PetscInt dims;
    PetscReal avg[3] = {0.0, 0.0, 0.0};
    ierr = DMSwarmGetField(particlesData->dm, DMSwarmPICField_coor, &dims, NULL, (void **)&coords);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    for (PetscInt p = 0; p < particleCount; p++) {
        for (PetscInt n = 0; n < dims; n++) {
            avg[n] += coords[p * dims + n] / particleCount;  // PetscReal
        }
    }
    ierr = DMSwarmRestoreField(particlesData->dm, DMSwarmPICField_coor, &dims, NULL, (void **)&coords);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = PetscPrintf(PETSC_COMM_WORLD,
                       "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g, %2.3g] ParticleCount: %d\n",
                       (int)step,
                       (double)crtime,
                       (double)ferrors[0],
                       (double)ferrors[1],
                       (double)ferrors[2],
                       particleCount,
                       (double)avg[0],
                       (double)avg[1]);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Avg Particle Location: [%2.3g, %2.3g, %2.3g]\n", (double)avg[0], (double)avg[1], (double)avg[2]);

    CHKERRABORT(PETSC_COMM_WORLD, ierr);

     ierr = VecViewFromOptions(u, NULL, "-vec_view");
     CHKERRABORT(PETSC_COMM_WORLD, ierr);

    //  ierr = VecView(u, NULL);CHKERRABORT(PETSC_COMM_WORLD, ierr);
    // ierr = PetscPrintf(PETSC_COMM_WORLD, "current time=%g\n", crtime);CHKERRQ(ierr);

    ierr = DMSetOutputSequenceNumber(particlesData->dm, step, crtime);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

     /*
     if (step == 0 ){
        ierr = VecViewFromOptions(u, NULL, "-vec_init");
        CHKERRABORT(PETSC_COMM_WORLD, ierr);
     } else {
        ierr = VecViewFromOptions(u, NULL, "-vec_view");
         CHKERRABORT(PETSC_COMM_WORLD, ierr);
     }
     */

    if (step == 0) {
        ierr = ParticleViewFromOptions(particlesData, NULL, "-particle_init");
        CHKERRABORT(PETSC_COMM_WORLD, ierr);
    } else {
        ierr = ParticleViewFromOptions(particlesData, NULL, "-particle_view");
        CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }

    PetscFunctionReturn(0);
}

static PetscErrorCode ParticleSetInitialVelocity(DM particleDm) {

    PetscFunctionBeginUser;

    PetscErrorCode  ierr;

    /*
    PetscScalar *vel;
    PetscInt dim, Npb, p, n;

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Box Particle Initialization Options", NULL);CHKERRQ(ierr);
    ierr = PetscOptionsInt("-Npb", "The initial number of particles per box dimension", NULL, Npb, &Npb, NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd();

    ierr = PetscPrintf(PETSC_COMM_WORLD, "***Npb=%D\n***", Npb);CHKERRQ(ierr);
    ierr = DMSwarmGetField(particleDm, ParticleVelocity, NULL, NULL, (void **) &vel);CHKERRQ(ierr);

    ierr = DMGetDimension(particleDm, &dim);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "***dim=%D\n***", dim);CHKERRQ(ierr);

    Npb *= dim;
    for (p = 0; p < Npb; ++p) {
        for (n = 0; n < dim; n++) {
            vel[p+n] = 0.0;
        }
    }
    ierr = DMSwarmRestoreField(particleDm, ParticleVelocity, NULL, NULL, (void **) &vel);CHKERRQ(ierr);
    */

    Vec vel;
    ierr = DMSwarmCreateGlobalVectorFromField(particleDm,ParticleVelocity, &vel);CHKERRQ(ierr);
    ierr = VecSet(vel, 0.0);CHKERRQ(ierr);
    //ierr = VecView(vel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = DMSwarmDestroyGlobalVectorFromField(particleDm,ParticleVelocity, &vel);CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

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

    DM dm;     // domain definition
    TS ts;     // time-stepper
    PetscErrorCode ierr;

    PetscBag parameterBag; // constant flow parameters
    Vec flowField;         // flow solution vector
    FlowData flowData;

    PetscReal t;
    PetscInt Dim=2;
    PetscReal dt=0.1;       // dt for time stepper
    PetscInt max_steps=2; // maximum time steps

    PetscReal Re=1.0;       // Reynolds number
    PetscReal St=1.0;       // Strouhal number
    PetscReal Pe=1.0;       // Peclet number
    PetscReal Mu=1.0;      // viscosity
    PetscReal K_input=1.0;  // thermal conductivity
    PetscReal Cp = 1.0;    // heat capacity
    // PetscInt Np=10;

    // initialize Petsc ...
    ierr = PetscInitialize(&argc, &argv, NULL, "");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "*** Set up a simple flow *** \n");CHKERRQ(ierr);

    // setup the ts
    ierr = TSCreate(PETSC_COMM_WORLD, &ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = CreateMesh(PETSC_COMM_WORLD, &dm, PETSC_TRUE, Dim);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetDM(ts, dm);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    //output the mesh
    ierr = DMViewFromOptions(dm, NULL, "-dm_view");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup the flow data
    ierr = FlowCreate(&flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // setup problem
    ierr = IncompressibleFlow_SetupDiscretization(flowData, dm);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    IncompressibleFlowParameters flowParameters;

    // changing non-dimensional parameters manually here ...
    flowParameters.strouhal = St;
    flowParameters.reynolds = Re;
    flowParameters.peclet = Pe;
    flowParameters.mu = Mu;
    flowParameters.k = K_input;
    flowParameters.cp = Cp;

    // print out the parameters
    ierr = PetscPrintf(PETSC_COMM_WORLD, "*** non-dimensional parameters ***\n");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "St=%g\n", flowParameters.strouhal);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Re=%g\n", flowParameters.reynolds);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Pe=%g\n", flowParameters.peclet);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Mu=%g\n", flowParameters.mu);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "K=%g\n", flowParameters.k);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Cp=%g\n", flowParameters.cp);CHKERRQ(ierr);

    // Start the problem setup
    PetscScalar constants[TOTAL_INCOMPRESSIBLE_FLOW_PARAMETERS];
    ierr = IncompressibleFlow_PackParameters(&flowParameters, constants);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = IncompressibleFlow_StartProblemSetup(flowData, TOTAL_INCOMPRESSIBLE_FLOW_PARAMETERS, constants);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Override problem with source terms, boundary, and set the exact solution
    PetscDS prob;
    ierr = DMGetDS(dm, &prob);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup Boundary Conditions
    // Note: DM_BC_ESSENTIAL is a Dirichlet BC.
    PetscInt id;
    id = 3;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "top wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    id = 1;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "bottom wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 2;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "right wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 4;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "left wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 3;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "top wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 1;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "bottom wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, parameterBag);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 2;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "right wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, parameterBag);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 4;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "left wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, parameterBag);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = IncompressibleFlow_CompleteProblemSetup(flowData, ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Name the flow field
    ierr = PetscObjectSetName(((PetscObject)flowData->flowField), "Numerical Solution");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = SetInitialConditions(ts, flowData->flowField);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = TSGetTime(ts, &t);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    //  *** added by Pedram for dt and number of time steps ***
    ierr = TSSetTimeStep(ts,dt);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = TSSetMaxSteps(ts,max_steps);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    // ***********************************

    ierr = DMSetOutputSequenceNumber(dm, 0, t);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // checks for convergence ...
    ierr = DMTSCheckFromOptions(ts, flowData->flowField);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ParticleData particles;

    // Setup the particle domain
    ierr = ParticleInertialCreate(&particles, Dim);
    //ierr = ParticleInertialCreate(&particles, Dim);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // link the flow to the particles
    ierr = ParticleInitializeFlow(particles, flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // name the particle domain
    ierr = PetscObjectSetOptionsPrefix((PetscObject)(particles->dm), "particles_");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscObjectSetName((PetscObject)particles->dm, "Particles");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // initialize the particles position and velocity
    ierr = ParticleInitialize(dm, particles->dm);CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = ParticleSetInitialVelocity(particles->dm);CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // setup the flow monitor to also check particles
    ierr = TSMonitorSet(ts, MonitorFlowAndParticleError, particles, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetFromOptions(ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup particle position integrator
    TS particleTs;
    ierr = TSCreate(PETSC_COMM_WORLD, &particleTs);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscObjectSetOptionsPrefix((PetscObject)particleTs, "particle_");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = ParticleInertialSetupIntegrator(particles, particleTs, flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Solve the one way coupled system
    ierr = TSSolve(ts, flowData->flowField);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Cleanup
    ierr = DMDestroy(&dm);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSDestroy(&ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSDestroy(&particleTs);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = FlowDestroy(&flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = ParticleInertialDestroy(&particles);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscFinalize();
    exit(ierr);

    // _lower 0.25,0.25 _upper 0.75,0.75
    // Npb 10
    // _ts_dt 0.1
    // _layout_type box
    // _ts_convergence_estimate
    // _num_refine 1
    // -particle_ts_monitor_cancel
    // -ts_monitor_cancel

    // -particle_layout_type box
    // -particle_lower 0.1,0.1
    // -particle_upper 0.3,0.3
    // -Npb 5
    // -particle_ts_max_steps 10
    // -particle_ts_dt 0.05

    /*
    -ksp_type
    fgmres
    -ksp_gmres_restart
    10
    -ksp_rtol
    1.0e-9
    -ksp_error_if_not_converged
    -pc_type
    fieldsplit
    -pc_fieldsplit_0_fields
    0,2
      -pc_fieldsplit_1_fields
    1
    -pc_fieldsplit_type
    schur
    -pc_fieldsplit_schur_factorization_type
    full
    -fieldsplit_0_pc_type
    lu
    -fieldsplit_pressure_ksp_rtol
    1e-10
    -fieldsplit_pressure_pc_type
    jacobi
     */

    /*
     * -Npb 1 -particle_lower 0.05,0.1 -particle_upper 0.3,0.9 -particle_layout_type box -vec_init hdf5:/Users/pedram/scratch/sol.h5 -vec_view hdf5:/Users/pedram/scratch/sol.h5::append -dm_view hdf5:/Users/pedram/scratch/sol.h5 -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 -ksp_type fgmres -ksp_rtol 1.0e-9 -ksp_atol 1.0e-12 -ksp_error_if_not_converged -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_0_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -dm_plex_separate_marker -snes_monitor -ksp_monitor -snes_converged_reason -ksp_converged_reason -particle_view hdf5:/Users/pedram/scratch/solP.h5::append -particle_init hdf5:/Users/pedram/scratch/solP.h5
     */

}

max_stars_repo_name: pakserep/ablateClient

max_issues_repo_name: pakserep/ablateClient

max_forks_repo_name: pakserep/ablateClient


text:
/* Author:  Lisandro Dalcin   */
/* Contact: dalcinl@gmail.com */

#ifndef PETSC4PY_H
#define PETSC4PY_H

#include <Python.h>
#include <petsc.h>
#include "petsc4py.PETSc_api.h"

static int import_petsc4py(void) {
  if (import_petsc4py__PETSc() < 0) goto bad;
  return 0;
 bad:
  return -1;
}

#endif /* !PETSC4PY_H */


text:
#if !defined(MAIN_H_INCLUDED)
#define MAIN_H_INCLUDED

#include "Exceptions.h"

#include <cstdlib>
#include <filesystem>
#include <gsl/gsl>
#include <iostream>
#include <stdexcept>

template <typename T>
int commonMain(size_t argCount, const char*const*const argList)
{
	using ::std::filesystem::path;
	using ::gsl::make_span;
	using ::std::cout;
	using ::std::endl;

	auto exitCode{EXIT_FAILURE};
	try
	{
		// Pass the argument list, not including the program name (0th element):
		T program{make_span(argList + 1, argCount - 1)};
		exitCode = program.run();
	}
	catch (const CmdLineError& ex)
	{
		exitCode = T::usage(cout, path{argList[0]}.stem().generic_string(), ex.what());
	}
	catch (const ::std::exception& ex)
	{
		cout << endl
			<< "Exception of type " << typeid(ex).name() << endl
			<< "   " << ex.what() << endl
			<< endl;
	}
	return exitCode;
}

#endif // MAIN_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


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

int
main (void)
{
  gsl_multiset * c;
  size_t i;

  printf ("All multisets of {0,1,2,3} by size:\n") ;
  for (i = 0; i <= 4; i++)
    {
      c = gsl_multiset_calloc (4, i);
      do
        {
          printf ("{");
          gsl_multiset_fprintf (stdout, c, " %u");
          printf (" }\n");
        }
      while (gsl_multiset_next (c) == GSL_SUCCESS);
      gsl_multiset_free (c);
    }

  return 0;
}


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

static int fdjac(const gsl_vector *x, const gsl_vector *wts,
                 gsl_multifit_function_fdf *fdf,
                 const gsl_vector *f, gsl_matrix *J);

/*
fdjac()
  Compute approximate Jacobian using forward differences

Inputs: x   - parameter vector
        wts - data weights
        fdf - fdf struct
        f   - (input) vector of function values f_i(x)
        J   - (output) Jacobian matrix

Return: success or error
*/

static int
fdjac(const gsl_vector *x, const gsl_vector *wts,
      gsl_multifit_function_fdf *fdf, const gsl_vector *f, gsl_matrix *J)
{
  int status = 0;
  size_t i, j;
  double h;
  const double epsfcn = 0.0;
  double eps = sqrt(GSL_MAX(epsfcn, GSL_DBL_EPSILON));

  for (j = 0; j < fdf->p; ++j)
    {
      double xj = gsl_vector_get(x, j);

      /* use column j of J as temporary storage for f(x + dx) */
      gsl_vector_view v = gsl_matrix_column(J, j);

      h = eps * fabs(xj);
      if (h == 0.0)
        h = eps;

      /* perturb x_j to compute forward difference */
      gsl_vector_set((gsl_vector *) x, j, xj + h);

      status += gsl_multifit_eval_wf (fdf, x, wts, &v.vector);
      if (status)
        return status;

      /* restore x_j */
      gsl_vector_set((gsl_vector *) x, j, xj);

      h = 1.0 / h;
      for (i = 0; i < fdf->n; ++i)
        {
          double fnext = gsl_vector_get(&v.vector, i);
          double fi = gsl_vector_get(f, i);

          gsl_matrix_set(J, i, j, (fnext - fi) * h);
        }
    }

  return status;
} /* fdjac() */

/*
gsl_multifit_fdfsolver_dif_df()
  Compute approximate Jacobian using finite differences

Inputs: x   - parameter vector
        wts - data weights (set to NULL if not needed)
        fdf - fdf
        f   - (input) function values f_i(x)
        J   - (output) approximate Jacobian matrix

Return: success or error
*/

int
gsl_multifit_fdfsolver_dif_df(const gsl_vector *x, const gsl_vector *wts,
                              gsl_multifit_function_fdf *fdf,
                              const gsl_vector *f, gsl_matrix *J)
{
  return fdjac(x, wts, fdf, f, J);
} /* gsl_multifit_fdfsolver_dif_df() */

#ifndef GSL_DISABLE_DEPRECATED

/*
gsl_multifit_fdfsolver_dif_fdf()
  Compute function values (analytic) and approximate Jacobian using finite
differences

Inputs: x      - parameter vector
        fdf    - fdf
        f      - (output) function values f_i(x)
        J      - (output) approximate Jacobian matrix

Return: success or error
*/

int
gsl_multifit_fdfsolver_dif_fdf(const gsl_vector *x,
                               gsl_multifit_function_fdf *fdf,
                               gsl_vector *f, gsl_matrix *J)
{
  int status = 0;

  status = gsl_multifit_eval_wf(fdf, x, NULL, f);
  if (status)
    return status;

  status = fdjac(x, NULL, fdf, f, J);
  if (status)
    return status;

  return status;
} /* gsl_multifit_fdfsolver_dif_fdf() */

#endif /* !GSL_DISABLE_DEPRECATED */


text:
#pragma once
#include "parser.h"
#include "format.h"
#include "nameutils.h"
#include "utils.h"
#include <sfun/string_utils.h>
#include <cmdlime/errors.h>
#include <gsl/gsl>
#include <algorithm>
#include <sstream>
#include <iomanip>
#include <functional>
#include <cassert>

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

template <FormatType formatType>
class PosixParser : public Parser<formatType>
{
    using Parser<formatType>::Parser;

    void processCommand(std::string command)
    {
        auto possibleNumberArg = command;
        command = str::after(command, "-");
        if (isParamOrFlag(command)){
            if (this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands){
                if (!foundParam_.empty())
                    throw ParsingError{"Parameter '-" + foundParam_ + "' value can't be empty"};
                if (argumentEncountered_)
                    throw ParsingError{"Flags and parameters must precede arguments"};
            }
            parseCommand(command);
        }
        else if (isNumber(possibleNumberArg)){
            this->readArg(possibleNumberArg);
            argumentEncountered_ = true;
        }
        else if (this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands)
            throw ParsingError{"Encountered unknown parameter or flag '-" + command + "'"};
    }

    void preProcess() override
    {
        checkNames();
        argumentEncountered_ = false;
        foundParam_.clear();
    }

    void process(const std::string& token) override
    {
        if (str::startsWith(token, "-") && token.size() > 1)
           processCommand(token);
        else if (!foundParam_.empty()){
            this->readParam(foundParam_, token);
            foundParam_.clear();
        }
        else{
            this->readArg(token);
            argumentEncountered_ = true;
        }
    }

    void postProcess() override
    {
        if (!foundParam_.empty())
            throw ParsingError{"Parameter '-" + foundParam_ + "' value can't be empty"};
    }

    void parseCommand(const std::string& command)
    {
        if (command.empty())
            throw ParsingError{"Flags and parameters must have a name"};
        auto paramValue = std::string{};
        for(auto ch : command){
            auto opt = std::string{ch};
            if (!foundParam_.empty())
                paramValue += opt;
            else if (this->findFlag(opt))
                this->readFlag(opt);
            else if (this->findParam(opt) || this->findParamList(opt))
                foundParam_ = opt;
            else if (this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands)
                throw ParsingError{"Unknown option '" + opt + "' in command '-" + command + "'"};
        }
        if (!foundParam_.empty() && !paramValue.empty()){
            this->readParam(foundParam_, paramValue);
            foundParam_.clear();
        }
    }

    void checkNames()
    {
        auto check = [](const OptionInfo& var, const std::string& varType){
            if (var.name().size() != 1)
                throw ConfigError{varType + "'s name '" + var.name() + "' can't have more than one symbol"};
            if (!std::isalnum(var.name().front()))
                throw ConfigError{varType + "'s name '" + var.name() + "' must be an alphanumeric character"};
        };
        this->forEachParamInfo([check](const OptionInfo& var){
            check(var, "Parameter");
        });
        this->forEachParamListInfo([check](const OptionInfo& var){
            check(var, "Parameter");
        });
        this->forEachFlagInfo([check](const OptionInfo& var){
            check(var, "Flag");
        });
    }

    bool isParamOrFlag(const std::string& str)
    {
        if (str.empty())
            return false;
        auto opt = str.substr(0,1);
        return this->findFlag(opt) ||
               this->findParam(opt) ||
               this->findParamList(opt);
    }

private:
    bool argumentEncountered_ = false;
    std::string foundParam_;
};

class PosixNameProvider{
public:
    static std::string name(const std::string& optionName)
    {
        Expects(!optionName.empty());
        return std::string{static_cast<char>(std::tolower(optionName.front()))};
    }

    static std::string shortName(const std::string& optionName)
    {
        Expects(!optionName.empty());
        return {};
    }

    static std::string fullName(const std::string& optionName)
    {
        Expects(!optionName.empty());
        return toKebabCase(optionName);
    }

    static std::string valueName(const std::string& typeName)
    {
        Expects(!typeName.empty());
        return toCamelCase(templateType(typeNameWithoutNamespace(typeName)));
    }
};

class PosixOutputFormatter{
public:
    static std::string paramUsageName(const IParam& param)
    {
        auto stream = std::stringstream{};
        if (param.isOptional())
            stream << "[" << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">]";
        else
            stream << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">";
        return stream.str();
    }

    static std::string paramListUsageName(const IParamList& param)
    {
        auto stream = std::stringstream{};
        if (param.isOptional())
            stream << "[" << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">...]";
        else
            stream << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">...";
        return stream.str();
    }

    static std::string paramDescriptionName(const IParam& param, int indent = 0)
    {
        auto stream = std::stringstream{};
        stream << std::setw(indent) << paramPrefix()
               << param.info().name() << " <" << param.info().valueName() << ">";
        return stream.str();
    }

    static std::string paramListDescriptionName(const IParamList& param, int indent = 0)
    {
        auto stream = std::stringstream{};
        stream << std::setw(indent) << paramPrefix()
               << param.info().name() << " <" << param.info().valueName() << ">";
        return stream.str();
    }

    static std::string paramPrefix()
    {
        return "-";
    }

    static std::string flagUsageName(const IFlag& flag)
    {
        auto stream = std::stringstream{};
        stream << "[" << flagPrefix() << flag.info().name() << "]";
        return stream.str();
    }

    static std::string flagDescriptionName(const IFlag& flag, int indent = 0)
    {
        auto stream = std::stringstream{};
        stream << std::setw(indent) << flagPrefix() << flag.info().name();
        return stream.str();
    }

    static std::string flagPrefix()
    {
        return "-";
    }

    static std::string argUsageName(const IArg& arg)
    {
        auto stream = std::stringstream{};
        stream << "<" << arg.info().name() << ">";
        return stream.str();
    }

    static std::string argDescriptionName(const IArg& arg, int indent = 0)
    {
        auto stream = std::stringstream{};
        if (indent)
            stream << std::setw(indent) << " ";
        stream << "<" << arg.info().name() << "> (" << arg.info().valueName() << ")";
        return stream.str();
    }

    static std::string argListUsageName(const IArgList& argList)
    {
        auto stream = std::stringstream{};
        if (argList.isOptional())
            stream << "[" << argList.info().name() << "...]";
        else
            stream << "<" << argList.info().name() << "...>";
        return stream.str();
    }

    static std::string argListDescriptionName(const IArgList& argList, int indent = 0)
    {
        auto stream = std::stringstream{};
        if (indent)
            stream << std::setw(indent) << " ";
        stream << "<" << argList.info().name() << "> (" << argList.info().valueName() << ")";
        return stream.str();
    }

};

template<>
struct Format<FormatType::POSIX>
{
    using parser = PosixParser<FormatType::POSIX>;
    using nameProvider = PosixNameProvider;
    using outputFormatter = PosixOutputFormatter;
    static constexpr bool shortNamesEnabled = false;
};

}

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


text:
/** \file IRLS.h
*
* `IRLS' is a C++ implementation of the IRLS algorithm for GLM
* Copyright (C) 2013 Xioaquan Wen, Timothee Flutre
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef _IRLS_H_
#define _IRLS_H_

#include <vector>

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

#include "LinkFunc.h"

class IRLS
{
private:
  gsl_vector * y;      // response vector
  gsl_matrix * X;      // design matrix
  gsl_vector * offset; // offset vector
  bool free_data;      // depends on load_data() or set_data()

  size_t n;            // sample size;
  size_t p;            // number of parameters (including intercept)
  size_t rank;         // of X (useful for p-values)

  gsl_vector * bv;     // vector of estimated effect sizes
  gsl_matrix * VB;     // covariance matrix of estimated effect sizes
  double psi;          // dispersion

public:

  LinkFunc * link;

  IRLS(const char * link_type);
  ~IRLS();
  void load_data(const std::vector<double> & yv,
		 const std::vector<std::vector<double> > &Xv,
		 const std::vector<double> & offv);
  void set_data(gsl_vector * yv,
		gsl_matrix * Xv,
		gsl_vector * offv);
  void fit_model();
  std::vector<double> get_coef();
  std::vector<double> get_stderr();
  size_t get_rank_X() { return rank; };
  double get_dispersion() { return psi; };

private:
  void compute_variance(gsl_vector *w);

}; // IRLS

#endif // _IRLS_H_

max_stars_repo_name: Lan-lab/Chrom-Lasso-

max_issues_repo_name: Lan-lab/Chrom-Lasso-

max_forks_repo_name: Lan-lab/Chrom-Lasso-


text:
/**
 * Cascade simulation based on pair production and bremsstrahlung interactions
 *
 * copyright  (C) 2007
 * the icecube collaboration
 *
 * $Id: I3CascadeSimulation.h 35916 2007-08-23 09:34:12Z bvoigt $
 *
 * @version $Revision: 35916 $
 * @date $LastChangedDate: 2007-08-23 03:34:12 -0600 (Thu, 23 Aug 2007) $
 * @author Bernhard Voigt <bernhard.voigt@desy.de>   Last changed by: $LastChangedBy: bvoigt $
 */

#ifndef I3_CASCADE_SIMULATION_H_INCLUDED
#define I3_CASCADE_SIMULATION_H_INCLUDED

//Standard C/C++ includes
#include <vector>
#include <stack>

// Gnu Scientific Library
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>

// local includes
#include "cmc/I3CascadeMCCommon.h"
#include "cmc/I3CascadeDevelopment.h"
#include "cmc/I3CascadeParametrization.h"
#include "cmc/I3MetropolisHastings.h"

// icetray includes
#include "phys-services/I3RandomService.h"

/**
 * @brief Local definition of a Particle struct
 *
 * The particle struct holds energy and postion informations
 * of particles in the simulation.
 */
struct Particle {
  static const short ELECTRON = 1;
  static const short PHOTON = 0;
  double energy;
  double x;
  short type;
  /**
   * Particle with energy, at position x (1-dimensional), of
   * type type (1 is electron, 0 is photon)
   *
   * @param energy
   * @param x
   * @param type
   */
  Particle(double energy, double x, short type) : energy(energy), x(x), type(type) { }
};

/**
 * @brief Local definition of beta distribution
 *
 * The struct defines a pdf function, ie. the beta distribtuion
 * not normed, though.
 * It provides a rvs method to draw random samples from
 * the beta distribution using gsl_rng
 *
 * The parameters used for the beta distribibution (alpha and beta)
 * are defined.
 *
 * Functions are passed to I3MetropolisHastings and need to be static.
 */
struct beta_dist {
  // define the logger name, use the same as for the simulation class
  SET_LOGGER("I3CascadeSimulation");

  // parameter to beta function
  static double alpha;
  static double beta;

  /**
   * Beta Function \f$ x^{\alpha-1} (1-x)^{\beta-1}\f$
   *
   * @param x
   */
  static double pdf(const double& x) {
    log_trace("alpha %.2f, beta %.2f", alpha, beta);
    return pow(x, alpha-1) * pow(1-x, beta-1);
  }

  /**
   * Returns a random variate drawn from the beta distribution
   *
   * Uses the gsl_ran_beta method
   *
   * @param random gsl_rng instance
   */
  static double rvs(const gsl_rng* random) {
    log_trace("alpha %.2f, beta %.2f", alpha, beta);
    return gsl_ran_beta(random, alpha, beta);
  }
};

/**
 *@brief This class implements a Cascade development simulation including
 *       the LPM effect
 *
 * Description of the simulation can be found here:
 * http://wiki.icecube.wisc.edu/index.php/Cascade_Simulation
 *
 * Cross section formulars used are defined in I3CascadeSimulationCrossSections
 *
 *
 * @author Bernhard Voigt
 */
class I3CascadeSimulation : public I3CascadeDevelopment, public I3CascadeMCCommon::Sampler {

  // define the logger name
  SET_LOGGER("I3CascadeSimulation");

 public:

  /**
   * @brief Constructor
   *
   * @param random random number service from icetray frame
   */
  I3CascadeSimulation(I3RandomServicePtr random);

  /**
   * @brief Destructor
   *
   * Frees the gsl objects
   */
  ~I3CascadeSimulation();

  virtual void SetRandomNumberGenerator(I3RandomServicePtr);

  /**
   * Simulates a cascade of the given energy
   *
   * @param energy
   */
  void Simulate(I3Particle::ParticleType type, double energy);

  /**
   * Sets the energy threshold down to which particles are tracked (Default is 50 TeV)
   *
   * This shouldn't be lower than 1 GeV, since this is the limit for the Bremsstrahlung
   * interaction
   *
   * @param thresholdEnergy
   */
  void SetThreshold(double thresholdEnergy) {
    threshold_ = thresholdEnergy;
  }

  /**
   * Returns the threshold energy down to which particles are tracked
   *
   * @return threshold energy
   */
  double GetThreshold() {
    return threshold_;
  }

  /**
   * @brief default energy threshold down to which particles are tracked
   */
  static const double DEFAULT_THRESHOLD;

 private:

  I3CascadeSimulation(const I3CascadeSimulation&);

  /**
   * Cacluates the pair productio mean free path for the given energy
   *
   * @param energy
   */
  double PairProductionMeanFreePath(double energy);

  /**
   * Cacluates the bremsstrahlung radiation length (cut off energy is defined as 1 GeV)
   *
   * @param energy
   */
  double BremsstrahlungMeanFreePath(double energy);

  /**
   * Samples the bremsstrahlung cross section to get a fractional energy
   * of a secondary particle produced in an interaction with a incident particle
   * of the given energy
   *
   * @see I3MetropolisHastings::Sample
   *
   * @param energy
   */
  double SampleBremsstrahlungCrossSection(double energy);

  /**
   * Samples the pair production cross section to get a fractional energy
   * of a secondary particle produced in an interaction with a incident particle
   * of the given energy
   *
   * @see I3MetropolisHastings::Sample
   *
   * @param energy
   */
  double SamplePairProductionCrossSection(double energy);

  /**
   * Inits the simulation
   *
   * Interpolation of the mean free paths is performed
   * The metropolis hastings random samplers are initialized
   */
  void Init();

  /**
   * Inits the metropolis hastings random samplers
   */
  void InitMetropolisHastings();

  /**
   * Draw and discard samples from the MH samplers
   * to make them independent of initial conditions
   */
  void BurnIn();

  /**
   * Calcuation and interpolation of the pair production mean free path
   * in an energy range specified by lower and upper (in log10)
   * The interpolation used the given number of points.
   * Default energy range is from 100 GeV to 10 EeV
   *
   * Uses the integration and interpolation routines from gsl
   *
   * @param lower
   * @param upper
   * @param points
   */
  void InterpolPairProductionMeanFreePath(double lower=2., double upper=13.0, int points=100);

  /**
   * Calcuation and interpolation of the bremsstrahlung radiation length
   * in an energy range specified by lower and upper (in log10)
   * The interpolation used the given number of points.
   * Default energy range is from 100 GeV to 10 EeV
   *
   * Uses the integration and interpolation routines from gsl
   *
   * The integration of the cross section is done down to a lower edge of
   * I3CascadeSimulation::BREMSSTRAHLUNGCUTOFF
   * this is due to the infrared limit, where the cross section raises to infinite values.
   *
   * @param lower
   * @param upper
   * @param points
   */
  void InterpolBremsstrahlungMeanFreePath(double lower=2., double upper=13.0, int points=100);

  /**
   * Returns the mean free path of a particle
   *
   * @param particle
   */
  double MeanFreePath(Particle particle);

  /**
   * Propagates a particle
   *
   * Draws the next interaction point
   * Calls the interaction method to produce secondaries
   *
   * @param particle
   */
  void Propagate(Particle particle);

  /**
   * Creates secondaries produced in an interaction of the given particle
   *
   * @param particle
   */
  void Interaction(Particle particle);

  /**
   * Produces a secondary gamma
   */
  void BremsstrahlungInteraction(Particle electron);

  /**
   * Produces a electron, positron pair
   */
  void PairProductionInteraction(Particle photon);

  /**
   * Deposits the particle
   *
   * Calculates the energy loss profile of the given particle
   * and adds it to the total energy loss profile of this simulation
   */
  void Deposit(Particle particle);

  /**
   * @brief stack holding particles that have to be propagated
   */
  std::stack<Particle> particles_;

  /**
   * @brief counter how many particles have been created
   */
  unsigned int particlesCreated_;

  /**
   * @brief counter how many particles have been deleted
   */
  unsigned int particlesDeleted_;

  /**
   * @brief struct holding beta distribution parameters, provides function
   * evaluation and random numbers
   */
  beta_dist betaDist_;

  /**
   * @brief random number sampler for low energy bremsstrahlung
   */
  I3MetropolisHastings lowBremsSampler_;

  /**
   * @brief random number sampler for high energy bremsstrahlung
   */
  I3MetropolisHastings highBremsSampler_;

  /**
   * @brief random number sampler for low energy pair production
   */
  I3MetropolisHastings lowPairSampler_;

  /**
   * @brief random number sampler for high energy pair production
   */
  I3MetropolisHastings highPairSampler_;

  /**
   * @brief parametrization of energy loss profile for low energy particles
   */
  I3CascadeParametrization parametrization_;

  /**
   * @brief gsl spline accelerator object used for spline evaluation (bremsstrahlung)
   */
  gsl_interp_accel* interpolAccBrems_;

  /**
   * @brief gsl spline - interpolation of Bresmstrahlung mean free path
   */
  gsl_spline* splineBrems_;

  /**
   * @brief gsl spline accelerator object used for spline evaluation (pair production)
   */
  gsl_interp_accel* interpolAccPair_;

  /**
   * @brief gsl spline - interpolation of Pairproduction mean free path
   */
  gsl_spline* splinePair_;

  /**
   * @brief threshold down to which particles are tracked, default is 50 TeV
   */
  double threshold_;

  /**
   * @brief repeat burn-in for each event
   *
   * This ensures that the events are truly independent and can be reproduced
   * given the state of the random number generator.
   */
  bool perEventBurnIn_;

  /**
   * @brief lowest energy of bremsstrahlung photons
   */
  static const double BREMSSTRAHLUNGCUTOFF;

};
#endif