📄 msvar1.3.c
字号:
nc[iloc],pars[iloc],1);
}
pars[iloc]->lik = lik_cal(pars[iloc]->base->tlf,
NULL,pars[iloc]);
priorcalc(pars[iloc],hpars);
pdiff = ldiff = tdiff = 0.0;
pdiff += pars[iloc]->prior - pars[iloc]->prior2;
ldiff += pars[iloc]->lik-pars[iloc]->lik2;
tdiff += pars[iloc]->tpr - pars[iloc]->tpf;
if(mtype == 22 && fabs(ldiff) > 0.00001){
printf("probs with mtype = 22; ldiff is %f\n",ldiff);
}
}
else if(twidtheta == 2){
choosepar(&ptlist,hpars,0,&mtype); /* This determines
type of updates for mutation and demographic
parameters. */
for(j=0;j<nloc;++j)pars[j]->lik2 = lik_cal(pars[j]->base->tlf,
NULL,pars[j]);
for(j=0;j<nloc;++j)copypars(pars[j]);
copyhpars(hpars);
/* if we are changing all of them need to change prior means,
but only want to set it once - so do it for the first */
for(j=0;j<nloc;++j)twidpars(&ptlist,pars[j],hpars,j==0 ? 1 : 0);
if(!ptlist.keeptime)
for(j=0;j<nloc;++j)tchange(pars[j],clist[j],mlist[j],&mm[j],
cmlist[j],&cmm[j]);
if((CHeck || rcheck)){
for(j=0;j<nloc;++j)checktree(list[j][0],nn[j],sno[j],clist[j],
nc[j],pars[j],1);
}
for(j=0;j<nloc;++j)pars[j]->lik = lik_cal(pars[j]->base->tlf,
NULL,pars[j]);
/* the reason for this is that when we change all loci, we change
the prior mean, and so everything should be the same */
for(j=0;j<nloc;++j)priorcalc(pars[j],hpars);
pdiff = ldiff = tdiff = 0.0;
for(j=0;j<nloc;++j)pdiff += pars[j]->prior - pars[j]->prior2;
if(fabs(pdiff) > 0.00001)printerr("main: abs(pdiff) too "
"large");
hpriorcalc(hpars);
pdiff += hpars->prior - hpars->prior2;
for(j=0;j<nloc;++j)ldiff += pars[j]->lik-pars[j]->lik2;
for(j=0;j<nloc;++j)tdiff += pars[j]->tpr - pars[j]->tpf;
if(mtype == 22 && fabs(ldiff) > 0.00001){
printf("probs with mtype = 22; ldiff is %f\n",ldiff);
}
}
else if(twidtheta == 3){
++htot;
for(j=0;j<nloc;++j)pars[j]->lik2 = lik_cal(pars[j]->base->tlf,
NULL,pars[j]);
for(j=0;j<nloc;++j)copypars(pars[j]);
copyhpars(hpars);
/*this is because it's not done in copyhpars*/
for(j=0;j<nloc;++j)pars[j]->prior2 = pars[j]->prior;
tdiff = twidhpars(&ptlist,hpars);
for(j=0;j<nloc;++j)twidpars_var(&ptlist,pars[j],hpars);
if(!ptlist.keeptime)
for(j=0;j<nloc;++j)tchange(pars[j],clist[j],mlist[j],&mm[j],
cmlist[j],&cmm[j]);
if((CHeck || rcheck)){
for(j=0;j<nloc;++j)checktree(list[j][0],nn[j],sno[j],clist[j],
nc[j],pars[j],1);
}
for(j=0;j<nloc;++j)pars[j]->lik = lik_cal(pars[j]->base->tlf,
NULL,pars[j]);
for(j=0;j<nloc;++j)priorcalc(pars[j],hpars);
pdiff = ldiff = 0.0; /*tdiff initialised in twidhpars */
for(j=0;j<nloc;++j)pdiff += pars[j]->prior - pars[j]->prior2;
hpriorcalc(hpars);
pdiff += hpars->prior - hpars->prior2;
for(j=0;j<nloc;++j)ldiff += pars[j]->lik-pars[j]->lik2;
for(j=0;j<nloc;++j)tdiff += pars[j]->tpr - pars[j]->tpf;
}
else{
if(!(mm[iloc] == 0 && cmm[iloc] == 0))printerr("main: mm/cmm != 0");
twiddle(list[iloc],&oldnn[iloc],&nn[iloc],clist[iloc],nc[iloc],mlist[iloc],&mm[iloc],
cmlist[iloc],&cmm[iloc],&mtype,pars[iloc]);
if(!Illegal && (CHeck || rcheck))checktree(list[iloc][0],nn[iloc],sno[iloc],clist[iloc],nc[iloc],
pars[iloc],1);
pdiff = 0.0;
ldiff = pars[iloc]->lik-pars[iloc]->lik2;
tdiff = pars[iloc]->tpr - pars[iloc]->tpf;
}
if(CHeck && !Illegal){
liknew = 0.0;
for(j=0;j<nloc;++j)liknew += lik_cal(list[j][0]->tlf,NULL,pars[j]);
ldiff2 = liknew - likold;
if(fabs(ldiff-ldiff2) > 1.0e-9){
printf("mtype is %d\n",mtype);
printerr("main: ldiff wrong");
}
}
l1 = ldiff + tdiff + pdiff;
for(j=0;j<nloc;++j)l1 = boundcheck(l1,nn[j],pars[j],hpars);
if(mtype >= 0 && mtype <= 30){
++mtlist[mtype];
++mtypecount;
}
if(!Illegal && ( l1 >= 0.0 || (l1 > -15.0 && gfsr8() < exp(l1)) )){
if(mtype >= 0 && mtype <= 30){
++mplist[mtype];
}
if(twidtheta == 0){
accept(mlist[iloc],&mm[iloc],cmlist[iloc],&cmm[iloc]);
if(CHeck || rcheck)checktree(list[iloc][0],nn[iloc],sno[iloc],
clist[iloc],nc[iloc],pars[iloc],0);
}
else if(twidtheta == 1){
accept(mlist[iloc],&mm[iloc],cmlist[iloc],&cmm[iloc]);
}
else for(j=0;j<nloc;++j)accept(mlist[j],&mm[j],cmlist[j],&cmm[j]);
if(CHeck)likold = liknew;
if(twidtheta == 3)hsucc += 1.0;
}
else{
if(twidtheta == 0){
untwiddle(mlist[iloc],&mm[iloc],cmlist[iloc],&cmm[iloc],pars[iloc],
&nn[iloc],oldnn[iloc]);
if(mtype == 0 || CHeck || rcheck)
checktree(list[iloc][0],nn[iloc],sno[iloc],clist[iloc],nc[iloc],pars[iloc],0);
}
else if(twidtheta == 1){
untwidpars(mlist[iloc],&mm[iloc],cmlist[iloc],
&cmm[iloc],pars[iloc]);
}
else if(twidtheta == 2){
untwidhpars(hpars);
for(j=0;j<nloc;++j)untwidpars(mlist[j],&mm[j],cmlist[j],&cmm[j],pars[j]);
}
else {
untwidhpars(hpars);
for(j=0;j<nloc;++j)pars[j]->prior = pars[j]->prior2;
for(j=0;j<nloc;++j)untwidpars(mlist[j],&mm[j],cmlist[j],&cmm[j],pars[j]);
}
mtype = 0;
}
if((iter+1)%Thin == 0){
++counter;iter=0;
if(counter%10 == 0){
mpout = fopen("mtype.dat","w");
for(j=0;j<30;++j){
fprintf(mpout,"%d %f %f\n", j,
(double)mtlist[j]/mtypecount,
mtlist[j] > 0 ? (double)mplist[j]/mtlist[j] : -1.0);
}
for(j=0;j<30;++j){mtlist[j] = 0;mplist[j] = 0;}
mtypecount = 0;
fclose(mpout);
}
lsum = 0.0;
for(j=0;j<nloc;++j){
lcheck1 = lik_cal(pars[j]->base->tlf,NULL,pars[j]);
lsum += lcheck1;
treesummary(nn[j],pars[j],list[j],&above,&mutbc,&blpc);
fprintf(lfile[j],"%d %d %e %f %f %f %f %f %f %f %f\n",counter,
nn[j]-(sno[j] + nc[j]),clist[j][nc[j]-1].p->time,log10(pars[j]->n0),
log10(pars[j]->n1),log10(pars[j]->mu),log10(pars[j]->tfa),lcheck1,
mutbc,blpc,above);
}
fprintf(hfile,"%d %f %f %f %f %f %f %f %f %f %f\n",
counter,hsucc/htot,lsum,hpars->ln0mean,hpars->ln0var,hpars->ln1mean,
hpars->ln1var,hpars->lmumean,hpars->lmuvar,hpars->ltfamean,
hpars->ltfavar);
hsucc = htot = 0.0;
if(counter%10 == 0){
for(j=0;j<nloc;++j)fflush(lfile[j]);
fflush(hfile);
closegfsr();
}
}
}
closegfsr();
}
void treesummary(int nn, Genpars *pars,Node **list,double *above, double
*mutbc,double *blpc)
{
Node *pp,*pt;
int mutno=0,coalno=0,cmutno=0,cabove=0,j,nm1;
double tsum1,tsuma;
pp = pars->base->tlf;
while(pp != NULL){
if(pp->type == 0)printerr("treesummary: pp->type == 0");
if(pp->type == 1){
++mutno;
if(coalno == 0)++cmutno;
}
if(pp->type == 2){
++coalno;
if(pp->time >= pars->tf)++cabove;
}
pp = pp->tlf;
}
/* *mutbc = cmutno/(double)mutno;*/
*above = cabove/(double)coalno;
tsum1 = 0.0;
nm1 = 0;
for(j=0;j<pars->sno;++j){
pt = list[j]->a;
while(1){
if(pt->type == 2)break;
++nm1;
pt = pt->a;
}
tsum1 += pt->time;
}
tsuma = 0.0;
for(j=pars->sno;j<nn;++j){
if(list[j]->type != 2)continue;
tsuma += 2.0*list[j]->time - list[j]->cp->l->time -
list[j]->cp->r->time;
}
if(mutno > 0)*mutbc = nm1/(double)mutno;
else *mutbc = -1.0;
*blpc = tsum1/tsuma;
return;
}
double boundcheck(double l1,int nn,Genpars *pars,Hierpars *hpars)
{
if(nn > 1995 )l1 = -20.0;
if(hpars->ln0var <= 0.0 || hpars->ln1var <= 0.0 || hpars->lmuvar <= 0.0 ||
hpars->ltfavar <= 0.0)l1 = -20.0;
return l1;
}
void read_params()
{
double up[10],sum;
FILE *initf;
int ic,j,kn0,kn1,kmu,ktfa,turno;
float ftest;
up[0] = 0.1;up[1] = 0.1;up[2] = 0.4;up[3] = 0.1;up[4] = 0.1;up[5] = 0.1,up[6] = 0.1;
initf = fopen("init_v_file","r");
if(initf == NULL)printerr("no init_v_file");
fscanf(initf,"%d",&turno);
if(turno < 0)turno = 0;
if(turno > 100000000)printerr("turno too big");
for(j=0;j<turno;++j)intrand();
ic = 0;
ic += fscanf(initf,"%d",&Ploidy);
if(Ploidy <= 0)printerr("ploidy is wrong");
if(Ploidy > 2)printf("Warning: ploidy is greater than 2 - "
"is this reasonable?\n");
ic += fscanf(initf,"%lf",&Init_gen);
for(j=0;j<Nloc;++j)ic += fscanf(initf,"%lf",&Init_n0[j]);
for(j=0;j<Nloc;++j)ic += fscanf(initf,"%lf",&Init_n1[j]);
for(j=0;j<Nloc;++j)ic += fscanf(initf,"%lf",&Init_mu[j]);
for(j=0;j<Nloc;++j)ic += fscanf(initf,"%lf",&Init_tfa[j]);
ic += fscanf(initf,"%d %d %d %d",&kn0,&kn1,&kmu,&ktfa);
if(kn0){
up[2] = 0.0;
}
if(kn1){
up[0] = 0.0;
up[2] = 0.0;
up[3] = 0.0;
up[4] = 0.0;
}
if(kmu){
up[0] = 0.0;
up[1] = 0.0;
up[2] = 0.0;
up[4] = 0.0;
up[5] = 0.0;
}
if(ktfa){
up[1] = 0.0;
up[2] = 0.0;
up[4] = 0.0;
up[6] = 0.0;
}
sum = 0.0;
for(j=0;j<7;++j)sum += up[j];
if(sum <= 0.000001){
for(j=0;j<7;++j)Usum[j] = 0.0;
}
else{
for(j=0;j<7;++j)up[j] /= sum;
Usum[0] = up[0];
for(j=1;j<7;++j)Usum[j] = Usum[j-1] + up[j];
}
ic += fscanf(initf,"%lf %lf",&Init_ln0mean,&Init_ln0var);
ic += fscanf(initf,"%lf %lf",&Init_ln1mean,&Init_ln1var);
ic += fscanf(initf,"%lf %lf",&Init_lmumean,&Init_lmuvar);
ic += fscanf(initf,"%lf %lf",&Init_ltfamean,&Init_ltfavar);
ic += fscanf(initf,"%lf %lf %lf %lf",&Init_ln0mm,&Init_ln0mv,&Init_ln0vm,&Init_ln0vv);
ic += fscanf(initf,"%lf %lf %lf %lf",&Init_ln1mm,&Init_ln1mv,&Init_ln1vm,&Init_ln1vv);
ic += fscanf(initf,"%lf %lf %lf %lf",&Init_lmumm,&Init_lmumv,&Init_lmuvm,&Init_lmuvv);
ic += fscanf(initf,"%lf %lf %lf %lf",&Init_ltfamm,&Init_ltfamv,&Init_ltfavm,&Init_ltfavv);
ic += fscanf(initf,"%d",&Mod_type);
if(!(Mod_type == 0 || Mod_type == 1))printerr("init_v_file: mod_type wrong");
ic += fscanf(initf,"%d",&Itno);
ic += fscanf(initf,"%d",&Thin);
if(ic != 4*Nloc+33){
printf("ic is %d\n",ic);
printerr("init_v_file probably too short");
}
if(Thin <= 0 )printerr("init_v_file: Thin seems funny");
fclose(initf);
}
void initpars(Genpars *pars,int j)
{
FILE *initf;
pars->pplus=pars->pplus2 = 0.5;
pars->n0_2 = pars->n0 = Init_n0[j];
pars->n1_2 = pars->n1 = Init_n1[j];
pars->mu2 = pars->mu = Init_mu[j];
pars->gen2 = pars->gen = Init_gen;
pars->tfa2 = pars->tfa = Init_tfa[j];
pars->which_model = Mod_type;
fillpars(pars);
}
void inithpars(Hierpars *hpars)
{
hpars->ln0mean = hpars->ln0mean2 = Init_ln0mean;
hpars->ln0mm = Init_ln0mm;
hpars->ln0mv = Init_ln0mv;
hpars->ln0var = hpars->ln0var2 = Init_ln0var;
hpars->ln0vm = Init_ln0vm;
hpars->ln0vv = Init_ln0vv;
hpars->ln1mean = hpars->ln1mean2 = Init_ln1mean;
hpars->ln1mm = Init_ln1mm;
hpars->ln1mv = Init_ln1mv;
hpars->ln1var = hpars->ln1var2 = Init_ln1var;
hpars->ln1vm = Init_ln1vm;
hpars->ln1vv = Init_ln1vv;
hpars->lmumean = hpars->lmumean2 = Init_lmumean;
hpars->lmumm = Init_lmumm;
hpars->lmumv = Init_lmumv;
hpars->lmuvar = hpars->lmuvar = Init_lmuvar;
hpars->lmuvm = Init_lmuvm;
hpars->lmuvv = Init_lmuvv;
hpars->ltfamean = hpars->ltfamean2 = Init_ltfamean;
hpars->ltfamm = Init_ltfamm;
hpars->ltfamv = Init_ltfamv;
hpars->ltfavar = hpars->ltfavar2 = Init_ltfavar;
hpars->ltfavm = Init_ltfavm;
hpars->ltfavv = Init_ltfavv;
}
void untwidhpars(Hierpars *hpars)
{
hpars->ln0mean = hpars->ln0mean2;
hpars->ln0var = hpars->ln0var2;
hpars->ln1mean = hpars->ln1mean2;
hpars->ln1var = hpars->ln1var2;
hpars->lmumean = hpars->lmumean2;
hpars->lmuvar = hpars->lmuvar2;
hpars->ltfamean = hpars->ltfamean2;
hpars->ltfavar = hpars->ltfavar2;
hpars->prior = hpars->prior2;
}
void copyhpars(Hierpars *hpars)
{
hpars->ln0mean2 = hpars->ln0mean;
hpars->ln0var2 = hpars->ln0var;
hpars->ln1mean2 = hpars->ln1mean;
hpars->ln1var2 = hpars->ln1var;
hpars->lmumean2 = hpars->lmumean;
hpars->lmuvar2 = hpars->lmuvar;
hpars->ltfamean2 = hpars->ltfamean;
hpars->ltfavar2 = hpars->ltfavar;
hpars->prior2 = hpars->prior;
}
void untwidpars(Node **mlist,int *mm,Cnode **cmlist,int *cmm,Genpars *pars)
{
int j;
pars->prior = pars->prior2;
pars->n0 = pars->n0_2;
pars->n1 = pars->n1_2;
pars->mu = pars->mu2;
pars->gen = pars->gen2;
pars->tfa = pars->tfa2;
pars->pplus = pars->pplus2;
pars->lwt_max = pars->lwt_max2;
pars->lwt_tot = pars->lwt_tot2;
pars->nwt_max = pars->nwt_max2;
pars->nwt_tot = pars->nwt_tot2;
fillpars(pars);
for(j=0;j<*mm;++j)restore(mlist[j]);
*mm = 0;
for(j=0;j<*cmm;++j)crestore(cmlist[j]);
*cmm = 0;
}
void choosepar(struct Gplist *ptlist,Hierpars *hpars,int oneonly,
int *mtype)
{
double dd,parwt,ddscale,rr,rr2,temp1;
*mtype = 0;
ptlist->keeptime
= ptlist->keepn0
= ptlist->keepn1
= ptlist->keepmu
= ptlist->keeptfa
= ptlist->keeppplus
= 1;
rr = gfsr8();
if(rr < 1.0){
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -