msvar1.3.c

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

		ddscale = 1.0;
		if(oneonly)ddscale *= 0.5;
	/*	ptlist->dd = ddscale*norm4(); */
		rr2 = gfsr4();
		parwt = 0.5;
		if(Usum[5] <= 0.0)parwt = 0.0;
		if(rr2 < parwt){
			rr2 = gfsr4(); 
			if(rr2 < Usum[0]){
			/*	keepn0 = 0; */
			/*	if(oneonly){
					ddscale = hpars->ln1var < hpars->lmuvar ? 
					          hpars->ln1var : hpars->lmuvar ;
				}*/
				ptlist->dd = ddscale*norm4();
				ptlist->keepn1 = 0; 
				ptlist->keepmu = 0;
				if(gfsr4() < 0.5)ptlist->keeptime = 0;
				*mtype = 20;
			}
			else if(rr2 < Usum[1]){
			/*	if(oneonly){
					ddscale = hpars->ltfavar < hpars->lmuvar ?
								hpars->ltfavar : hpars->lmuvar;
				}*/
				ptlist->dd = ddscale*norm4();
				ptlist->keepmu = 0;
				ptlist->keeptfa = 0;
				if(gfsr4() < 0.5)ptlist->keeptime = 0;
				*mtype = 21;
			}
			else if (rr2 < Usum[2]){
			/*	if(oneonly){
					temp1 = hpars->ln0var < hpars->ln1var ?
							hpars->ln0var : hpars->ln1var;
					ddscale = hpars->ltfavar < hpars->lmuvar ?
							hpars->ltfavar : hpars->lmuvar;
					ddscale = ddscale < temp1 ? ddscale : temp1;
				}*/
				ptlist->dd = ddscale*norm4();
				ptlist->keepn0 = 0;
				ptlist->keepn1 = 0;
				ptlist->keeptfa = 0;
				ptlist->keepmu = 0;
				*mtype = 22;
			}
			else if (rr2 < Usum[3]){
			/*	if(oneonly)ddscale = hpars->ln1var ;*/
				ptlist->dd = ddscale*norm4();
				ptlist->keepn1 = 0;
				if(gfsr4() < 0.5)ptlist->keeptime = 0;
				*mtype = 23;
			}
			else if(rr2 < Usum[4]){
			/*	if(oneonly){
					ddscale = hpars->ln1var < hpars->lmuvar ?
							  hpars->ln1var : hpars->lmuvar;
					ddscale = ddscale < hpars->ltfavar ? ddscale :
								hpars->ltfavar;	
				} */
				ptlist->dd = ddscale*norm4();
			/* I think an equivalent change is to just change
			   n0 rather than the other 3, so I do that here
				ptlist->keepn1 = 0;
				ptlist->keepmu = 0;
				ptlist->keeptfa = 0; */
				ptlist->keepn0 = 0;
				if(gfsr4() < 0.5)ptlist->keeptime = 0;
				*mtype = 29;
			}
			else if (rr2 < Usum[5]){
			/*	if(oneonly)ddscale = hpars->lmuvar;*/
				ptlist->dd = ddscale*norm4();
				ptlist->keepmu = 0;
				ptlist->keeptime = 0;
				*mtype = 24;
			}
			else {
			/*	if(oneonly)ddscale = hpars->ltfavar;*/
				ptlist->dd = ddscale*norm4();
				if(!oneonly)ptlist->dd = 5*ptlist->dd; /* to try to leap around more */
				ptlist->keeptfa = 0;
				
				if(gfsr4() < 0.5)ptlist->keeptime = 0;
				*mtype = 25;
			}
		}
		else {
			ptlist->keeptime=0;
			*mtype = 26;
		}
	}
	else{
		ptlist->keeppplus = 0;
		ptlist->ppval = gfsr4();
		*mtype = 29;
	}
	

}

void twidpars(struct Gplist *ptlist,Genpars *pars,Hierpars *hpars,int doprior)
{
	double dd;
	static double ltrans=M_LOG10E;
	dd = ptlist->dd;
	
	pars->tpr = pars->tpf = 0.0; 
	if(!ptlist->keepn0){
		pars->n0 = pars->n0*exp(dd);
		if(doprior)hpars->ln0mean += ltrans*dd;	
	}
	if(!ptlist->keepn1){
		pars->n1 = pars->n1*exp(dd);
		if(doprior)hpars->ln1mean += ltrans*dd;
	}
	if(!ptlist->keepmu){
		pars->mu = pars->mu*exp(-dd);
		if(doprior)hpars->lmumean -= ltrans*dd;
	}
	if(!ptlist->keeptfa){
		pars->tfa = pars->tfa*exp(dd); 
		if(doprior)hpars->ltfamean += ltrans*dd;
	}
	if(!ptlist->keeppplus){
		pars->pplus = ptlist->ppval;
	}
	fillpars(pars);

}

void twidpars_var(struct Gplist *ptlist,Genpars *pars,Hierpars *hpars)
{
	double dd,xx;
	static double ltrans=M_LOG10E;
	dd = ptlist->dd;
	
	pars->tpr = pars->tpf = 0.0; 
	if(!ptlist->keepn0){
		xx = exp(dd)*(log10(pars->n0) - hpars->ln0mean) + hpars->ln0mean;
		pars->n0 = /*pars->n0*/pow(10.0,xx);
	}
	if(!ptlist->keepn1){
		xx = exp(dd)*(log10(pars->n1) - hpars->ln1mean) + hpars->ln1mean;
		pars->n1 = /*pars->n1*/pow(10.0,xx);
	}
	if(!ptlist->keepmu){
		xx = exp(dd)*(log10(pars->mu) - hpars->lmumean) + hpars->lmumean;
		pars->mu = /*pars->mu*/pow(10.0,xx);
	}
	if(!ptlist->keeptfa){
		xx = exp(dd)*(log10(pars->tfa) - hpars->ltfamean) + hpars->ltfamean;
		pars->tfa = /*pars->tfa*/pow(10.0,xx);
	}
	if(!ptlist->keeppplus){
		pars->pplus = ptlist->ppval;
	}
	fillpars(pars);
	pars->tpr = dd;
	
	/* question - do we need to do anything for tpr/tpf? */

}

double twidhpars(struct Gplist *ptlist,Hierpars *hpars)
{
	/* using gaussian updates means that the variances need to be checked 
	with boundcheck */
	double dd,rr;
	double scale=0.5;
	

	ptlist->keeptime 
		= ptlist->keepn0 
		= ptlist->keepn1 
		= ptlist->keepmu 
		= ptlist->keeptfa 
		= ptlist->keeppplus
		= 1;


	rr = gfsr4();
	dd = norm4()*scale;
	ptlist->dd = dd;
	if(rr < 0.25){
		hpars->ln0var *= exp(dd);
		ptlist->keepn0 = 0;
	}
	else if(rr < 0.5){
		hpars->ln1var *= exp(dd);
		ptlist->keepn1 = 0;
		
	}
	else if(rr < 0.75){
		hpars->lmuvar *= exp(dd);
		ptlist->keepmu = 0;
	}
	else {
		hpars->ltfavar *= exp(dd);
		ptlist->keeptfa = 0;
	}
	if(gfsr4() < 0.5)ptlist->keeptime = 0;
	
	return dd;

}

void priorcalc(Genpars *pars,Hierpars *hpars)
{
/*The idea behind this is that, say, for each locus the  mutation rate is different, but
it is drawn from some distribution common to all loci with  mean lmumean, and 
variance lmuvar. If we have a number of loci then we can also say something about the
distribution of lmumean and lmuvar. For example there may be little variation between loci 
and therefore lmuvar will be small. 
 
 */
	pars->prior = log_normdev(log10(pars->n0),hpars->ln0mean,hpars->ln0var);
	pars->prior += log_normdev(log10(pars->n1),hpars->ln1mean,hpars->ln1var);
	pars->prior += log_normdev(log10(pars->mu),hpars->lmumean,hpars->lmuvar);
	pars->prior += log_normdev(log10(pars->tfa),hpars->ltfamean,hpars->ltfavar);
	
}

void hpriorcalc(Hierpars *hpars)
{
/*The idea behind this is that, say, for each locus the  mutation rate is different, but
it is drawn from some distribution common to all loci with  mean lmumean, and 
variance lmuvar. If we have a number of loci then we can also say something about the
distribution of lmumean and lmuvar. For example there may be little variation between loci 
and therefore lmuvar will be small. 
 
 */
	hpars->prior = log_normdev(hpars->ln0mean,hpars->ln0mm,hpars->ln0mv);
	hpars->prior += log_normdev(hpars->ln0var,hpars->ln0vm,hpars->ln0vv);
	hpars->prior += log_normdev(hpars->ln1mean,hpars->ln1mm,hpars->ln1mv);
	hpars->prior += log_normdev(hpars->ln1var,hpars->ln1vm,hpars->ln1vv);
	hpars->prior += log_normdev(hpars->lmumean,hpars->lmumm,hpars->lmumv);
	hpars->prior += log_normdev(hpars->lmuvar,hpars->lmuvm,hpars->lmuvv);
	hpars->prior += log_normdev(hpars->ltfamean,hpars->ltfamm,hpars->ltfamv);
	hpars->prior += log_normdev(hpars->ltfavar,hpars->ltfavm,hpars->ltfavv);
}


void copypars(Genpars *pars)
{
	pars->prior2 = pars->prior;
	pars->n0_2 = pars->n0;
	pars->n1_2 = pars->n1;
	pars->mu2 = pars->mu;
	pars->gen2 = pars->gen;
	pars->tfa2 = pars->tfa;
	pars->pplus2 = pars->pplus;
	
	pars->lwt_tot2 = pars->lwt_tot;
	pars->nwt_tot2 = pars->nwt_tot;
	pars->lwt_max2 = pars->lwt_max;
	pars->nwt_max2 = pars->nwt_max;
	
	fillpars(pars);
}

void fillpars(Genpars *pars)
{
	pars->theta = Ploidy*pars->n0*2.0*pars->mu;
	pars->theta2 = Ploidy*pars->n0_2*2.0*pars->mu2;
	pars->r = pars->n0/pars->n1;
	pars->r2 = pars->n0_2/pars->n1_2;
	pars->tf = pars->tfa/(pars->gen*Ploidy*pars->n0);
	pars->tf2 = pars->tfa2/(pars->gen2*Ploidy*pars->n0_2);
}
		


void untwiddle(Node **mlist,int *mm,Cnode **cmlist,int *cmm,Genpars *pars,int *nn,int oldnn)
{
	int j;
	for(j=0;j<(*mm);++j){
		restore(mlist[j]);
	}
	(*mm) = 0;
	
	for(j=0;j<(*cmm);++j){
		crestore(cmlist[j]);
	}
	(*cmm) = 0;
	
	pars->evswap_indic = pars->evswap_indic2;
	pars->linmut_indic = pars->linmut_indic2;
	pars->nodemut_indic = pars->nodemut_indic2;
	pars->linswap_indic = pars->linswap_indic2;
	pars->nwt_tot = pars->nwt_tot2;
	pars->nwt_max = pars->nwt_max2;
	pars->lwt_tot = pars->lwt_tot2;
	pars->lwt_max = pars->lwt_max2;
	pars->pspace = pars->pspace2;
	
	*nn = oldnn;
}

void accept(Node **mlist,int *mm,Cnode **cmlist,int *cmm)
{
	int j;
	
	for(j=0;j<*mm;++j)mlist[j]->mark = -1;
	(*mm) = 0;
	for(j=0;j<*cmm;++j){cmlist[j]->mark[0] = -1;cmlist[j]->mark[1] = -1;}
	(*cmm) = 0;
	
}


int cal_nodeprob(int mlin[3][3],int indic)
{
	int prob;
	if(indic){
		prob = mlin[0][0]*mlin[1][0] + mlin[0][1]*mlin[1][1] > 0;
	}
	else{
		prob = mlin[0][0]*mlin[1][0]*mlin[2][1] + mlin[0][1]*mlin[1][1]*mlin[2][0] > 0;
	}
	return prob;
}

int cal_linprob(int mlin[3][3],int rl)
{
	int prob;
	prob = mlin[rl][0]*mlin[rl][1] > 0;
	return prob;
}

double nwt_cal(Node *pp)
{
	double tim1;
	
	if(pp == NULL || pp->type != 2)printerr("nwt_cal: wrong input");
/*	return 1.0; */
	if(pp->a == NULL)tim1 = (pp->time - pp->cp->l->time)*(pp->time -
	pp->cp->r->time);
	else tim1 = (pp->cp->u->time - pp->time)*
		(pp->time - pp->cp->l->time)*(pp->time - pp->cp->r->time);
	return tim1;

}

double lwt_cal(Node *pp,int irl)
{
	if(pp == NULL || pp->type != 2)printerr("lwt_cal: wrong input pp");
	if(!(irl == 0 || irl == 1))printerr("lwt_cal: wrong input irl");
/*	return 1.0; */
	if(irl == 0)return (pp->time - pp->cp->l->time)*(pp->time -
	pp->cp->l->time);
	else return (pp->time - pp->cp->r->time)*(pp->time - pp->cp->r->time);
	
}


int findinclist(Node *pp)
{
	Cnode *cp;
	if(pp == NULL)return -1;
	if(pp->type == 0)return -1;
	cp = pp->cp;
	if(cp == NULL)return -1;
	return cp->cpos;
}

void checktree(Node *pt,int ulim,int sno,Cnode *clist,int nc,Genpars *pars,int ismark)
{
	int nl,blim,ic,lastval1,lastval2,j,nminus,nplus,n3sum,pairsum,ll;
	int ncheck,  sumswn,  swn,ltot;
	double pscheck,lwtot,lwmax,nwtot,nwmax,tim1;
	Node *pp,*pv,*startpp;
	Cnode *cp,*cpt,*cpp;
	blim = pt->pos;
	ulim = ulim-1;
	ic = 0;
	if(pt->type != 0){
		printerr("treecheck: error 000");
	}
	else{
		if(blim == 0 && pt->pos > sno)printerr("treecheck: error 00");
		nl = pt->nlin;
		pt = pt->tlf;
		startpp = pt;
	}
	ltot = 0;
	sumswn = 0;
	while(pt != NULL){
		if(!ismark && pt->mark != -1)printerr("treecheck: error a00");
		if(pt->type == 0)printerr("treecheck: error 00");
		if(pt->pos < blim || pt->pos > ulim)printerr("treecheck: error 0");
		if(pt->a != NULL && (pt->a->pos < blim || pt->pos > ulim))
			printerr("treecheck: error 0a");
		if(pt->nlin != nl)printerr("treecheck: error 1");
		if(pt->type == 1){
			if(pt->dr != NULL)printerr("treecheck: relate error 1");
			if(pt->dl == NULL)printerr("treecheck: relate error 2");
			if(pt->dl->a != pt)printerr("treecheck: relate error 3");
			if(pt->a == NULL)printerr("treecheck: relate error 4");
			if(pt->a->type == 1 && pt->a->dl != pt)
				printerr("treecheck: relate error 5");
			if(pt->a->type == 2 && !(pt->a->dl == pt || pt->a->dr == pt))
				printerr("treecheck: relate error 6");
			if(pt->time < pt->dl->time)printerr("treecheck: error 2");
			if(pt->val == pt->dl->val)printerr("treecheck: error 3");
			if(pt->dl->pos < blim || pt->dl->pos > ulim)
						printerr("treecheck: error 3a");
			if(pt->lno < 0 || pt->lno >= 2*sno-1)
					printerr("treecheck: lno error 0");
			if(pt->lno != pt->dl->lno)printerr("treecheck: lno error 1");
		}
		else if(pt->type == 2){
			if(pt->dr == NULL || pt->dl == NULL)
				printerr("treecheck: relate error 7");
			if(pt->lno == pt->dl->lno || pt->lno == pt->dr->lno)
					printerr("treecheck: lno error 2");
			if(pt->lno < 0 || pt->lno >= 2*sno-1)printerr("treecheck: lno error 3");
			if(!(pt->dr->a == pt && pt->dl->a == pt))
				printerr("treecheck: relate error 8");
			if(pt->a != NULL){
				if(pt->a->type == 1 && pt->a->dl != pt)
					printerr("treecheck: relate error 9");
				if(pt->a->type == 2 && !(pt->a->dl == pt || pt->a->dr == pt))
					printerr("treecheck: relate error 10");
			}
			else{
				if(pt->tlf != NULL)printerr("treecheck: mrca->tlf != NULL");
			}
			if(pt->time < pt->dl->time)printerr("treecheck: error 4");
			if(pt->time < pt->dr->time)printerr("treecheck: error 5");
			if(pt->val != pt->dl->val)printerr("treecheck: error 6");
			if(pt->val != pt->dr->val)printerr("treecheck: error 7");
			if(pt->dl->pos < blim || pt->dl->pos > ulim)
						printerr("treecheck: error 7a");
			if(pt->dr->pos < blim || pt->dr->pos > ulim)
						printerr("treecheck: error 7b");
			--nl;
		}
		else printerr("treecheck: error 7ba");
		if(pt->mark < -2 || pt->mark == 0 )
						printerr("treecheck: error 7c");
		ll = linswap_cal(pt);
		if(ll != pt->linswap)printerr("treecheck: linswap error");
		ltot += ll;
		sumswn += swn = swapable(pt);
		if(swn != pt->evswap)printerr("treecheck: evswap error");
		
		pt = pt->tlf;
		++ic;
		if(ic > 5000)printerr("treecheck: loop too long");
	}
	if(nl != 1)printerr("treecheck: error 8");
	if(blim == 0){
		if(ic != ulim+1-sno)printerr("treecheck: error 9");
	}
	else if(ic != ulim+1-(blim+sno))printerr("treecheck: error 10");
	if(pars->linswap_indic != ltot)printerr("treecheck: pars->linswap_indic wrong");
	if(sumswn != pars->evswap_indic)printerr("treecheck: pars->evswap_indic wrong");
	
	lwmax = -1.0;
	nwmax = -1.0;
	lwtot = 0.0;
	nwtot = 0.0;
	
	lastval1 = 0;
	lastval2 = 0;
	n3sum = 0;
	pairsum = 0;
	for(j=0;j<nc;++j){  /* this loop checks mutation parameters */
		cp = &(clist[j]);
		if(!ismark && cp->mark[0] != -1 && cp->mark[1] != -1)
				printerr("treecheck: mlin check, error 00");