msvar1.3.c

MCMC（马尔可夫－盟特卡罗方法）实现的程序

/* This program performs a Hierarchical Bayesian analysis using multilocus
microsatellite data to infer population growth and decline


*/


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#ifdef __BORLANDC__
#include <float.h>
#endif
#include "myutil.h"
#include "myutil.c"

#define M_LOG10E        0.43429448190325182765

/* original set */
#define MUTLIN_ADD_P 0.1
#define MUTLIN_DEL_P 0.1
#define MUTNODE_ADD_P 0.2
#define MUTNODE_DEL_P 0.2
#define ESWAP_PROB  0.2 
#define LSWAP_PROB  0.2 
#define CUM_MUTLIN_DEL_P 0.2
#define CUM_MUTNODE_ADD_P 0.4
#define CUM_MUTNODE_DEL_P 0.6
#define CUM_ESWAP_PROB 0.8
#define CUM_LSWAP_PROB 1.0

#define MCONS 12

struct node{
	int val;
	int val2;
	int type;   /* 0=sample,1=mutation,2=coalescence,3=ancestor */
	double time;
	double time2;
	int nlin;
	int nlin2;
	int pos;
	struct cnode *cp;
	int mark;
	int lno;
	int lno2;
	int linswap;
	int linswap2;
	int evswap;
	int evswap2;
	struct node *tlr;
	struct node *tlr2;
	struct node *tlf;
	struct node *tlf2;
	struct node *a;
	struct node *a2;
	struct node *dl;
	struct node *dl2;
	struct node *dr;
	struct node *dr2;
};


typedef struct node Node;

struct cnode{
	Node *p;
	Node *u;
	Node *u2;
	Node *l;
	Node *l2;
	Node *r;
	Node *r2;
	Node *tf;
	Node *tf2;
	Node *tr;
	Node *tr2;
	int mlin[3][3];
	int mlin2[3][3];
	int linmut[2];
	int linmut2[2];
	int nodemut;
	int nodemut2;
	double nwt;
	double nwt2;
	double lwt[2];
	double lwt2[2];
	int cpos;
	int mark[2];
};
	
typedef struct cnode Cnode;

struct hierpars{
	double ln0mean;
	double ln0var;
	double ln1mean;
	double ln1var;
	double lmumean;
	double lmuvar;
	double ltfamean;
	double ltfavar;

	double ln0mean2;
	double ln0var2;
	double ln1mean2;
	double ln1var2;
	double lmumean2;
	double lmuvar2;
	double ltfamean2;
	double ltfavar2;
	
	double ln0mm;
	double ln0mv;
	
	double ln0vm;
	double ln0vv;
	
	double ln1mm;
	double ln1mv;
	
	double ln1vm;
	double ln1vv;
	
	double lmumm;
	double lmumv;
	
	double lmuvm;
	double lmuvv;
	
	double ltfamm;
	double ltfamv;
	
	double ltfavm;
	double ltfavv;
	
	double prior;
	double prior2;


};

typedef struct hierpars Hierpars;

struct genpars{
	int evswap_indic;
	int evswap_indic2;
	int linmut_indic;
	int linmut_indic2;
	int nodemut_indic;
	int nodemut_indic2;
	int linswap_indic;
	int linswap_indic2;
	double nwt_tot;
	double nwt_tot2;
	double lwt_tot;
	double lwt_tot2;
	double nwt_max;
	double nwt_max2;
	double lwt_max;
	double lwt_max2;
	double pspace;
	double pspace2;
	double tpf;
	double tpr;
	double lik;
	double lik2;
	
	double theta;
	double theta2;
	double r;
	double r2;
	double tf;
	double tf2;
	
	
	int which_model;
	
	
	double n0;
	double n0_2;
	double n1;
	double n1_2;
	double mu;
	double mu2;
	double gen;
	double gen2;
	double tfa;
	double tfa2;
	
	Node *base;
	double prior;
	double prior2;
	double pplus;
	double pplus2;
	int sno ;
	
};

typedef struct genpars Genpars;

struct Gplist{
	int keeptime;
	int keepn0;
	int keepn1;
	int keepmu;
	int keeptfa;
	int keeppplus;
	int scale_r;
	int scale_tf;
	int tfuni;
	double dd;
	double ppval;
};

void twidpars_var(struct Gplist *ptlist,Genpars *pars,Hierpars *hpars);
void untwidhpars(Hierpars *hpars);
double twidhpars(struct Gplist *ptlist,Hierpars *hpars);
void copyhpars(Hierpars *hpars);
void inithpars(Hierpars *hpars);
void treesummary(int nn,Genpars *pars,Node **list,double *above, double
*mutbc, double *blpc);
void choosepar(struct Gplist *ptlist,Hierpars *hpars,int oneonly,
int *mtype) ;
void copypars(Genpars *pars);
void priorcalc(Genpars *pars,Hierpars *hpars);
void hpriorcalc(Hierpars *hpars);
void twidpars(struct Gplist *ptlist,Genpars *pars,Hierpars *hpars,int doprior);
void untwidpars(Node **mlist,int *mm,Cnode **cmlist,int *cmm,Genpars *pars);
double boundcheck(double l1,int nn,Genpars *pars,Hierpars *hpars);
void initpars(Genpars *pars,int j);
void fillpars(Genpars *pars);
double lwt_cal(Node *pp,int irl);
double nwt_cal(Node *pp);
double rejprob(Node *pp,int nc);
double tdens(double t1,double t2,int nl,double theta,double r,double tf,int which_model);
double tdensL(double t1,double t2,int nl,double theta,double r,double tf);
double tdensE(double t1,double t2,int nl,double theta,double r,double tf);
double cosim(double t1,int ni,double rval,double tf,int which_model);
double cosimL(double t1,int ni,double rval,double tf);
double cosimE(double t1,int ni,double rval,double tf);
void assimlin(Cnode *cpt,int vec[],int lr);
void linswap(Node **list,int nn,Node **mlist,int *mm,Cnode **cmlist,int *cmm,Genpars *pars);
double lik_cal(Node *pt,Node *out,Genpars *pars);
double lik_calL(Node *pt,Node *out,Genpars *pars);
double lik_calE(Node *pt,Node *out,Genpars *pars);
void proptchange(struct Gplist *ptlist,Genpars *pars,Cnode *clist,Node **mlist,int *mm,
	Cnode **cmlist,int *cmm);
void tchange(Genpars *pars,Cnode *clist,Node **mlist,int *mm,Cnode **cmlist,int *cmm);
int swapable(Node *p);
void eventswap(Node **list,int nn,Node **mlist,int *mm,Cnode **cmlist,int *cmm,
	Genpars *pars);
int cal_nodeprob(int mlin[3][3],int indic);
int cal_linprob(int mlin[3][3],int rl);
int findinclist(Node *pp);
void checktree(Node *pt,int ulim,int sno,Cnode *clist,int nc,Genpars *pars,int ismark);
void read_data(int ***freq,int ***val, int **noall,int *nloc,int **sums);
void fill_tree(Node **list,int *freq,int *val,int noall,int sno,int *nn,int
	*nc,Cnode *clist,Genpars *pars);
void fill_tree_old(Node **list,int *freq,int *val,int noall,int sno,int *nn,int
	*nc,Cnode *clist,Genpars *pars);
void getevent(double t0,int nlin, int *event_type,double *event_time,Genpars *pars);
Node *mutate(Node **list,int *nn,Node *n1,double tt,int mutval,int lno);
Node *coalesce(Node **list,int *nn,Node *n1,Node *n2,double tt);
void twiddle(Node **list,int *oldnn, int *nn, Cnode *clist,int nc,Node **mlist,
	int *mm,Cnode **cmlist,int *cmm,int *mtype,
	Genpars *pars);
void mutlin_add2(Node **list,int *nn, Cnode *clist,int nc,Node **mlist,int *mm,
	Cnode **cmlist,int *cmm,Genpars *pars);
void mutlin_del2(Node **list,int *nn, Cnode *clist,int nc,int *ifail,Node **mlist,
	int *mm,Cnode **cmlist,int *cmm,
	Genpars *pars);
void mutnode_add3(Node **list,int *nn, Cnode *clist,int nc,Node **mlist,int *mm,
	Cnode **cmlist,int *cmm,int *mtype,Genpars *pars);
void mutnode_del3(Node **list,int *nn, Cnode *clist,int nc,int *ifail,Node **mlist,
	int *mm,Cnode **cmlist,int *cmm,int *mtype,
	Genpars *pars);
Node *insertmut(Node **list,int *nn,double mt1,Node *pt,int mutval,
	Node **mlist,int *mm,int lno);
void deletemut(Node *pp,Node **list,int *nn,Node **mlist,int *mm);
Node *cfind_down(Node *pt,int indic);
Node *cfind_up(Node *pt,int *rl);
void markup(Node *pt,char *str,Node *target,Node **mlist,int *mm);
void cmarkup(Cnode *pt,char *str,Cnode *target,Cnode **cmlist,int *cmm);
int make_mark(int mark,int val);
int inside(int ic,int mark);
void untwiddle(Node **mlist,int *mm,Cnode **cmlist,int *cmm,Genpars *pars,int *nn,int oldnn);
void restore(Node *pp);
void crestore(Cnode *pp);
int decomp(int *mark);
void resit(Node *pp,int ic);
void cresit(Cnode *pp,int ic,int ip);
void accept(Node **mlist,int *mm,Cnode **cmlist,int *cmm);
Node *getmnode(Node *pp,char *ss,int val);
double log_normdev(double x,double mu,double s);
int linswap_cal(Node *p1);
double normdev(double x,double mu,double s);
double lndev(double x,double mu,double s);
void read_params();


int CHeck=0;

int Mod_type;
int Itno;
double Init_n0[100],Init_n1[100],
		Init_mu[100],Init_gen,Init_tfa[100];
double Init_ln0mean,Init_ln0var,Init_ln1mean,Init_ln1var,
       Init_lmumean,Init_lmuvar,Init_ltfamean,Init_ltfavar,
       Init_ln0mm,Init_ln0mv,Init_ln0vm,Init_ln0vv,
       Init_ln1mm,Init_ln1mv,Init_ln1vm,Init_ln1vv,
       Init_lmumm,Init_lmumv,Init_lmuvm,Init_lmuvv,
       Init_ltfamm,Init_ltfamv,Init_ltfavm,Init_ltfavv;
       
int Thin;
int Nloc;
double Usum[10];


Node **Block;
FILE *checkout;

#ifdef __MWERKS__
int MTspeed = 500;
float MaxMTspeed = 10000.0;
float MinMTspeed = 25.0;
#endif
int Illegal;
int Ploidy;
int PredictDiffpair_tot,PredictNode3_tot,PredictCswaptot,PredictEswaptot;
double PredictPspace;
int main(void)
{
	Node **list[50],**mlist[50];
	Cnode *clist[50],**cmlist[50];
	Genpars *pars[50],realpars[50];
	Hierpars *hpars,realhpars;
	int j,k,i,**freq,**val,*noall,*sno,nloc,nn[50],nc[50],iter,mm[50],cmm[50],
		oldnn[50],mtype,mtlist[30],mplist[30],mtypecount,
		twidtheta;
	double l1,likold,liknew,ldiff,ldiff2,tdiff,pdiff,mutbc,above,blpc;
	double Scale2=100.0;
	int rcheck;
	double lcheck1,lsum;
	int iloc,counter;
	FILE *lfile[50],*mpout,*hfile;
	char fstring[50][10],ncharstr[100];
	struct Gplist ptlist;
	double hsucc=0.0,htot=0.0;
	
	int oneonly,bl,ul;
	
	#ifdef __BORLANDC__

	_control87(MCW_EM,MCW_EM);

	#endif



	opengfsr();
	
	read_data(&freq,&val,&noall,&nloc,&sno);
	Nloc = nloc;
	read_params();
	
	for(j=0;j<nloc;++j){
		sprintf(ncharstr,"%d",j+1);
		strcpy(fstring[j],"out");
		lfile[j] = fopen(strcat(fstring[j],ncharstr),"w");
	}
	hfile = fopen("hpars.dat","w");

	Block = (Node **)malloc(nloc*sizeof(Node *));
	for(j=0;j<nloc;++j)Block[j] = (Node *)malloc(2000*sizeof(Node));



	for(j=0;j<nloc;++j)pars[j] = &realpars[j];
	hpars = &realhpars;
	
	
	for(j=0;j<nloc;++j)list[j] = (Node **)malloc(2000*sizeof(Node *));
	for(j=0;j<nloc;++j)mlist[j] = (Node **)malloc(2000*sizeof(Node *));
	
	for(j=0;j<30;++j){mtlist[j] = 0;mplist[j] = 0;}
	mtypecount = 0;
	
	
	
	for(k=0;k<nloc;++k){
		for(j=0;j<2000;++j){
			list[k][j] = &Block[k][j];
		}
	}
	
	for(j=0;j<nloc;++j){pars[j]->sno = sno[j];}
	for(j=0;j<nloc;++j)clist[j] = (Cnode *)malloc((sno[j]-1)*sizeof(Cnode));
	for(j=0;j<nloc;++j)cmlist[j] = (Cnode **)malloc((sno[j]-1)*sizeof(Cnode *));
	
	
	for(j=0;j<nloc;++j)initpars(pars[j],j);
	
	inithpars(hpars);
	
	for(j=0;j<nloc;++j)fill_tree(list[j],freq[j],val[j],
		noall[j],sno[j],&nn[j],&nc[j],clist[j],pars[j]);
	for(j=0;j<nloc;++j)checktree(list[j][0],nn[j],sno[j],clist[j],nc[j],pars[j],0);
		
	for(j=0;j<nloc;++j)pars[j]->base = list[j][0];
	for(j=0;j<nloc;++j)priorcalc(pars[j],hpars);
	hpriorcalc(hpars);
	hpars->prior2 = hpars->prior;
	for(j=0;j<nloc;++j)pars[j]->prior2 = pars[j]->prior;
	likold = 0.0;
	for(j=0;j<nloc;++j)likold += lik_cal(list[j][0]->tlf,NULL,pars[j]);
	
	for(j=0;j<nloc;++j)cmm[j] = mm[j] = 0;
	counter = 0;
	
	#ifdef __MWERKS__
	firstime = clock()/(3600.0*CLOCKS_PER_SEC);	
	itcounter = 100;
	firstiter = 0;
	itdiff = Itno/100;
	itcounter = itdiff;
	#endif
	printf("starting MCMC\n");
	for(iter=0;iter<Thin;++iter){
		Illegal = 0;
		if(counter == Itno)break;
	
		#ifdef __MWERKS__
		if((Rk++)%MTspeed==0)MyEventLoop(&MyTick);
		if((iter+1) == itcounter){
			ruptime = clock()/(3600.0*CLOCKS_PER_SEC);
			printf("iter is %d      running for %f hours\n",iter+1,ruptime-firstime);
			togo = Itno-iter;
			idone = iter-firstiter;
			estogo = (ruptime-firstime)/idone*togo;
			printf("estimate %f hours to go\n",estogo);
			itcounter = itcounter+itdiff;
		}
		#endif
		
		if((iter+1)%500 == 0)rcheck = 1;
		else rcheck = 0;
		twidtheta = 0; /* change G */
		if(gfsr8() < 0.01)twidtheta = 1; /* change params for one locus
										without changing hierarchical params */
		else if(gfsr8() < 0.01)twidtheta = 2; /* change params for all loci and
		                                       change prior mean */
		else if(gfsr8() < 0.01)twidtheta = 3; /* Change params for all loci and 
										change prior variance*/
		iloc = disrand(0,nloc-1);
		if(twidtheta == 1){ /* change params for one locus
										without changing hierarchical params */
			choosepar(&ptlist,hpars,1,&mtype); /* This determines
						type of updates for mutation and demographic
						parameters. */
			pars[iloc]->lik2 = lik_cal(pars[iloc]->base->tlf,
				NULL,pars[iloc]);
			copypars(pars[iloc]);
			
			twidpars(&ptlist,pars[iloc],hpars,0);
			if(!ptlist.keeptime)
				tchange(pars[iloc],clist[iloc],mlist[iloc],&mm[iloc],
					cmlist[iloc],&cmm[iloc]);
			
			if((CHeck || rcheck)){
				checktree(list[iloc][0],nn[iloc],sno[iloc],clist[iloc],