📄 msvar1.3.c
字号:
pv = clist[j].p;
if(pv == NULL || pv->type != 2)printerr("treecheck:mlin check, error 0");
tim1 = lwt_cal(pv,0);
if(cp->lwt[0] != tim1)printerr("treecheck: lwt[0] wrong");
if(lwmax < tim1)lwmax = tim1;
lwtot += tim1;
tim1 = lwt_cal(pv,1);
if(cp->lwt[1] != tim1)printerr("treecheck: lwt[1] wrong");
if(lwmax < tim1)lwmax = tim1;
lwtot += tim1;
tim1 = nwt_cal(pv);
if(cp->nwt != tim1)printerr("treecheck: nwt wrong");
if(nwmax < tim1)nwmax = tim1;
nwtot += tim1;
pp = pv->a;
if(pp == NULL){
if(!(clist[j].mlin[2][0] == 0 && clist[j].mlin[2][1] == 0 &&
clist[j].mlin[2][2] == 0)){
printerr("treecheck: mlin check, error 1");
}
}
else{
nminus = nplus = 0;
while(pp->type == 1){
if(pp->val - pp->dl->val == -1){
++nminus;
}
else ++nplus;
pp = pp->a;
}
if(!(cp->mlin[2][0] == nminus && cp->mlin[2][1] == nplus &&
cp->mlin[2][2] == nplus+nminus)){
printerr("treecheck: mlin check, error 2"); /* cf below */
}
}
pp = pv->dl;
nminus = nplus = 0;
while(pp->type == 1){
if(pp->val - pp->dl->val == -1){
++nminus;
}
else ++nplus;
pp = pp->dl;
}
if(!(cp->mlin[0][0] == nminus && cp->mlin[0][1] == nplus &&
cp->mlin[0][2] == nplus+nminus)){
printerr("treecheck: mlin check, error 3");
/* need to change last bit if more than
1-step mutation */
}
pp = pv->dr;
nminus = nplus = 0;
while(pp->type == 1){
if(pp->val - pp->dl->val == -1){
++nminus;
}
else ++nplus;
pp = pp->dl;
}
if(!(cp->mlin[1][0] == nminus && cp->mlin[1][1] == nplus && /* cf above */
cp->mlin[1][2] == nplus+nminus)){
printerr("treecheck: mlin check, error 4");
}
lastval1 = cp->mlin[0][0]*cp->mlin[0][1] > 0;
if(lastval1 != clist[j].linmut[0])printerr("treecheck: mlin check, error 5");
pairsum += lastval1;
lastval1 = cp->mlin[1][0]*cp->mlin[1][1] > 0;
if(lastval1 != clist[j].linmut[1])printerr("treecheck: mlin check, error 6");
pairsum += lastval1;
lastval2 = cal_nodeprob(cp->mlin,clist[j].u == NULL);
if(lastval2 != clist[j].nodemut)printerr("treecheck: mlin check, error 7");
n3sum += lastval2;
}
if(pars->nodemut_indic != n3sum)printerr("treecheck:mlin check, error, 8");
if(pars->linmut_indic != pairsum)printerr("treecheck: mlin check, error, 9");
if(lwmax > pars->lwt_max)printerr("treecheck: pars->lwt_max wrong");
if(fabs(pars->lwt_tot - lwtot)/lwtot > 1.0e-5)printerr("treecheck: pars->lwt_tot wrong");
if(fabs(pars->nwt_tot - nwtot)/nwtot > 1.0e-5)printerr("treecheck: pars->nwtot wrong");
if(nwmax > pars->nwt_max)printerr("treecheck: pars->nwt_max wrong");
/* for rounding errors - also correction of nwt_max, lwt_max*/
pars->lwt_max = lwmax;
pars->lwt_tot = lwtot;
pars->nwt_tot = nwtot;
pars->nwt_max = nwmax;
ncheck = 0;
cpt = &clist[nc-1];
while(1){
if(cpt->cpos < 0 || cpt->cpos >= nc)printerr("treecheck: clist check, error 3");
if(cpt->p == NULL)printerr("treecheck: clist check, error 4");
if(cpt->cpos == nc-1){
if(cpt->u != NULL)printerr("treecheck: clist check, error 5");
if(cpt->tf != NULL)printerr("treecheck: clist check, error 5b");
}
if(cpt->cpos != nc-1){
if(cpt->u == NULL)printerr("treecheck: clist check, error 6");
if(cpt->u->cp == NULL)printerr("treecheck: clist check, error 7");
if(cpt->u->cp->l == cpt->u->cp->r)
printerr("treecheck: clist check, error 8");
if(!(cpt->u->cp->l == cpt->p || cpt->u->cp->r == cpt->p))
printerr("treecheck: clist check, error 9");
if(cpt->tf == NULL)printerr("trecheck: clist check, error 9b");
if(cpt->tf->cp == NULL)printerr("treecheck: clist check, error 9c");
if(cpt->tf->cp->tr != cpt->p)printerr("treecheck: clist check, error 9d");
if(cpt->tf->time <= cpt->p->time)
printerr("treecheck: clist check, error 9e");
if(cpt->tr->time >= cpt->p->time)
printerr("treecheck: clist check, error 9f");
}
if(cpt->l == NULL)printerr("treecheck: clist check, error 9e");
if(cpt->l->type != 0){
if(cpt->l->cp == NULL)printerr("treecheck: clist check, error 10");
if(cpt->l->cp->u != cpt->p)printerr("treecheck: clist check, error11");
}
if(cpt->r == NULL)printerr("treecheck: clist check, error 11b");
if(cpt->r->type != 0){
if(cpt->r->cp == NULL)printerr("treecheck: clist check, error 12");
if(cpt->r->cp->u != cpt->p)printerr("treecheck: clist check, error13");
}
++ncheck;
if(cpt->tr == NULL)printerr("treecheck: clist check, error 13b");
if(cpt->tr->type != 0){
if(cpt->tr->cp == NULL)printerr("treecheck: clist check, error 13c");
if(cpt->tr->cp->tf != cpt->p)printerr("treecheck: clist check, error 14");
cpt = cpt->tr->cp;
}
else break;
}
if(ncheck != nc)printerr("treecheck: clist check, error 14b");
pscheck = MUTLIN_ADD_P + MUTNODE_ADD_P;
if(pars->evswap_indic)pscheck += ESWAP_PROB;
if(pars->linmut_indic)pscheck += MUTLIN_DEL_P;
if(pars->nodemut_indic)pscheck += MUTNODE_DEL_P;
if(pars->linswap_indic)pscheck += LSWAP_PROB;
if(fabs(pscheck - pars->pspace) > 1.0e-9)printerr("treecheck: clist check, error 17");
/* Extra checks */
pp = startpp;
while(pp->type != 2)pp = pp->tlf;
if(pp->cp->tr->type != 0)printerr("treecheck extras: not first coalescence");
cpp = pp->cp;
while(cpp != NULL){
pp = pp->tlf;
if(pp==NULL)break;
while(pp->type != 2)pp=pp->tlf;
if(pp->cp->p != cpp->tf)printerr("treecheck extras: pp->cp != cpp->tf");
if(pp->cp->tr != cpp->p)printerr("trecheck extras: pp->cp->tr != cpp");
cpp = cpp->tf->cp;
}
}
double tdens(double t1,double t2,int nl,double theta,double r,double tf,int which_model)
{
if(which_model == 0)return tdensL(t1,t2,nl,theta,r,tf);
else if(which_model == 1)return tdensE(t1,t2,nl,theta,r,tf);
else printerr("tdens: which_model wrong");
}
double tdensL(double t1,double t2,int nl,double theta,double r,double tf)
{
/* this gives the log density. It is slightly inefficient for last two cases because
usually we want the ratio and coeff cancels out */
double fx,coeff;
if(t1 >= t2)printerr("tdensL: t1 >= t2");
if(fabs(r-1.0) < 0.00001) tf = 0.0; /* to avoid division by zero */
if(tf > 1.e10){r=1;tf=0.0;} /* to trap rounding errors with large tf */
if(t1 < tf && t2 < tf){
fx = nl*theta*(t2-t1)*0.5 + (nl*(nl-1.0)*r*tf/(2.0*(1.0-r)))*
log((r*tf+t2-r*t2)/(r*tf+t1-r*t1));
coeff = (nl*theta*0.5 + (nl-1.0)*nl*r*tf*0.5/(r*tf+t2-r*t2));
}
else if(t1 < tf && t2 >= tf){
fx = nl*theta*(tf-t1)*0.5 + (nl*(nl-1.0)*r*tf/(2.0*(1.0-r)))*
log(tf/(r*tf+t1-r*t1)) + nl*(r*(nl-1)+theta)*0.5*(t2-tf);
coeff = nl*(r*(nl-1)+theta)*0.5;
}
else if(t1 >= tf && t2 >= tf){
coeff = nl*(r*(nl-1)+theta)*0.5;
fx = coeff*(t2-t1);
}
return log(coeff) - fx;
}
double tdensE(double t1,double t2,int nl,double theta,double r,double tf)
{
/* this gives the log density. It is slightly inefficient for last two cases
because
usually we want the ratio and coeff cancels out */
double fx,coeff,logr,rtf,powrt2tf;
logr = log(r);
rtf = r*tf;
if(t1 >= t2)printerr("tdens: t1 >= t2");
if(fabs(r-1.0) < 0.00001) tf = 0.0; /* to avoid division by zero */
if(tf > 1.e10){r=1;tf=0.0;} /* to trap rounding errors with large tf */
if(t1 < tf && t2 < tf){
powrt2tf = pow(r,t2/tf);
fx = 0.5*nl*(theta*(t1-t2) + (nl-1)*tf/logr*(pow(r,t1/tf)-
powrt2tf));
/* this is -probfunc -> hence different from fx in linear model */
coeff = powrt2tf*nl*(nl-1)/2+nl*theta/2;
}
else if(t1 < tf && t2 >= tf){
fx = (nl*(rtf*(logr-1)*(nl-1) + tf*pow(r,t1/tf)*(nl-1) +
logr*(theta*t1
+t2*(r*(1 -nl) - theta))))/(2*logr);
/* this is -probfunc -> hence different from fx in linear model */
coeff = nl*(r*(nl-1)+theta)*0.5;
}
else if(t1 >= tf && t2 >= tf){
coeff = nl*(r*(nl-1)+theta)*0.5;
fx = -coeff*(t2-t1); /* note different sign from linear model */
}
return log(coeff) + fx; /* note different sign from linear model */
}
double lik_cal(Node *pt,Node *out,Genpars *pars)
{
if(Illegal){
return -1.0e100;
}
if(pars->which_model == 0)return lik_calL(pt,out,pars);
else if(pars->which_model == 1)return lik_calE(pt,out,pars);
else printerr("lik_cal: which_model wrong");
}
double lik_calL(Node *pt,Node *out,Genpars *pars)
{
int nl,nlcheck,mdir;
double lik,t1,t2,pcdfx1,pmdfx1,fx,complik,l1,pmval;
double logmut,logplus,logminus,logr;
double theta,r,tf;
Node *stp;
if(pt == NULL)printerr("lik_calL: pt == NULL");
if(pt->type == 0)pt = pars->base->tlf; /* lik_calL assumes
we start from the first event, not baseline */
if(out == NULL)stp = NULL;
else stp = out->tlf;
lik = 0.0;
complik = 0.0;
nlcheck = pt->nlin;
theta = pars->theta;
r = pars->r;
tf = pars->tf;
if(fabs(r-1.0) < 0.0001)tf = 0.0;
if(tf > 1.e10){r=1;tf=0.0;} /* to trap rounding errors with large tf */
logmut = log(theta*0.5);
if(pars->pplus > 0.0)logplus = log(pars->pplus);
else logplus = -10000.0;
if(1.0 - pars->pplus > 0.0)logminus = log(1.0-pars->pplus);
logr = log(r);
while(1){
nl = pt->nlin;
if(nl != nlcheck)printerr("lik_calL: nlcheck wrong");
if(pt->time < pt->tlr->time)printerr("lik_calL: time wrong");
if(pt->time == pt->tlr->time){ /* bodge to fix times */
printerr2("lik_calL: time wrong - equal to, try to fix",0);
if(pt->tlf == NULL){
pt->time += 0.00001;
}
else if(pt->tlf == stp)printerr
("lik_calL: time wrong - equal to, can't fix");
/* because this might make likelihoods inconsistent */
else{
pt->time = pt->tlr->time +
0.5*(pt->tlf->time - pt->tlr->time);
if(pt->time == pt->tlr->time)printerr
("lik_calL: time wrong - equal to, can't fix");
}
}
t1 = pt->tlr->time;
t2 = pt->time;
if(t1 < tf && t2 < tf){
fx = nl*theta*(t2-t1)*0.5 + (nl*(nl-1.0)*r*tf/(2.0*(1.0-r)))*
log((r*tf+t2-r*t2)/(r*tf+t1-r*t1));
if(pt->type == 1){
pmdfx1 = theta*0.5;
mdir = pt->dl->val - pt->val;
if(mdir == 1)pmval = logplus;
else if(mdir == -1)pmval = logminus;
else printerr("lik_calL: mdir bigger than 1 step");
l1 = logmut + pmval -fx;
/* pmval chance of mut in that direction, 1/nl choose that
lineage */
lik += l1;
}
else if(pt->type == 2){
pcdfx1 = r*tf/(r*tf + t2 - r*t2);
l1 = log(pcdfx1) - fx;
/*last bit - prob choosing that
pair of lins */
lik += l1;
--nlcheck;
}
}
else if(t1 < tf && t2 >= tf){
fx = nl*theta*(tf-t1)*0.5 + (nl*(nl-1.0)*r*tf/(2.0*(1.0-r)))*
log(tf/(r*tf+t1-r*t1)) + nl*(r*(nl-1)+theta)*0.5*(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_calL: mdir bigger than 1 step");
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;
/*last bit - prob choosing that
pair of lins */
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_calL: 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 lik_calE(Node *pt,Node *out,Genpars *pars)
{
int nl,nlcheck,mdir;
double lik,t1,t2,pcdfx1,pmdfx1,fx,complik,l1,pmval,rtf,logrvaltf;
double logmut,logplus,logminus,logr;
double theta,r,tf;
Node *stp;
if(pt == NULL)printerr("lik_cal: pt == NULL");
if(pt->type == 0)pt = pars->base->tlf; /* lik_cal assumes
we start from the first event, not baseline */
if(out == NULL)stp = NULL;
else stp = out->tlf;
lik = 0.0;
complik = 0.0;
nlcheck = pt->nlin;
theta = pars->theta;
r = pars->r;
tf = pars->tf;
if(fabs(r-1.0) < 0.0001)tf = 0.0;
if(tf > 1.e10){r=1;tf=0.0;} /* to trap rounding errors with large tf */
logmut = log(theta*0.5);
if(pars->pplus > 0.0)logplus = log(pars->pplus);
else logplus = -10000.0;
if(1.0 - pars->pplus > 0.0)logminus = log(1.0-pars->pplus);
logr = log(r);
rtf = r*tf;
while(1){
nl = pt->nlin;
if(nl != nlcheck)printerr("old_lik_cal: nlcheck wrong");
if(pt->time < pt->tlr->time)printerr("old_lik_cal: time wrong");
if(pt->time == pt->tlr->time){ /* bodge to fix times */
printerr2("lik_calL: time wrong - equal to, try to fix",0);
if(pt->tlf == NULL){
pt->time += 0.00001;
}
else if(pt->tlf == stp)printerr
("lik_calL: time wrong - equal to, can't fix");
/* because this might make likelihoods inconsistent */
else{
pt->time = pt->tlr->time +
0.5*(pt->tlf->time - pt->tlr->time);
if(pt->time == pt->tlr->time)printerr
("lik_calL: time wrong - equal to, can't fix");
}
}
t1 = pt->tlr->time;
t2 = pt->time;
if(t1 < tf && t2 < tf){
/* fx = (nl*theta*t1)/2 - (nl*theta*t2)/2 +
((-1 + nl)*nl*pow(r,t1/tf)*tf)/(2*logr) -
((-1 + nl)*nl*pow(r,t2/tf)*tf)/(2*logr); */
fx = 0.5*nl*(theta*(t1-t2) + (nl-1)*tf/logr*(pow(r,t1/tf)-
pow(r,t2/tf)));
/* this is -probfunc -> hence different from fx in linear model */
if(pt->type == 1){
/* pmdfx1 = theta*0.5; */
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 = logmut + pmval +fx;
/* pmval chance of mut in that direction, 1/nl choose that
lineage */
lik += l1;
}
else if(pt->type == 2){
/* pcdfx1 = pow(r,t2/tf); can avoid pow by multiplying log with
t2/tf */
l1 = t2/tf*logr + fx;
lik += l1;
--nlcheck;
}
}
else if(t1 < tf && t2 >= tf){
/* fx = (nl*(rtf - nl*rtf - rpow*tf + nl*rpow*tf - rtf*logr +
nl*rtf*logr +
theta*t1*logr + r*t2*logr - nl*r*t2*logr -
theta*t2*logr))/(2*logr); */
fx = (nl*(rtf*(logr-1)*(nl-1) + tf*pow(r,t1/tf)*(nl-1) +
logr*(theta*t1
+t2*(r*(1 -nl) - theta))))/(2*logr);
/* this is -probfunc -> hence different from fx in linear model */
/* rtf = r*tf
rpow = pow(r,t1/tf)
logr = log(r) */
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");
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -