New upstream version 3.2.1+dfsg

/*  

 *  SweeD

 *

 *  Copyright September 2012 by Nikolaos Alachiotis and Pavlos Pavlidis

 *

 *  This program is free software; you may redistribute it and/or modify its

 *  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.

 *

 *  For any other enquiries send an email to

 *  Pavlos Pavlidis (pavlidisp@gmail.com) or

 *  Nikolaos Alachiotis (n.alachiotis@gmail.com)

 *  

 */

#ifndef _SweeD_H

#define _SweeD_H

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <assert.h>

#include <math.h>

#include <time.h>

#include <sys/time.h>

#include <limits.h>

#include "SweeD_BFGS.h"

#ifdef _USE_PTHREADS

#include <pthread.h>

#include <unistd.h>

#define EXIT 127

#define BUSYWAIT 0

#define LIKSFS 1

#define COMPRVLUT 2

#define COMPALPHA 3

#endif

#ifdef _DO_CHECKPOINTS

#include "dmtcpaware.h"

#endif

#ifdef _ANALYTICAL_SFS

#include <mpfr.h>

#include <float.h>

#define ACC 10000

#define APPROX 0

#define TOLTIME 1e-7

#define RESCALE 1

#endif

#define ENTRIES_PER_INT 32

#define DENOMINATOR_OFFSET 0.00001

#define REGISTER_WIDTH 32

#define MAXSIZE 10000

#define INFILENAMESIZE 1000

#define SEQNAMESIZE 1000

#define DECREASE 0.99

#define MINSNPPERWINDOW 5

#define STRINGLENGTH 100

#define VCF_HLENGTH 9 // number of fields in the VCF header line

#define MAX_CHROM_NAME_VCF 100

#define MAX_STATES_VCF 5

#define RESULTS_DEFAULT 0

#define RESULTS_ALL 1

#define RSQUARE 0

#define DOM 1

#define ABSDOM 2

#define JUSTD 3

#define ABSD 4

#define ABSDOM2 5

#define ZERO '0'

#define ONE  '1'

#define GAP '-'

#define AD 'A'

#define CY 'C'

#define GU 'G'

#define TH 'T'

#define UN 'N'

#define BINARY 2

#define BINARY_WITH_GAPS 3

#define DNA 4

#define DNA_WITH_GAPS 5

#define STATESALL 8

#define OTHER_FORMAT -1

#define MS_FORMAT 0

#define FASTA_FORMAT 1

#define MACS_FORMAT 2

#define VCF_FORMAT 3

#define SF_FORMAT 4

#define INTOKS 8

#define N_REF 2

#define N_ALT 3

#define N_INFO 6

#define MAX_STATES 100

#define MAXINT 2147483647

char bits_in_16bits [0x1u << 16];

int linkage_disequilibrium;

int borderTol;

int VCF_header_lines;

char VCF_alignment_name [MAX_CHROM_NAME_VCF];

int VCF_first_SNP;

int nxtVCFalignment;

char runName[INFILENAMESIZE];

typedef float cor_t;

cor_t ABS (cor_t input);

double mainTime0;

int map[INTOKS];

int nofTokens;

int nofSamples;

typedef int al_t;

typedef double t_sfs;

struct pattern

{

	int x;

	int n;

	unsigned char f;

	int c;

	struct pattern * nxt;		

};

typedef struct 

{

 	al_t	*n; // number of valid bases in every site

	al_t	*n_t; // number of valid bases in every site

	al_t *x; // number of derived bases in every site

	al_t *x_t; // number of derived bases in every site

	float 	       * positions;

	int 	       * positionsInd;

	int 	       * p_t;

	unsigned char *folded;

	unsigned char *f_t;

	unsigned char userSetFolded;

	unsigned char userSFS;

	char* outgroupSequence;

  	char* outgroupName;

	int outgroupIndex;

	int sequences;

  	int foldedFormat; // if outgroup exists we set this to zero

	int 		 segsites;	

	int		 length;

	int minn;

	int minx;	

  	int maxn;

	int maxx;

	al_t * stateA;

	al_t * stateC;

	al_t * stateG;

	al_t * stateT;

	t_sfs * SFS;

	t_sfs * tmpSFS;

	int startSFS;

	int endSFS;

	int * ng;

	t_sfs * lb;

	t_sfs * ub;

	int params;

	al_t * n_pat;

	al_t * x_pat;

	unsigned char * f_pat;

	int *  c_pat;

	int patterns;

	t_sfs *** gridProbs;

	t_sfs * gridADs;

	t_sfs interval;

	t_sfs logAD0;	

	t_sfs * baseLikelihood;

	struct pattern * patternList;

	float sweepStep;

	int discSites;

#ifdef _USE_PTHREADS

	t_sfs * threadScoresLIKSFS;

#endif

} alignment_struct;

alignment_struct * alignment;

double * rv;

double ** rvLUT;

double ** rvLUT2;

#define GRIDSIZE 300

#define LOCALGRIDSIZE 100

double div6; 

typedef struct

{

	double * val;

	int * sc;

} factorial_LUT;

factorial_LUT * factLUT;

typedef struct

{

  int 		sfRealPos; 

  t_sfs		likelihood; 

  t_sfs		alpha; 

} clr_struct;

clr_struct * clr;

#ifdef _ANALYTICAL_SFS

typedef struct

{

	double t;

	double d;

} event;

int eventsTotal;

event * eventList;

mpfr_t ** analyticalSFS_LUT;

void computeAnalyticalSFS(int sequences, FILE * fp);

void computeAnalyticalSFSExpo(int sequences, double R, FILE * fpSFSo);

long double xExponential_Integral_Ei( long double x );

void mpfr_Exponential_Integral_Ei( mpfr_t EImp, mpfr_t xmp );

#endif

void printHeading (FILE * fp);

void computeCorrelationMatrixPairwise(alignment_struct * alignment, clr_struct * omega, int omegaIndex, int firstRowIndex, void * threadData, cor_t ** myCorrelationMatrix, char * lookuptable);

void applyCorrelationMatrixAdditions (clr_struct * omega, int omegaIndex, int firstRowIndex, cor_t ** correlationMatrix);

void overlapCorrelationMatrixAdditions (alignment_struct * alignment, clr_struct * omega, int lvw_i, int cvw_i, int * firstRowToCopy, int * firstRowToCompute, int * firstRowToAdd);

void shiftCorrelationMatrixValues (clr_struct * omega, int lvw_i, int cvw_i, int firstRowToCopy, cor_t ** correlationMatrix);

void commandLineParser(int argc, char** argv, 

		       char * infile,

		       char * sfsfile,

                       char * sfsofile,

		       char * sfofile, 

		       int * grid, 

                       int * length,                    

                       unsigned int * seed,

		       int * fileFormat,	

		       int * results,

		       char * outgroupName,

		       unsigned char *userSetFolded,

			unsigned char * userSFS,

			unsigned char * onlySFS,

			unsigned char * onlySF,

		       int * monomorphic,

			int * threads,

			int * sequences,

		       int * analyticalSFS,

		       double *growthRate);

void removeMonomorphicSites (int monomorphic, FILE * fp);

t_sfs computeBaseLikelihood();

t_sfs getAlpha (int sweepPosition, t_sfs * likelihood);

int isBinary(char input);

void printIntArray(int *array, int n);

void printUCharArray(unsigned char *array, int n);

void printFloatArray(float *array, int n);

void printCharArray(char *array, int n);

int isEndOfLine(char ent);

void goToLineEnd(FILE *fp);

int findFirstAlignment(FILE *fp, int format);

int findNextAlignment(FILE *fp, int fileFormat);

void freeAlignment();

void readAlignment(FILE *fp, int format, FILE * fpInfo, FILE * sfO);

void compressAlignment(alignment_struct *alignment);

void knuthShuffle(int orgArray[], int arraySize, int start, int end);

int parseSample(FILE *fp, int targetTok, char **word, int *wordLength, int curTok);

void unGetChars(FILE *fp, char* chars, int n);

int isValidDNACharacter(char input);

void parseIndividual(char *word, char** genotype, char delimiter, int token, char* states, char phased, char unphased, int *ploidy, int *max_ploidy);

void ignoreAll(FILE * fp, char * ent);

int isSpace(char ent);

double XchooseY_ln(int x, int y);

void checkSNIPPositions (FILE* fp, alignment_struct * alignment, int index);

void createSFS (FILE * fpSFS, FILE * fpSFSo, int alignmentIndex, int * SFSsize);

void createCLR (int grid);

void createPROBS (int probGrid);

void createPatterns(void);

void computeFactLUT(void);

double gettime(void);

t_sfs likelihoodSFS_SNP(t_sfs * tmpSFS, int x, int n, unsigned char f, int c);

void updateSF(FILE * fpSFo);

#ifdef _USE_PTHREADS

typedef struct

{

	int threadID;

	int threadTOTAL;

	int threadBARRIER;

	int threadOPERATION;

}threadData_t;

pthread_t * workerThreadL;

threadData_t * threadDataL;

int gridP;

void likelihoodSFS_thread(int tid, int threads);

void parallelComputationRVLUT_thread(int tid, int threads);

void computeAlpha_parallel(int grid);

void parallelComputationAlpha_thread(int tid, int threads);

void syncThreadsBARRIER();

#endif

#ifdef _DO_CHECKPOINTS

double lastCheckpointTime;

double checkPointInterval;

extern void writeCheckpoint();

#endif

#endif