-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathAFS.h
147 lines (105 loc) · 5.32 KB
/
AFS.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#ifndef AFS_H
#define AFS_H
#include <gsl/gsl_matrix.h>
#include <glib.h>
//afsObject represents a two population coalescent
typedef struct afsObject
{
int npops,nalleles;
int aCounts[2]; //this keeps track of popn specific allele counts
gsl_matrix_int **popMats;
char *matString;
}
afsObject;
typedef struct afsStateSpace
{
int nstates;
afsObject **states;
//below are only to hasten the process of state
//space enumeration
//int **nstateArray;
//afsObject ****tmpStorage;
GHashTable *stateHash;
}
afsStateSpace;
//for use in programs like cmc_island_mle
struct island_lik_params{
int snpNumber, maxSampleSize,nstates,n1,n2, hStates,nParams;
gsl_vector *weights, *sampleSizeVector,*rates;
//paramVector for HMM == [N2,thetaState1,thetaState2,...,thetaStaten, m1State1,m2State1,m1State2,m2State2,...,m1Staten,m2Staten] <n is number of hmm states>;
gsl_vector *paramVector,*paramCIUpper,*paramCILower, *mlParams;
//For straight MLE the paramVector takes the form [N1,N2,m1,m2]
gsl_matrix *log_afs, *posts, *transMat,*topol,*expAFS, *visitMat;
gsl_matrix_int *moveType;
afsStateSpace *stateSpace;
struct site *data;
double *cs_work;
int *dim1, *dim2; //these pointers ease the construction of sparse matrices
int nzCount;
double maxLik;
int optimized;
}island_lik_params;
#define MAXSTATES 100000000
#define MAXSIZE 20
afsObject *afsObjectNew(int n1, int n2);
void afsObjectFree(afsObject *v);
void afsObjectInit(afsObject *v,int n1, int n2);
void afsObjectPrint(afsObject *v);
afsObject *afsObjectNewFrom(afsObject *v);
afsObject *afsObjectDelta(afsObject *a, afsObject *b);
void afsObjectDeltaPre(afsObject *a, afsObject *b,afsObject *tmp);
int afsObjectCount(afsObject *a);
afsStateSpace *afsStateSpaceNew();
void afsStateSpaceFree(afsStateSpace *a);
int afsStateSpaceContainsState(afsStateSpace *S, afsObject *aState);
void afsStateSpaceRemoveAbsorbing(afsStateSpace *S);
int afsStateSpaceContainsStateTmpStorage(afsStateSpace *S, afsObject *aState);
void afsStateSpaceMapPopn(afsStateSpace *S, int *map);
void afsStateSpaceMapAndReducePopn(afsStateSpace *S, int *map, afsStateSpace *reducedSpace, int *reverseMap);
void afsStateSpaceRemovePopulation(afsStateSpace *S, int popn);
afsStateSpace *afsStateSpaceCopy(afsStateSpace *S);
void mcMatsImportFromFile(const char *fileName,int *nnz,double *topol, int *moveA, int *dim1, int *dim2);
int nDescPop1(afsObject *v);
int nDescPop2(afsObject *v);
int compatibleAfsGivenInitial(afsObject *v, int n1, int n2);
int afsObjectsEqual(afsObject *a, afsObject *b);
void acceptOrRejectProposedState(afsStateSpace *S, afsObject *child, int n1, int n2);
int afsObjectQsortCompare(const void *p1, const void *p2);
afsStateSpace *afsStateSpaceImportFromFile(const char *fileName);
void nonZeroEntries(gsl_matrix_int *m, int **coords, int *count);
void entriesGreaterThan(gsl_matrix_int *m, int **coords, int *count, int cut);
void entriesLessThan(gsl_matrix_int *m, int **coords, int *count, int cut);
int countNegativeEntries(gsl_matrix_int *m);
double xchoosey(int n, int k);
int matrixAllGTZ(gsl_matrix_int *m);
int matrixSum(gsl_matrix_int *m);
void gsl_matrix_prettyPrint(gsl_matrix *m);
double matrixSumDouble(gsl_matrix *m);
void preorderEvolveState(afsStateSpace *S, afsObject *aState, int n1, int n2);
void coalMarkovChainTopologyMatrix(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int n1, int n2);
int coalMarkovChainTopologyMatrix_2(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int *dim1, int *dim2);
int coalMarkovChainTopologyMatrix_sparse(afsStateSpace *S,double *topol, int *moveType, int *dim1, int *dim2);
void readTopolFile(const char *fileName, double *topA, int *moveA, int *dim1, int *dim2);
void fillTransitionMatrix(afsStateSpace *S, gsl_matrix *topol, gsl_matrix_int *moveType, gsl_matrix *transMat, gsl_vector *rates,
double theta1, double theta2, double mig1, double mig2);
struct cs_di_sparse *fillTransitionMatrix_2(afsStateSpace *S, gsl_matrix *topol, gsl_matrix_int *moveType, gsl_vector *rates, int nzCount,
int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2);
void fillExpectedAFS(afsStateSpace *S, gsl_matrix *transMat,gsl_vector *rates, gsl_matrix *expAFS, int n1, int n2);
void fillExpectedAFS_visits(afsStateSpace *S, gsl_matrix *visitMat,gsl_vector *rates, gsl_matrix *expAFS);
void fillExpectedAFS_visits_2(afsStateSpace *S, double *visitMat,gsl_vector *rates, gsl_matrix *expAFS);
void overlayPopnSplitting(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int n1, int n2);
void maximize_island_lik_const(double *lik, void * p);
void maximize_island_lik_hmm_const(double *lik, void * p);
void maximize_island_lik_CILower(double *lik, void * p);
void maximize_island_lik_CIUpper(double *lik, void * p);
void island_lik_constraint_CI(size_t n, const double *pGuess, void * p, double *result);
void island_lik_constraint(size_t n, const double *pGuess, void * p, double *result);
void island_lik_constraint_hmm(size_t n, const double *pGuess, void * p, double *result);
void fillLogAFS_ptr_sp(void * p);
void fillLogAFS_ptr_sp2(void * p);
void fillLogAFS_ptr_sp_hmm(int state,void * p);
void itoa(int n, char s[]);
void reverse(char s[]);
void NumericalGradient_island_lik(size_t n, const double *v, void *params, double *df);
void island_lik_fdf(const size_t n, const double *x,void *fparams,double *fval,double *grad);
#endif