msvar1.3.c

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

				l1 = logmut + pmval + fx	; 
				/* pmval chance of mut in that direction, 1/nl choose that
				      						lineage */
				lik +=  l1;
			}
			else if(pt->type == 2){
				l1 = logr +fx;
				lik += l1;
				--nlcheck;
			}
		}
		else if(t1 >= tf && t2 >= tf){
			if(pt->type == 1){
				
				mdir = pt->dl->val - pt->val;
				if(mdir == 1)pmval = logplus;
				else if(mdir == -1)pmval = logminus;
				else printerr("lik_cal: mdir bigger than 1 step");
				l1 = -0.5*nl*(t2 - t1)*(theta+r*(nl-1))+ pmval + logmut;
				lik +=  l1;
			}
			else if(pt->type == 2){
				l1 = -0.5*nl*(t2 - t1)*(theta+r*(nl-1))+
							logr; 
				lik += l1;
				--nlcheck;
			}
		}
			
		pt = pt->tlf;
		if(pt == stp)break; 
	}
	tf = 1.0;
	return lik;
} 



double cosim(double t1,int ni,double rval,double tf,int which_model)
{
	switch(which_model){
		case 0:
			return cosimL(t1,ni,rval,tf);
			break;
		case 1:
			return cosimE(t1,ni,rval,tf);
			break;
		default:
			printerr("cosim: which_model wrong");
	}
			
}

double cosimL(double t1,int ni,double rval,double tf)
{
	double rtf,cutoff,t2,deltat,dtop,tdash;
	
	dtop = (ni*(ni-1)*0.5);
	tdash = -log(gfsr8())/dtop;
	tdash += t1;
	
	
	rtf = rval*tf;
	deltat = tdash - t1;
	if(fabs(rval-1.0) < 0.00001 || tf > 1.0e10){ /* to avoid division by zero */
		return deltat + t1;
	}
	
	t2 = (exp(deltat*(1-rval)/rtf)*(rtf+t1-rval*t1) - rtf)/(1.0 - rval);
	if(t1 < tf && t2 < tf)return t2;
	else if(t1 < tf && t2 > tf){
		cutoff = rtf/(1-rval)*log(tf/(rtf+t1*(1.0-rval)));
		t2 = (deltat - cutoff)/rval + tf;
	}
	else t2 = deltat/rval + t1;
	return t2;
}



double cosimE(double t1,int ni,double rval,double tf)
{
	double rtf,cutoff,t2,deltat,dtop,tdash,logr;
	
	dtop = (ni*(ni-1)*0.5);
	tdash = -log(gfsr8())/dtop;
	tdash += t1;
	
	
	if(fabs(1.0-rval) < 1.0e-5 || tf > 1.0e10){ /*to avoid division by 0 */
		return tdash;
	}
	deltat = tdash - t1;
	
	if(t1 < tf){
		logr = log(rval);
		cutoff = (rval - pow(rval,t1/tf))*tf/logr;
		/* this is the value of (tdash - t1) corresponding to a real time of
	   		tf - when t2 == tf */
		if(deltat < cutoff){
			t2 = log(deltat*logr/tf+pow(rval,t1/tf))*tf/logr;

			return t2;
		}
		else {
			t2 = (deltat - cutoff)/rval + tf;
			return t2;
		}
	}
	else {
		t2 = deltat/rval + t1;
		return t2;
	}
}



void read_data(int ***freq,int ***val, int **noall,int *nloc,int **sums)
{
	FILE *inp;	
	int j,k,l,ic,imin,i;
	char c;
	
	inp = fopen("infile","r");
	if(inp == 0){
		printf("no infile\n");
		exit(1);
	}
	fscanf(inp,"%d",nloc);
	while(!((c=getc(inp)) == '\n' || c == '\f' || c == '\r'));
	
	(*noall) = (int *)malloc((*nloc)*sizeof(int));
	(*sums) = (int *)malloc((*nloc)*sizeof(int ));
	(*freq) = (int **)malloc((*nloc)*sizeof(int *));
	(*val) = (int **)malloc((*nloc)*sizeof(int *));
	for(j=0;j<(*nloc);++j){
		fscanf(inp,"%d",&(*noall)[j]);
		(*freq)[j] = (int *)malloc((*noall)[j]*sizeof(int ));
		(*val)[j] = (int *)malloc((*noall)[j]*sizeof(int ));
		while(!((c=getc(inp)) == '\n' || c == '\f' || c == '\r'));
		for(k=0;k<(*noall)[j];++k){
			ic = fscanf(inp,"%d",&(*freq)[j][k]);
			if(ic <= 0 || ic == EOF){
				printerr("error reading data");
			}
		}
		for(k=0;k<(*noall)[j];++k){
			ic = fscanf(inp,"%d",&(*val)[j][k]);
			if(ic <= 0 || ic == EOF){
				printerr("error reading data");
			}
		}
		imin = 10000000;
		for(i=0;i<(*noall)[j];++i){
			if((*val)[j][i] < imin)imin = (*val)[j][i];
		}
		for(i=0;i<(*noall)[j];++i){
			(*val)[j][i] -= imin;
		}
		(*sums)[j] = 0;
		for(i=0;i<(*noall)[j];++i){
			(*sums)[j] += (*freq)[j][i];
		}
	}
	fclose(inp);
	return;
}



	

void fill_tree_old(Node **list,int *freq,int *val,int noall,int sno,int *nn,int
				*nc,Cnode *clist,Genpars *pars)
{
	Node *lastp,*p1,*pt,*pn,*p2,*pp;
	Cnode *cp;
	int i,k,cumt,up,numlin,mutval,ic,nplus,nminus,j,*tvec,itemp,it1,it2,it3;
	double tt,tim1;
	int icount;
	double rr,rcheck;
	
	tvec = (int *)malloc(sno*sizeof(int));
	for(i=0,k=0,cumt=freq[0];i<sno;++i){
		while(i>=cumt){
			++k;
			if(k>=noall)printerr("fill_tree: k>=noall");
			cumt += freq[k];
		}
		tvec[i] = val[k];
	}
	for(i=0;i<sno-1;++i){
		itemp = tvec[i];
		ic = disrand(i+1,sno-1);
		tvec[i] = tvec[ic];
		tvec[ic] = itemp;
	}
	for(i=0;i<sno;++i){
		list[i]->val = tvec[i];
		list[i]->type = 0;
		list[i]->time = 0.0;
		list[i]->time2 = 0.0;
		list[i]->nlin = 0;
		list[i]->tlr = NULL;
		list[i]->tlf = NULL;
		list[i]->a = NULL;
		list[i]->dl = NULL;
		list[i]->dr = NULL;
		list[i]->pos = i;
		list[i]->mark = -1;
		list[i]->lno = i;
		list[i]->cp = NULL;
	}
	free(tvec);
	*nn = sno;
	tt = 0.0;
	list[0]->nlin = sno;
	lastp = list[0];
	*nc = 0;
	pars->linmut_indic = 0;
	while(*nc < sno-1){
		numlin = sno - *nc;
		icount = 0;
		while(1){
			it3 = it1 = disrand(0,numlin-1);
			for(i=0;i<*nn;++i){
				if(list[i]->a == NULL){
					if(it1 == 0)break;
					--it1;
				}
			}
			if(i == *nn)printerr("fill_tree: fallen through loop 1");
			p1 = list[i];
			it2 = disrand(0,numlin-2);
			if(it2 == it3)it2 = numlin-1;
			for(i=0;i<*nn;++i){
				if(list[i]->a == NULL){
					if(it2 == 0)break;
					--it2;
				}
			}
			if(i == *nn)printerr("fill_tree: fallen through loop 1");
			p2 = list[i];
			if(p1 == p2)printerr("fill_tree: p1 == p2");
			rr = pow(0.1,sqrt((double)abs(p1->val-p2->val)));
			rcheck = gfsr4();
			if(rcheck <= rr || *nc >= sno-2)break;
			if(icount > 1000000)printerr("fill_tree: icount > 1000000");
			++icount;
		}
		while(p1->val != p2->val){
			it1 = disrand(0,1);
			if(it1){
				pn = p1;
				up = p1->val < p2->val;
			}
			else{
				pn = p2;
				up = p2->val < p1->val;
			}
			tt += expdev()/(numlin*((numlin-1)+pars->theta)*0.5);
			if(up)mutval = pn->val+1;
			else mutval = pn->val-1;
			pn = mutate(list,nn,pn,tt,mutval,pn->lno);
			pn->dl->a = pn;
			pn->nlin = numlin;
			pn->tlr = lastp;
			lastp->tlf = pn;
			lastp = pn;
			if(it1)p1 = pn;
			else p2 = pn;
		}
		
		tt += expdev()/(numlin*((numlin-1)+pars->theta)*0.5);
		p1 = coalesce(list,nn,p1,p2,tt);
		p1->lno = sno + *nc;
		p1->nlin = numlin;
		p1->tlr = lastp;
		lastp->tlf = p1;
		lastp = p1;
		++(*nc);
		
/*--------------------------------------------------------
This bit added to load mlin, cumdpmlist, and diffpairtot 
--------------------------------------------------------*/
		/* left branch */
		ic = (*nc)-1;/* ic is the clist entry */
		clist[ic].mark[0] = -1;
		clist[ic].mark[1] = -1;
		nplus = 0;
		nminus = 0;
		pt = p1->dl;
		while(pt->type == 1){
			if(pt->val-pt->dl->val < 0)++nminus;
			else ++nplus;
			pt = pt->dl;
		}
		clist[ic].mlin[0][0] = nminus; 
		 	/* NB: minus (-1) refers back in time : ie the ancestor is 
				shorter; the descendent is longer */
		clist[ic].mlin[0][1] = nplus;
		clist[ic].mlin[0][2] = nminus + nplus;
		clist[ic].p = p1;
		p1->cp = &(clist[ic]);
		clist[ic].cpos = ic;
		clist[ic].l = pt;
		if(pt->type == 2){/*pt should be either type 0 or type 2 */
			pt->cp->u = p1;  
		}
		
		pt = p1->tlr;
		while(!(pt->type == 0  || pt->type == 2))pt=pt->tlr;
		clist[ic].tr = pt;
		if(pt ->type != 0)pt->cp->tf = p1;
		clist[ic].linmut[0] = (nminus*nplus > 0);
		pars->linmut_indic += clist[ic].linmut[0];
		
		
		/* right branch */
		nplus = 0;
		nminus = 0;
		pt = p1->dr;
		while(pt->type == 1){
			if(pt->val-pt->dl->val < 0)++nminus;
			else ++nplus;
			pt = pt->dl;
		}
		clist[ic].r = pt;
		if(pt->type == 2){
			pt->cp->u = p1;
		}
		clist[ic].mlin[1][0] = nminus;
		clist[ic].mlin[1][1] = nplus;
		clist[ic].mlin[1][2] = nplus + nminus;
		clist[ic].linmut[1] = (nminus*nplus > 0);
		pars->linmut_indic += clist[ic].linmut[1];
		
	}
	lastp->tlf = NULL;
	if(*nc != sno-1)printerr("fill_tree: nc != sno-1");
	clist[*nc-1].u = NULL;
	clist[*nc-1].tf = NULL;
	/* this bit fills in the upper lineage for mlin, 
			and gives node3_tot and cumnode3list */
	pars->nodemut_indic = 0;
	pars->linswap_indic = 0;
	pars->evswap_indic = 0;
	for(j=0;j<*nn;++j){
		if(linswap_cal(list[j])){
			list[j]->linswap = 1;
			++pars->linswap_indic;
		}
		else list[j]->linswap = 0;
		list[j]->evswap = swapable(list[j]);
		pars->evswap_indic += list[j]->evswap;
	}
	
	pars->lwt_tot = 0.0;
	pars->nwt_tot = 0.0;
	pars->lwt_max = -1.0;
	pars->nwt_max = -1.0;
	
	for(j=0;j<*nc;++j){
		cp = &(clist[j]);
		pp = cp->p;

		if(pars->lwt_max < (tim1 = lwt_cal(pp,0)))
			pars->lwt_max = tim1;
		cp->lwt[0] = tim1;
		pars->lwt_tot += tim1;
		if(pars->lwt_max < (tim1 = lwt_cal(pp,1)))
			pars->lwt_max = tim1;
		cp->lwt[1] = tim1;
		pars->lwt_tot += tim1;
		
		if(pars->nwt_max < (tim1 = nwt_cal(pp)))
			pars->nwt_max = tim1;
		pars->nwt_tot += tim1;
		cp->nwt = tim1;
		
		pt = cp->p->a;
		if(pt==NULL){
			cp->mlin[2][0] = 0;
			cp->mlin[2][1] = 0;
			cp->mlin[2][2] = 0;
			cp->nodemut = (cp->mlin[0][0]*cp->mlin[1][0] +
								cp->mlin[0][1]*cp->mlin[1][1] > 0);
			pars->nodemut_indic += cp->nodemut;
		}
		else{
			nplus = nminus = 0;
			while(pt->type == 1){
				if(pt->val-pt->dl->val < 0)++nminus;
				else ++nplus;
				pt = pt->a;
			}
			cp->mlin[2][0] = nminus;
			cp->mlin[2][1] = nplus;
			cp->mlin[2][2] = nminus+nplus;
			cp->nodemut =  (
								cp->mlin[0][0]*cp->mlin[1][0]*cp->mlin[2][1] +
								cp->mlin[0][1]*cp->mlin[1][1]*cp->mlin[2][0] > 0);
			pars->nodemut_indic += cp->nodemut;
		}
		
	}

	pars->pspace = MUTLIN_ADD_P + MUTNODE_ADD_P;
	if(pars->evswap_indic)pars->pspace += ESWAP_PROB;
	if(pars->linmut_indic)pars->pspace += MUTLIN_DEL_P;
	if(pars->nodemut_indic)pars->pspace += MUTNODE_DEL_P;
	if(pars->linswap_indic)pars->pspace += LSWAP_PROB;
}

void fill_tree(Node **list,int *freq,int *val,int noall,int sno,int *nn,int
				*nc,Cnode *clist,Genpars *pars)
{
	Node *lastp,*p1,*pt,*pn,*p2,*pp;
	Cnode *cp;
	int i,k,cumt,up,numlin,mutval,ic,nplus,nminus,j,*tvec,itemp,it1,it2,it3;
	double tt,tim1;
	int *ix,event_type,mut_target;
	double rr,rcheck,event_time;
	
	tvec = (int *)malloc(sno*sizeof(int));
	ix = (int *)malloc(2000*sizeof(int));
	mut_target = 0;
	for(i=0,k=0,cumt=freq[0];i<sno;++i){
		while(i>=cumt){ /* cumt for this only*/
			++k;
			if(k>=noall)printerr("fill_tree: k>=noall");
			cumt += freq[k];
		}
		tvec[i] = val[k]; /* this fills individual vector from freq. vectors */
		mut_target += val[k];
	}
	mut_target = (float)mut_target/sno + 0.5;
	
	for(i=0;i<sno-1;++i){ /* creates a random permutation - prob unnecessary */
		itemp = tvec[i];
		ic = disrand(i+1,sno-1);
		tvec[i] = tvec[ic];
		tvec[ic] = itemp;
	}
	for(i=0;i<sno;++i){
		list[i]->val = tvec[i];
		list[i]->type = 0;
		list[i]->time = 0.0;
		list[i]->time2 = 0.0;
		list[i]->nlin = 0;
		list[i]->tlr = NULL;
		list[i]->tlf = NULL;
		list[i]->a = NULL;
		list[i]->dl = NULL;
		list[i]->dr = NULL;
		list[i]->pos = i;
		list[i]->mark = -1;
		list[i]->lno = i;
		list[i]->cp = NULL;
	}
	free(tvec);
	
	*nn = sno;
	tt = 0.0;
	list[0]->nlin = sno;
	lastp = list[0];
	*nc = 0;
	pars->linmut_indic = 0;
	event_time = 0.0;
	while(*nc < sno-1){
		if(*nn > 1995)printerr("fill_tree: