|
version 1.264, 2017/04/26 06:01:29
|
version 1.266, 2017/05/13 07:26:12
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.266 2017/05/13 07:26:12 brouard |
| |
Summary: Version 0.99r13 (improvements and bugs fixed) |
| |
|
| |
Revision 1.265 2017/04/26 16:22:11 brouard |
| |
Summary: imach 0.99r13 Some bugs fixed |
| |
|
| Revision 1.264 2017/04/26 06:01:29 brouard |
Revision 1.264 2017/04/26 06:01:29 brouard |
| Summary: Labels in graphs |
Summary: Labels in graphs |
| |
|
|
Line 2654 Earliest age to start was %d-%d=%d, ncvl
|
Line 2660 Earliest age to start was %d-%d=%d, ncvl
|
| max=vector(1,nlstate); |
max=vector(1,nlstate); |
| meandiff=vector(1,nlstate); |
meandiff=vector(1,nlstate); |
| |
|
| dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms; |
dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms; |
| oldm=oldms; savm=savms; |
oldm=oldms; savm=savms; |
| |
|
| /* Starting with matrix unity */ |
/* Starting with matrix unity */ |
| for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
| } |
} |
| |
|
|
Line 2733 Earliest age to start was %d-%d=%d, ncvl
|
Line 2739 Earliest age to start was %d-%d=%d, ncvl
|
| /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, ageminpar, agemaxpar, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ |
/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, ageminpar, agemaxpar, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ |
| /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ |
/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ |
| out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij)); /* Bug Valgrind */ |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij)); /* Bug Valgrind */ |
| |
/* if((int)age == 70){ */ |
| |
/* printf(" Backward prevalim age=%d agefin=%d \n", (int) age, (int) agefin); */ |
| |
/* for(i=1; i<=nlstate+ndeath; i++) { */ |
| |
/* printf("%d newm= ",i); */ |
| |
/* for(j=1;j<=nlstate+ndeath;j++) { */ |
| |
/* printf("%f ",newm[i][j]); */ |
| |
/* } */ |
| |
/* printf("oldm * "); */ |
| |
/* for(j=1;j<=nlstate+ndeath;j++) { */ |
| |
/* printf("%f ",oldm[i][j]); */ |
| |
/* } */ |
| |
/* printf(" pmmij "); */ |
| |
/* for(j=1;j<=nlstate+ndeath;j++) { */ |
| |
/* printf("%f ",pmmij[i][j]); */ |
| |
/* } */ |
| |
/* printf("\n"); */ |
| |
/* } */ |
| |
/* } */ |
| savm=oldm; |
savm=oldm; |
| oldm=newm; |
oldm=newm; |
| |
|
| for(j=1; j<=nlstate; j++){ |
for(j=1; j<=nlstate; j++){ |
| max[j]=0.; |
max[j]=0.; |
| min[j]=1.; |
min[j]=1.; |
|
Line 2785 Oldest age to start was %d-%d=%d, ncvloo
|
Line 2810 Oldest age to start was %d-%d=%d, ncvloo
|
| double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) |
double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) |
| { |
{ |
| /* According to parameters values stored in x and the covariate's values stored in cov, |
/* According to parameters values stored in x and the covariate's values stored in cov, |
| computes the probability to be observed in state j being in state i by appying the |
computes the probability to be observed in state j (after stepm years) being in state i by appying the |
| model to the ncovmodel covariates (including constant and age). |
model to the ncovmodel covariates (including constant and age). |
| lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc] |
lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc] |
| and, according on how parameters are entered, the position of the coefficient xij(nc) of the |
and, according on how parameters are entered, the position of the coefficient xij(nc) of the |
|
Line 2794 double **pmij(double **ps, double *cov,
|
Line 2819 double **pmij(double **ps, double *cov,
|
| j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel |
j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel |
| Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation, |
Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation, |
| sums on j different of i to get 1-pii/pii, deduces pii, and then all pij. |
sums on j different of i to get 1-pii/pii, deduces pii, and then all pij. |
| Outputs ps[i][j] the probability to be observed in j being in j according to |
Outputs ps[i][j] or probability to be observed in j being in i according to |
| the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij] |
the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij] |
| |
Sum on j ps[i][j] should equal to 1. |
| */ |
*/ |
| double s1, lnpijopii; |
double s1, lnpijopii; |
| /*double t34;*/ |
/*double t34;*/ |
|
Line 2859 double **pmij(double **ps, double *cov,
|
Line 2885 double **pmij(double **ps, double *cov,
|
| /* |
/* |
| for(i=1; i<= npar; i++) printf("%f ",x[i]); |
for(i=1; i<= npar; i++) printf("%f ",x[i]); |
| goto end;*/ |
goto end;*/ |
| return ps; |
return ps; /* Pointer is unchanged since its call */ |
| } |
} |
| |
|
| /*************** backward transition probabilities ***************/ |
/*************** backward transition probabilities ***************/ |
|
Line 2868 double **pmij(double **ps, double *cov,
|
Line 2894 double **pmij(double **ps, double *cov,
|
| /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ |
/* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ |
| double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) |
double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) |
| { |
{ |
| /* Computes the backward probability at age agefin and covariate ij |
/* Computes the backward probability at age agefin and covariate combination ij. In fact cov is already filled and x too. |
| * and returns in **ps as well as **bmij. |
* Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij. |
| */ |
*/ |
| int i, ii, j,k; |
int i, ii, j,k; |
| |
|
|
Line 2886 double **pmij(double **ps, double *cov,
|
Line 2912 double **pmij(double **ps, double *cov,
|
| |
|
| agefin=cov[2]; |
agefin=cov[2]; |
| /* bmij *//* age is cov[2], ij is included in cov, but we need for |
/* bmij *//* age is cov[2], ij is included in cov, but we need for |
| the observed prevalence (with this covariate ij) */ |
the observed prevalence (with this covariate ij) at beginning of transition */ |
| dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); |
/* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
| |
pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */ |
| |
/* outputs pmmij which is a stochastic matrix */ |
| /* We do have the matrix Px in savm and we need pij */ |
/* We do have the matrix Px in savm and we need pij */ |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| sumnew=0.; /* w1 p11 + w2 p21 only on live states */ |
sumnew=0.; /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */ |
| for (ii=1;ii<=nlstate;ii++){ |
for (ii=1;ii<=nlstate;ii++){ |
| sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; |
/* sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; */ |
| |
sumnew+=pmmij[ii][j]*prevacurrent[(int)agefin][ii][ij]; /* Yes prevalence at beginning of transition */ |
| } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ |
} /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ |
| for (ii=1;ii<=nlstate+ndeath;ii++){ |
if(sumnew >= 1.e-10){ |
| if(sumnew >= 1.e-10){ |
for (ii=1;ii<=nlstate+ndeath;ii++){ |
| /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ |
/* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ |
| /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
/* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
| /* }else if(agefin >= agemaxpar+stepm/YEARM){ */ |
/* }else if(agefin >= agemaxpar+stepm/YEARM){ */ |
| /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
/* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
| /* }else */ |
/* }else */ |
| doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
| }else{ |
} /*End ii */ |
| ; |
}else{ /* We put the identity matrix */ |
| /* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */ |
for (ii=1;ii<=nlstate+ndeath;ii++){ |
| } |
doldm[ii][j]=(ii==j ? 1. : 0.0); |
| } /*End ii */ |
} /*End ii */ |
| } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ |
/* printf("Problem internal bmij A: sum_i w_i*p_ij=N.j/N.. <1.e-10 i=%d, j=%d, sumnew=%lf,agefin=%d\n",ii,j,sumnew, (int)agefin); */ |
| /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */ |
} |
| bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ |
} /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] or identity*/ |
| |
/* left Product of this diag matrix by dsavm=Px (dnewm=dsavm*doldm) */ |
| |
/* bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /\* Bug Valgrind *\/ */ |
| |
bbmij=matprod2(dnewm, pmmij,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ |
| /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ |
/* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ |
| /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
/* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
| /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
/* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
| /* left Product of this matrix by diag matrix of prevalences (savm) */ |
/* left Product of this matrix by diag matrix of prevalences (savm) */ |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| |
sumnew=0.; |
| for (ii=1;ii<=nlstate+ndeath;ii++){ |
for (ii=1;ii<=nlstate+ndeath;ii++){ |
| |
sumnew+=prevacurrent[(int)agefin][ii][ij]; |
| dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); |
dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); |
| } |
} |
| } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */ |
/* if(sumnew <0.9){ */ |
| ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ |
/* printf("Problem internal bmij B: sum on i wi <0.9: j=%d, sum_i wi=%lf,agefin=%d\n",j,sumnew, (int)agefin); */ |
| |
/* } */ |
| |
} /* End j, At the end dsavm is diag[(w_i)] */ |
| |
/* What if dsavm doesn't sum ii to 1? */ |
| |
/* ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /\* Bug Valgrind *\/ */ |
| |
ps=matprod2(ps, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ |
| /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ |
/* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ |
| /* end bmij */ |
/* end bmij */ |
| return ps; |
return ps; /*pointer is unchanged */ |
| } |
} |
| /*************** transition probabilities ***************/ |
/*************** transition probabilities ***************/ |
| |
|
|
Line 3148 double ***hpxij(double ***po, int nhstep
|
Line 3187 double ***hpxij(double ***po, int nhstep
|
| /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */ |
/* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */ |
| double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) |
double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) |
| { |
{ |
| /* Computes the transition matrix starting at age 'age' over |
/* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over |
| 'nhstepm*hstepm*stepm' months (i.e. until |
'nhstepm*hstepm*stepm' months (i.e. until |
| age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
| nhstepm*hstepm matrices. |
nhstepm*hstepm matrices. |
|
Line 3156 double ***hbxij(double ***po, int nhstep
|
Line 3195 double ***hbxij(double ***po, int nhstep
|
| (typically every 2 years instead of every month which is too big |
(typically every 2 years instead of every month which is too big |
| for the memory). |
for the memory). |
| Model is determined by parameters x and covariates have to be |
Model is determined by parameters x and covariates have to be |
| included manually here. |
included manually here. Then we use a call to bmij(x and cov) |
| |
The addresss of po (p3mat allocated to the dimension of nhstepm) should be stored for output |
| */ |
*/ |
| |
|
| int i, j, d, h, k; |
int i, j, d, h, k; |
| double **out, cov[NCOVMAX+1]; |
double **out, cov[NCOVMAX+1], **bmij(); |
| double **newm; |
double **newm, ***newmm; |
| double agexact; |
double agexact; |
| double agebegin, ageend; |
double agebegin, ageend; |
| double **oldm, **savm; |
double **oldm, **savm; |
| |
|
| oldm=oldms;savm=savms; |
newmm=po; /* To be saved */ |
| |
oldm=oldms;savm=savms; /* Global pointers */ |
| /* Hstepm could be zero and should return the unit matrix */ |
/* Hstepm could be zero and should return the unit matrix */ |
| for (i=1;i<=nlstate+ndeath;i++) |
for (i=1;i<=nlstate+ndeath;i++) |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
|
Line 3180 double ***hbxij(double ***po, int nhstep
|
Line 3220 double ***hbxij(double ***po, int nhstep
|
| newm=savm; |
newm=savm; |
| /* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
| cov[1]=1.; |
cov[1]=1.; |
| agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ |
agexact=age-((h-1)*hstepm + (d))*stepm/YEARM; /* age just before transition, d or d-1? */ |
| /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ |
/* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ |
| cov[2]=agexact; |
cov[2]=agexact; |
| if(nagesqr==1) |
if(nagesqr==1) |
| cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
| for (k=1; k<=cptcovn;k++) |
for (k=1; k<=cptcovn;k++){ |
| cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ |
| /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ |
/* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ |
| |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
| |
/* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
| |
|
| |
} |
| for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ |
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ |
| /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
| cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
|
Line 3196 double ***hbxij(double ***po, int nhstep
|
Line 3240 double ***hbxij(double ***po, int nhstep
|
| cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
| /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ |
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ |
| |
|
| |
|
| /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
| /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
| /* Careful transposed matrix */ |
/* Careful transposed matrix */ |
| /* age is in cov[2] */ |
/* age is in cov[2], prevacurrent at beginning of transition. */ |
| /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ |
/* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ |
| /* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ |
/* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ |
| out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ |
out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ |
|
Line 4308 void freqsummary(char fileres[], double
|
Line 4351 void freqsummary(char fileres[], double
|
| int firstpass, int lastpass, int stepm, int weightopt, char model[]) |
int firstpass, int lastpass, int stepm, int weightopt, char model[]) |
| { /* Some frequencies as well as proposing some starting values */ |
{ /* Some frequencies as well as proposing some starting values */ |
| |
|
| int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0; |
int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0, s1=1, s2=1; |
| int iind=0, iage=0; |
int iind=0, iage=0; |
| int mi; /* Effective wave */ |
int mi; /* Effective wave */ |
| int first; |
int first; |
|
Line 4387 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4430 Title=%s <br>Datafile=%s Firstpass=%d La
|
| k2cpt=0; |
k2cpt=0; |
| |
|
| if(cptcoveff == 0 ) |
if(cptcoveff == 0 ) |
| nl=1; /* Constant model only */ |
nl=1; /* Constant and age model only */ |
| else |
else |
| nl=2; |
nl=2; |
| |
|
| |
/* if a constant only model, one pass to compute frequency tables and to write it on ficresp */ |
| |
/* Loop on nj=1 or 2 if dummy covariates j!=0 |
| |
* Loop on j1(1 to 2**cptcoveff) covariate combination |
| |
* freq[s1][s2][iage] =0. |
| |
* Loop on iind |
| |
* ++freq[s1][s2][iage] weighted |
| |
* end iind |
| |
* if covariate and j!0 |
| |
* headers Variable on one line |
| |
* endif cov j!=0 |
| |
* header of frequency table by age |
| |
* Loop on age |
| |
* pp[s1]+=freq[s1][s2][iage] weighted |
| |
* pos+=freq[s1][s2][iage] weighted |
| |
* Loop on s1 initial state |
| |
* fprintf(ficresp |
| |
* end s1 |
| |
* end age |
| |
* if j!=0 computes starting values |
| |
* end compute starting values |
| |
* end j1 |
| |
* end nl |
| |
*/ |
| for (nj = 1; nj <= nl; nj++){ /* nj= 1 constant model, nl number of loops. */ |
for (nj = 1; nj <= nl; nj++){ /* nj= 1 constant model, nl number of loops. */ |
| if(nj==1) |
if(nj==1) |
| j=0; /* First pass for the constant */ |
j=0; /* First pass for the constant */ |
| else |
else{ |
| j=cptcoveff; /* Other passes for the covariate values */ |
j=cptcoveff; /* Other passes for the covariate values */ |
| |
} |
| first=1; |
first=1; |
| for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */ |
for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all covariates combination of the model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */ |
| posproptt=0.; |
posproptt=0.; |
| /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
| scanf("%d", i);*/ |
scanf("%d", i);*/ |
| for (i=-5; i<=nlstate+ndeath; i++) |
for (i=-5; i<=nlstate+ndeath; i++) |
| for (jk=-5; jk<=nlstate+ndeath; jk++) |
for (s2=-5; s2<=nlstate+ndeath; s2++) |
| for(m=iagemin; m <= iagemax+3; m++) |
for(m=iagemin; m <= iagemax+3; m++) |
| freq[i][jk][m]=0; |
freq[i][s2][m]=0; |
| |
|
| for (i=1; i<=nlstate; i++) { |
for (i=1; i<=nlstate; i++) { |
| for(m=iagemin; m <= iagemax+3; m++) |
for(m=iagemin; m <= iagemax+3; m++) |
|
Line 4421 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4489 Title=%s <br>Datafile=%s Firstpass=%d La
|
| /* dateintsum=0; */ |
/* dateintsum=0; */ |
| /* k2cpt=0; */ |
/* k2cpt=0; */ |
| |
|
| /* For that combination of covariate j1, we count and print the frequencies in one pass */ |
/* For that combination of covariates j1 (V4=1 V3=0 for example), we count and print the frequencies in one pass */ |
| for (iind=1; iind<=imx; iind++) { /* For each individual iind */ |
for (iind=1; iind<=imx; iind++) { /* For each individual iind */ |
| bool=1; |
bool=1; |
| if(j !=0){ |
if(j !=0){ |
|
Line 4437 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4505 Title=%s <br>Datafile=%s Firstpass=%d La
|
| /* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
/* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
| /* }else */ |
/* }else */ |
| if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */ |
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */ |
| /* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */ |
/* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */ |
| bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ |
bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ |
| /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", |
/* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", |
| bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), |
bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), |
|
Line 4448 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4516 Title=%s <br>Datafile=%s Firstpass=%d La
|
| } /* cptcovn > 0 */ |
} /* cptcovn > 0 */ |
| } /* end any */ |
} /* end any */ |
| }/* end j==0 */ |
}/* end j==0 */ |
| if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */ |
if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ |
| /* for(m=firstpass; m<=lastpass; m++){ */ |
/* for(m=firstpass; m<=lastpass; m++){ */ |
| for(mi=1; mi<wav[iind];mi++){ /* For that wave */ |
for(mi=1; mi<wav[iind];mi++){ /* For that wave */ |
| m=mw[mi][iind]; |
m=mw[mi][iind]; |
|
Line 4510 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4578 Title=%s <br>Datafile=%s Firstpass=%d La
|
| |
|
| |
|
| /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
| pstamp(ficresp); |
if(cptcoveff==0 && nj==1) /* no covariate and first pass */ |
| |
pstamp(ficresp); |
| if (cptcoveff>0 && j!=0){ |
if (cptcoveff>0 && j!=0){ |
| |
pstamp(ficresp); |
| printf( "\n#********** Variable "); |
printf( "\n#********** Variable "); |
| fprintf(ficresp, "\n#********** Variable "); |
fprintf(ficresp, "\n#********** Variable "); |
| fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
|
Line 4539 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4609 Title=%s <br>Datafile=%s Firstpass=%d La
|
| fprintf(ficlog, "**********\n"); |
fprintf(ficlog, "**********\n"); |
| } |
} |
| fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
| |
if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ |
| |
fprintf(ficresp, " Age"); |
| |
if(nj==2) for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| for(i=1; i<=nlstate;i++) { |
for(i=1; i<=nlstate;i++) { |
| fprintf(ficresp, " Age Prev(%d) N(%d) N ",i,i); |
if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i); |
| fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); |
fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); |
| } |
} |
| fprintf(ficresp, "\n"); |
if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); |
| fprintf(ficresphtm, "\n"); |
fprintf(ficresphtm, "\n"); |
| |
|
| /* Header of frequency table by age */ |
/* Header of frequency table by age */ |
| fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">"); |
fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">"); |
| fprintf(ficresphtmfr,"<th>Age</th> "); |
fprintf(ficresphtmfr,"<th>Age</th> "); |
| for(jk=-1; jk <=nlstate+ndeath; jk++){ |
for(s2=-1; s2 <=nlstate+ndeath; s2++){ |
| for(m=-1; m <=nlstate+ndeath; m++){ |
for(m=-1; m <=nlstate+ndeath; m++){ |
| if(jk!=0 && m!=0) |
if(s2!=0 && m!=0) |
| fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m); |
fprintf(ficresphtmfr,"<th>%d%d</th> ",s2,m); |
| } |
} |
| } |
} |
| fprintf(ficresphtmfr, "\n"); |
fprintf(ficresphtmfr, "\n"); |
|
Line 4577 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4650 Title=%s <br>Datafile=%s Firstpass=%d La
|
| fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage); |
fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage); |
| fprintf(ficlog,"Age %d", iage); |
fprintf(ficlog,"Age %d", iage); |
| } |
} |
| for(jk=1; jk <=nlstate ; jk++){ |
for(s1=1; s1 <=nlstate ; s1++){ |
| for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) |
for(m=-1, pp[s1]=0; m <=nlstate+ndeath ; m++) |
| pp[jk] += freq[jk][m][iage]; |
pp[s1] += freq[s1][m][iage]; |
| } |
} |
| for(jk=1; jk <=nlstate ; jk++){ |
for(s1=1; s1 <=nlstate ; s1++){ |
| for(m=-1, pos=0; m <=0 ; m++) |
for(m=-1, pos=0; m <=0 ; m++) |
| pos += freq[jk][m][iage]; |
pos += freq[s1][m][iage]; |
| if(pp[jk]>=1.e-10){ |
if(pp[s1]>=1.e-10){ |
| if(first==1){ |
if(first==1){ |
| printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
printf(" %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); |
| } |
} |
| fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]); |
| }else{ |
}else{ |
| if(first==1) |
if(first==1) |
| printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
printf(" %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); |
| fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1); |
| } |
} |
| } |
} |
| |
|
| for(jk=1; jk <=nlstate ; jk++){ |
for(s1=1; s1 <=nlstate ; s1++){ |
| /* posprop[jk]=0; */ |
/* posprop[s1]=0; */ |
| for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ |
for(m=0, pp[s1]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ |
| pp[jk] += freq[jk][m][iage]; |
pp[s1] += freq[s1][m][iage]; |
| } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ |
} /* pp[s1] is the total number of transitions starting from state s1 and any ending status until this age */ |
| |
|
| for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ |
for(s1=1,pos=0, pospropta=0.; s1 <=nlstate ; s1++){ |
| pos += pp[jk]; /* pos is the total number of transitions until this age */ |
pos += pp[s1]; /* pos is the total number of transitions until this age */ |
| posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state |
posprop[s1] += prop[s1][iage]; /* prop is the number of transitions from a live state |
| from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
| pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state |
pospropta += prop[s1][iage]; /* prop is the number of transitions from a live state |
| from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
| } |
} |
| for(jk=1; jk <=nlstate ; jk++){ |
|
| |
/* Writing ficresp */ |
| |
if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ |
| |
if( iage <= iagemax){ |
| |
fprintf(ficresp," %d",iage); |
| |
} |
| |
}else if( nj==2){ |
| |
if( iage <= iagemax){ |
| |
fprintf(ficresp," %d",iage); |
| |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| |
} |
| |
} |
| |
for(s1=1; s1 <=nlstate ; s1++){ |
| if(pos>=1.e-5){ |
if(pos>=1.e-5){ |
| if(first==1) |
if(first==1) |
| printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
printf(" %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos); |
| fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos); |
| }else{ |
}else{ |
| if(first==1) |
if(first==1) |
| printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
printf(" %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); |
| fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1); |
| } |
} |
| if( iage <= iagemax){ |
if( iage <= iagemax){ |
| if(pos>=1.e-5){ |
if(pos>=1.e-5){ |
| fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ |
| /* fprintf(ficresp, "%d %d %d %.5f %.0f %.0f",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)],iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); */ |
fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); |
| |
}else if( nj==2){ |
| fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); |
| /*probs[iage][jk][j1]= pp[jk]/pos;*/ |
} |
| /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ |
fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[s1][iage]/pospropta, prop[s1][iage],pospropta); |
| } |
/*probs[iage][s1][j1]= pp[s1]/pos;*/ |
| else{ |
/*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/ |
| fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); |
} else{ |
| fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta); |
if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); |
| |
fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[s1][iage],pospropta); |
| } |
} |
| } |
} |
| pospropt[jk] +=posprop[jk]; |
pospropt[s1] +=posprop[s1]; |
| } /* end loop jk */ |
} /* end loop s1 */ |
| /* pospropt=0.; */ |
/* pospropt=0.; */ |
| for(jk=-1; jk <=nlstate+ndeath; jk++){ |
for(s1=-1; s1 <=nlstate+ndeath; s1++){ |
| for(m=-1; m <=nlstate+ndeath; m++){ |
for(m=-1; m <=nlstate+ndeath; m++){ |
| if(freq[jk][m][iage] !=0 ) { /* minimizing output */ |
if(freq[s1][m][iage] !=0 ) { /* minimizing output */ |
| if(first==1){ |
if(first==1){ |
| printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); |
printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); |
| } |
} |
| /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); */ |
/* printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); */ |
| fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); |
fprintf(ficlog," %d%d=%.0f",s1,m,freq[s1][m][iage]); |
| } |
} |
| if(jk!=0 && m!=0) |
if(s1!=0 && m!=0) |
| fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]); |
fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[s1][m][iage]); |
| } |
} |
| } /* end loop jk */ |
} /* end loop s1 */ |
| posproptt=0.; |
posproptt=0.; |
| for(jk=1; jk <=nlstate; jk++){ |
for(s1=1; s1 <=nlstate; s1++){ |
| posproptt += pospropt[jk]; |
posproptt += pospropt[s1]; |
| } |
} |
| fprintf(ficresphtmfr,"</tr>\n "); |
fprintf(ficresphtmfr,"</tr>\n "); |
| if(iage <= iagemax){ |
fprintf(ficresphtm,"</tr>\n"); |
| fprintf(ficresp,"\n"); |
if((cptcoveff==0 && nj==1)|| nj==2 ) { |
| fprintf(ficresphtm,"</tr>\n"); |
if(iage <= iagemax) |
| |
fprintf(ficresp,"\n"); |
| } |
} |
| if(first==1) |
if(first==1) |
| printf("Others in log...\n"); |
printf("Others in log...\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| } /* end loop age iage */ |
} /* end loop age iage */ |
| |
|
| fprintf(ficresphtm,"<tr><th>Tot</th>"); |
fprintf(ficresphtm,"<tr><th>Tot</th>"); |
| for(jk=1; jk <=nlstate ; jk++){ |
for(s1=1; s1 <=nlstate ; s1++){ |
| if(posproptt < 1.e-5){ |
if(posproptt < 1.e-5){ |
| fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt); |
fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[s1],posproptt); |
| }else{ |
}else{ |
| fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt); |
fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[s1]/posproptt,pospropt[s1],posproptt); |
| } |
} |
| } |
} |
| fprintf(ficresphtm,"</tr>\n"); |
fprintf(ficresphtm,"</tr>\n"); |
|
Line 4687 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4775 Title=%s <br>Datafile=%s Firstpass=%d La
|
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| if(j!=0){ |
if(j!=0){ |
| printf("#Freqsummary: Starting values for combination j1=%d:\n", j1); |
printf("#Freqsummary: Starting values for combination j1=%d:\n", j1); |
| for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,s1=1; i <=nlstate; i++){ |
| for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
| if (k != i) { |
if (k != i) { |
| for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */ |
for(jj=1; jj <=ncovmodel; jj++){ /* For counting s1 */ |
| if(jj==1){ /* Constant case (in fact cste + age) */ |
if(jj==1){ /* Constant case (in fact cste + age) */ |
| if(j1==1){ /* All dummy covariates to zero */ |
if(j1==1){ /* All dummy covariates to zero */ |
| freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */ |
freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */ |
| freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */ |
freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */ |
| printf("%d%d ",i,k); |
printf("%d%d ",i,k); |
| fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
| printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3])); |
printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3])); |
| fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
| pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); |
pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); |
| } |
} |
| }else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */ |
}else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */ |
| for(iage=iagemin; iage <= iagemax+3; iage++){ |
for(iage=iagemin; iage <= iagemax+3; iage++){ |
| x[iage]= (double)iage; |
x[iage]= (double)iage; |
| y[iage]= log(freq[i][k][iage]/freq[i][i][iage]); |
y[iage]= log(freq[i][k][iage]/freq[i][i][iage]); |
| /* printf("i=%d, k=%d, jk=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,jk,j1,jj, iage, y[iage]); */ |
/* printf("i=%d, k=%d, s1=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,s1,j1,jj, iage, y[iage]); */ |
| } |
} |
| linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */ |
linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */ |
| pstart[jk]=b; |
pstart[s1]=b; |
| pstart[jk-1]=a; |
pstart[s1-1]=a; |
| }else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ |
}else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ |
| printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); |
printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); |
| printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); |
printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]); |
| pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])); |
pstart[s1]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])); |
| printf("%d%d ",i,k); |
printf("%d%d ",i,k); |
| fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
| printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4])); |
printf("s1=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",s1,i,k,s1,p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4])); |
| }else{ /* Other cases, like quantitative fixed or varying covariates */ |
}else{ /* Other cases, like quantitative fixed or varying covariates */ |
| ; |
; |
| } |
} |
| /* printf("%12.7f )", param[i][jj][k]); */ |
/* printf("%12.7f )", param[i][jj][k]); */ |
| /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ |
/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ |
| jk++; |
s1++; |
| } /* end jj */ |
} /* end jj */ |
| } /* end k!= i */ |
} /* end k!= i */ |
| } /* end k */ |
} /* end k */ |
| } /* end i, jk */ |
} /* end i, s1 */ |
| } /* end j !=0 */ |
} /* end j !=0 */ |
| } /* end selected combination of covariate j1 */ |
} /* end selected combination of covariate j1 */ |
| if(j==0){ /* We can estimate starting values from the occurences in each case */ |
if(j==0){ /* We can estimate starting values from the occurences in each case */ |
| printf("#Freqsummary: Starting values for the constants:\n"); |
printf("#Freqsummary: Starting values for the constants:\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,s1=1; i <=nlstate; i++){ |
| for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
| if (k != i) { |
if (k != i) { |
| printf("%d%d ",i,k); |
printf("%d%d ",i,k); |
| fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
| for(jj=1; jj <=ncovmodel; jj++){ |
for(jj=1; jj <=ncovmodel; jj++){ |
| pstart[jk]=p[jk]; /* Setting pstart to p values by default */ |
pstart[s1]=p[s1]; /* Setting pstart to p values by default */ |
| if(jj==1){ /* Age has to be done */ |
if(jj==1){ /* Age has to be done */ |
| pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); |
pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]); |
| printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
| fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])); |
| } |
} |
| /* printf("%12.7f )", param[i][jj][k]); */ |
/* printf("%12.7f )", param[i][jj][k]); */ |
| /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ |
/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */ |
| jk++; |
s1++; |
| } |
} |
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
|
Line 4755 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4843 Title=%s <br>Datafile=%s Firstpass=%d La
|
| } |
} |
| printf("#Freqsummary\n"); |
printf("#Freqsummary\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| for(jk=-1; jk <=nlstate+ndeath; jk++){ |
for(s1=-1; s1 <=nlstate+ndeath; s1++){ |
| for(m=-1; m <=nlstate+ndeath; m++){ |
for(s2=-1; s2 <=nlstate+ndeath; s2++){ |
| /* param[i]|j][k]= freq[jk][m][iagemax+3] */ |
/* param[i]|j][k]= freq[s1][s2][iagemax+3] */ |
| printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); |
printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); |
| fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); |
fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); |
| /* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */ |
/* if(freq[s1][s2][iage] !=0 ) { /\* minimizing output *\/ */ |
| /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ |
/* printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ |
| /* fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */ |
/* fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */ |
| /* } */ |
/* } */ |
| } |
} |
| } /* end loop jk */ |
} /* end loop s1 */ |
| |
|
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
|
Line 4899 void prevalence(double ***probs, double
|
Line 4987 void prevalence(double ***probs, double
|
| } else{ |
} else{ |
| if(first==1){ |
if(first==1){ |
| first=0; |
first=0; |
| printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,j1,probs[i][jk][j1]); |
printf("Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); |
| |
fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); |
| |
}else{ |
| |
fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); |
| } |
} |
| } |
} |
| } |
} |
|
Line 6137 void varprob(char optionfilefiname[], do
|
Line 6228 void varprob(char optionfilefiname[], do
|
| fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
| fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
| |
|
| fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov); |
fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. %s</li>\n",optionfilehtmcov,optionfilehtmcov); |
| fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov); |
fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov); |
| fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \ |
fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \ |
| and drawn. It helps understanding how is the covariance between two incidences.\ |
and drawn. It helps understanding how is the covariance between two incidences.\ |
|
Line 6354 To be simple, these graphs help to under
|
Line 6445 To be simple, these graphs help to under
|
| fprintf(ficgp,"\nset parametric;unset label"); |
fprintf(ficgp,"\nset parametric;unset label"); |
| fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); |
fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); |
| fprintf(ficgp,"\nset ter svg size 640, 480"); |
fprintf(ficgp,"\nset ter svg size 640, 480"); |
| fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\ |
fprintf(fichtmcov,"\n<p><br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\ |
| :<a href=\"%s_%d%1d%1d-%1d%1d.svg\"> \ |
:<a href=\"%s_%d%1d%1d-%1d%1d.svg\"> \ |
| %s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\ |
%s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\ |
| subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \ |
subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \ |
|
Line 6365 To be simple, these graphs help to under
|
Line 6456 To be simple, these graphs help to under
|
| fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
| fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
| fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
| mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ |
mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \ |
| mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); |
mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */ |
| }else{ |
}else{ |
| first=0; |
first=0; |
| fprintf(fichtmcov," %d (%.3f),",(int) age, c12); |
fprintf(fichtmcov," %d (%.3f),",(int) age, c12); |
| fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
| fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
| fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
| mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ |
mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \ |
| mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); |
mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); |
| }/* if first */ |
}/* if first */ |
| } /* age mod 5 */ |
} /* age mod 5 */ |
| } /* end loop age */ |
} /* end loop age */ |
|
Line 6665 true period expectancies (those weighted
|
Line 6756 true period expectancies (those weighted
|
| } |
} |
| |
|
| /******************* Gnuplot file **************/ |
/******************* Gnuplot file **************/ |
| void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){ |
void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear){ |
| |
|
| char dirfileres[132],optfileres[132]; |
char dirfileres[132],optfileres[132]; |
| char gplotcondition[132], gplotlabel[132]; |
char gplotcondition[132], gplotlabel[132]; |
|
Line 6675 void printinggnuplot(char fileresu[], ch
|
Line 6766 void printinggnuplot(char fileresu[], ch
|
| int vpopbased; |
int vpopbased; |
| int ioffset; /* variable offset for columns */ |
int ioffset; /* variable offset for columns */ |
| int nres=0; /* Index of resultline */ |
int nres=0; /* Index of resultline */ |
| |
int istart=1; /* For starting graphs in projections */ |
| |
|
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
| /* printf("Problem with file %s",optionfilegnuplot); */ |
/* printf("Problem with file %s",optionfilegnuplot); */ |
|
Line 6773 void printinggnuplot(char fileresu[], ch
|
Line 6865 void printinggnuplot(char fileresu[], ch
|
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1)); |
/* fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1)); */ |
| |
|
| |
fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_")); |
| |
if(cptcoveff ==0){ |
| |
fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3", 2+(cpt-1), cpt ); |
| |
}else{ |
| |
kl=0; |
| |
for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
| |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[k]][lv]; |
| |
kl++; |
| |
/* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ |
| |
/*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ |
| |
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
| |
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ |
| |
if(k==cptcoveff){ |
| |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Observed prevalence in state %d' w l lt 2",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ |
| |
2+cptcoveff*2+3*(cpt-1), cpt ); /* 4 or 6 ?*/ |
| |
}else{ |
| |
fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]); |
| |
kl++; |
| |
} |
| |
} /* end covariate */ |
| |
} /* end if no covariate */ |
| |
|
| if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
| /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
| fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */ |
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */ |
|
Line 7178 set ter svg size 640, 480\nunset log y\n
|
Line 7297 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); |
fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
| set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
| for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ |
|
| |
/* for (i=1; i<= nlstate+1 ; i ++){ /\* nlstate +1 p11 p21 p.1 *\/ */ |
| |
istart=nlstate+1; /* Could be one if by state, but nlstate+1 is w.i projection only */ |
| |
/*istart=1;*/ /* Could be one if by state, but nlstate+1 is w.i projection only */ |
| |
for (i=istart; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ |
| /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
/*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| if(i==1){ |
if(i==istart){ |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); |
| }else{ |
}else{ |
| fprintf(ficgp,",\\\n '' "); |
fprintf(ficgp,",\\\n '' "); |
|
Line 7195 set ter svg size 640, 480\nunset log y\n
|
Line 7318 set ter svg size 640, 480\nunset log y\n
|
| /*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ |
/*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ |
| /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
| fprintf(ficgp," u %d:(", ioffset); |
fprintf(ficgp," u %d:(", ioffset); |
| if(i==nlstate+1) |
if(i==nlstate+1){ |
| fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ", \ |
fprintf(ficgp," $%d/(1.-$%d)):5 t 'pw.%d' with line lc variable ", \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| else |
fprintf(ficgp,",\\\n '' "); |
| |
fprintf(ficgp," u %d:(",ioffset); |
| |
fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \ |
| |
offyear, \ |
| |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| |
}else |
| fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ |
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| }else{ /* more than 2 covariates */ |
}else{ /* more than 2 covariates */ |
|
Line 7230 set ter svg size 640, 480\nunset log y\n
|
Line 7358 set ter svg size 640, 480\nunset log y\n
|
| /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
| /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ |
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ |
| if(i==nlstate+1){ |
if(i==nlstate+1){ |
| fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):5 t 'p.%d' with line lc variable", gplotcondition, \ |
| |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| |
fprintf(ficgp,",\\\n '' "); |
| |
fprintf(ficgp," u %d:(",ioffset); |
| |
fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ |
| |
offyear, \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| |
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/ |
| }else{ |
}else{ |
| fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
|
Line 7451 set ter svg size 640, 480\nunset log y\n
|
Line 7585 set ter svg size 640, 480\nunset log y\n
|
| int mobilavrange, mob; |
int mobilavrange, mob; |
| int iage=0; |
int iage=0; |
| |
|
| double sum=0.; |
double sum=0., sumr=0.; |
| double age; |
double age; |
| double *sumnewp, *sumnewm; |
double *sumnewp, *sumnewm, *sumnewmr; |
| double *agemingood, *agemaxgood; /* Currently identical for all covariates */ |
double *agemingood, *agemaxgood; |
| |
double *agemingoodr, *agemaxgoodr; |
| |
|
| |
|
| /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ |
/* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ |
|
Line 7462 set ter svg size 640, 480\nunset log y\n
|
Line 7597 set ter svg size 640, 480\nunset log y\n
|
| |
|
| sumnewp = vector(1,ncovcombmax); |
sumnewp = vector(1,ncovcombmax); |
| sumnewm = vector(1,ncovcombmax); |
sumnewm = vector(1,ncovcombmax); |
| |
sumnewmr = vector(1,ncovcombmax); |
| agemingood = vector(1,ncovcombmax); |
agemingood = vector(1,ncovcombmax); |
| |
agemingoodr = vector(1,ncovcombmax); |
| agemaxgood = vector(1,ncovcombmax); |
agemaxgood = vector(1,ncovcombmax); |
| |
agemaxgoodr = vector(1,ncovcombmax); |
| |
|
| for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
| sumnewm[cptcod]=0.; |
sumnewm[cptcod]=0.; sumnewmr[cptcod]=0.; |
| sumnewp[cptcod]=0.; |
sumnewp[cptcod]=0.; |
| agemingood[cptcod]=0; |
agemingood[cptcod]=0, agemingoodr[cptcod]=0; |
| agemaxgood[cptcod]=0; |
agemaxgood[cptcod]=0, agemaxgoodr[cptcod]=0; |
| } |
} |
| if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */ |
if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */ |
| |
|
| if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
if(mobilav==-1 || mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
| if(mobilav==1) mobilavrange=5; /* default */ |
if(mobilav==1 || mobilav==-1) mobilavrange=5; /* default */ |
| else mobilavrange=mobilav; |
else mobilavrange=mobilav; |
| for (age=bage; age<=fage; age++) |
for (age=bage; age<=fage; age++) |
| for (i=1; i<=nlstate;i++) |
for (i=1; i<=nlstate;i++) |
|
Line 7486 set ter svg size 640, 480\nunset log y\n
|
Line 7624 set ter svg size 640, 480\nunset log y\n
|
| */ |
*/ |
| for (mob=3;mob <=mobilavrange;mob=mob+2){ |
for (mob=3;mob <=mobilavrange;mob=mob+2){ |
| for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
| for (i=1; i<=nlstate;i++){ |
for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
| for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
sumnewm[cptcod]=0.; |
| |
for (i=1; i<=nlstate;i++){ |
| mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
| for (cpt=1;cpt<=(mob-1)/2;cpt++){ |
for (cpt=1;cpt<=(mob-1)/2;cpt++){ |
| mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; |
mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; |
| mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; |
mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; |
| } |
} |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; |
| } |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| } |
} /* end i */ |
| |
if(sumnewm[cptcod] >1.e-3) mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/sumnewm[cptcod]; /* Rescaling to sum one */ |
| |
} /* end cptcod */ |
| }/* end age */ |
}/* end age */ |
| }/* end mob */ |
}/* end mob */ |
| }else |
}else{ |
| |
printf("Error internal in movingaverage, mobilav=%d.\n",mobilav); |
| return -1; |
return -1; |
| for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
} |
| |
|
| |
for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ /* for each combination */ |
| /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */ |
/* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */ |
| if(invalidvarcomb[cptcod]){ |
if(invalidvarcomb[cptcod]){ |
| printf("\nCombination (%d) ignored because no cases \n",cptcod); |
printf("\nCombination (%d) ignored because no cases \n",cptcod); |
| continue; |
continue; |
| } |
} |
| |
|
| agemingood[cptcod]=fage-(mob-1)/2; |
for (age=fage-(mob-1)/2; age>=bage+(mob-1)/2; age--){ /*looking for the youngest and oldest good age */ |
| for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */ |
|
| sumnewm[cptcod]=0.; |
sumnewm[cptcod]=0.; |
| |
sumnewmr[cptcod]=0.; |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| |
sumnewmr[cptcod]+=probs[(int)age][i][cptcod]; |
| |
} |
| |
if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good without smoothing */ |
| |
agemingoodr[cptcod]=age; |
| } |
} |
| if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
| agemingood[cptcod]=age; |
agemingood[cptcod]=age; |
| }else{ /* bad */ |
} |
| for (i=1; i<=nlstate;i++){ |
} /* age */ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ /*looking for the youngest and oldest good age */ |
| } /* i */ |
|
| } /* end bad */ |
|
| }/* age */ |
|
| sum=0.; |
|
| for (i=1; i<=nlstate;i++){ |
|
| sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
|
| } |
|
| if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
|
| printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod); |
|
| /* for (i=1; i<=nlstate;i++){ */ |
|
| /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ |
|
| /* } /\* i *\/ */ |
|
| } /* end bad */ |
|
| /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
|
| /* From youngest, finding the oldest wrong */ |
|
| agemaxgood[cptcod]=bage+(mob-1)/2; |
|
| for (age=bage+(mob-1)/2; age<=fage; age++){ |
|
| sumnewm[cptcod]=0.; |
sumnewm[cptcod]=0.; |
| |
sumnewmr[cptcod]=0.; |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| |
sumnewmr[cptcod]+=probs[(int)age][i][cptcod]; |
| |
} |
| |
if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good without smoothing */ |
| |
agemaxgoodr[cptcod]=age; |
| } |
} |
| if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
| agemaxgood[cptcod]=age; |
agemaxgood[cptcod]=age; |
| }else{ /* bad */ |
} |
| for (i=1; i<=nlstate;i++){ |
} /* age */ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
/* Thus we have agemingood and agemaxgood as well as goodr for raw (preobs) */ |
| } /* i */ |
/* but they will change */ |
| |
for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, filling up to the youngest */ |
| |
sumnewm[cptcod]=0.; |
| |
sumnewmr[cptcod]=0.; |
| |
for (i=1; i<=nlstate;i++){ |
| |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| |
sumnewmr[cptcod]+=probs[(int)age][i][cptcod]; |
| |
} |
| |
if(mobilav==-1){ /* Forcing raw ages if good else agemingood */ |
| |
if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good without smoothing */ |
| |
agemaxgoodr[cptcod]=age; /* age min */ |
| |
for (i=1; i<=nlstate;i++) |
| |
mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
| |
}else{ /* bad we change the value with the values of good ages */ |
| |
for (i=1; i<=nlstate;i++){ |
| |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgoodr[cptcod]][i][cptcod]; |
| |
} /* i */ |
| |
} /* end bad */ |
| |
}else{ |
| |
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
| |
agemaxgood[cptcod]=age; |
| |
}else{ /* bad we change the value with the values of good ages */ |
| |
for (i=1; i<=nlstate;i++){ |
| |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
| |
} /* i */ |
| |
} /* end bad */ |
| |
}/* end else */ |
| |
sum=0.;sumr=0.; |
| |
for (i=1; i<=nlstate;i++){ |
| |
sum+=mobaverage[(int)age][i][cptcod]; |
| |
sumr+=probs[(int)age][i][cptcod]; |
| |
} |
| |
if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
| |
printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, bage); |
| |
} /* end bad */ |
| |
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
| |
if(fabs(sumr - 1.) > 1.e-3) { /* bad */ |
| |
printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, bage); |
| } /* end bad */ |
} /* end bad */ |
| }/* age */ |
}/* age */ |
| sum=0.; |
|
| for (i=1; i<=nlstate;i++){ |
for (age=bage+(mob-1)/2; age<=fage; age++){/* From youngest, finding the oldest wrong */ |
| sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
sumnewm[cptcod]=0.; |
| } |
sumnewmr[cptcod]=0.; |
| if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
for (i=1; i<=nlstate;i++){ |
| printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod); |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| /* for (i=1; i<=nlstate;i++){ */ |
sumnewmr[cptcod]+=probs[(int)age][i][cptcod]; |
| /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ |
} |
| /* } /\* i *\/ */ |
if(mobilav==-1){ /* Forcing raw ages if good else agemingood */ |
| } /* end bad */ |
if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good */ |
| |
agemingoodr[cptcod]=age; |
| |
for (i=1; i<=nlstate;i++) |
| |
mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
| |
}else{ /* bad we change the value with the values of good ages */ |
| |
for (i=1; i<=nlstate;i++){ |
| |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingoodr[cptcod]][i][cptcod]; |
| |
} /* i */ |
| |
} /* end bad */ |
| |
}else{ |
| |
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
| |
agemingood[cptcod]=age; |
| |
}else{ /* bad */ |
| |
for (i=1; i<=nlstate;i++){ |
| |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
| |
} /* i */ |
| |
} /* end bad */ |
| |
}/* end else */ |
| |
sum=0.;sumr=0.; |
| |
for (i=1; i<=nlstate;i++){ |
| |
sum+=mobaverage[(int)age][i][cptcod]; |
| |
sumr+=mobaverage[(int)age][i][cptcod]; |
| |
} |
| |
if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
| |
printf("Moving average B1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you decrease fage=%d?\n",cptcod, sum, (int) age, fage); |
| |
} /* end bad */ |
| |
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
| |
if(fabs(sumr - 1.) > 1.e-3) { /* bad */ |
| |
printf("Moving average B2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase fage=%d\n",cptcod,sumr, (int)age, fage); |
| |
} /* end bad */ |
| |
}/* age */ |
| |
|
| |
|
| for (age=bage; age<=fage; age++){ |
for (age=bage; age<=fage; age++){ |
| /* printf("%d %d ", cptcod, (int)age); */ |
/* printf("%d %d ", cptcod, (int)age); */ |
|
Line 7571 set ter svg size 640, 480\nunset log y\n
|
Line 7775 set ter svg size 640, 480\nunset log y\n
|
| } |
} |
| /* printf("\n"); */ |
/* printf("\n"); */ |
| /* } */ |
/* } */ |
| |
|
| /* brutal averaging */ |
/* brutal averaging */ |
| for (i=1; i<=nlstate;i++){ |
/* for (i=1; i<=nlstate;i++){ */ |
| for (age=1; age<=bage; age++){ |
/* for (age=1; age<=bage; age++){ */ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
/* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ |
| /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ |
/* /\* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); *\/ */ |
| } |
/* } */ |
| for (age=fage; age<=AGESUP; age++){ |
/* for (age=fage; age<=AGESUP; age++){ */ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
/* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; */ |
| /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ |
/* /\* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); *\/ */ |
| } |
/* } */ |
| } /* end i status */ |
/* } /\* end i status *\/ */ |
| for (i=nlstate+1; i<=nlstate+ndeath;i++){ |
/* for (i=nlstate+1; i<=nlstate+ndeath;i++){ */ |
| for (age=1; age<=AGESUP; age++){ |
/* for (age=1; age<=AGESUP; age++){ */ |
| /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/ |
/* /\*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*\/ */ |
| mobaverage[(int)age][i][cptcod]=0.; |
/* mobaverage[(int)age][i][cptcod]=0.; */ |
| } |
/* } */ |
| } |
/* } */ |
| }/* end cptcod */ |
}/* end cptcod */ |
| free_vector(sumnewm,1, ncovcombmax); |
free_vector(agemaxgoodr,1, ncovcombmax); |
| free_vector(sumnewp,1, ncovcombmax); |
|
| free_vector(agemaxgood,1, ncovcombmax); |
free_vector(agemaxgood,1, ncovcombmax); |
| free_vector(agemingood,1, ncovcombmax); |
free_vector(agemingood,1, ncovcombmax); |
| |
free_vector(agemingoodr,1, ncovcombmax); |
| |
free_vector(sumnewmr,1, ncovcombmax); |
| |
free_vector(sumnewm,1, ncovcombmax); |
| |
free_vector(sumnewp,1, ncovcombmax); |
| return 0; |
return 0; |
| }/* End movingaverage */ |
}/* End movingaverage */ |
| |
|
|
Line 7699 set ter svg size 640, 480\nunset log y\n
|
Line 7907 set ter svg size 640, 480\nunset log y\n
|
| for(j=1; j<=nlstate+ndeath;j++) { |
for(j=1; j<=nlstate+ndeath;j++) { |
| ppij=0.; |
ppij=0.; |
| for(i=1; i<=nlstate;i++) { |
for(i=1; i<=nlstate;i++) { |
| if (mobilav==1) |
/* if (mobilav>=1) */ |
| ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k]; |
ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k]; |
| else { |
/* else { */ /* even if mobilav==-1 we use mobaverage */ |
| ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; |
/* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ |
| } |
/* } */ |
| if (h*hstepm/YEARM*stepm== yearp) { |
if (h*hstepm/YEARM*stepm== yearp) { |
| fprintf(ficresf," %.3f", p3mat[i][j][h]); |
fprintf(ficresf," %.3f", p3mat[i][j][h]); |
| } |
} |
|
Line 7715 set ter svg size 640, 480\nunset log y\n
|
Line 7923 set ter svg size 640, 480\nunset log y\n
|
| } /* end h */ |
} /* end h */ |
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| } /* end agec */ |
} /* end agec */ |
| |
/* diffyear=(int) anproj1+yearp-ageminpar-1; */ |
| |
/*printf("Prevforecast %d+%d-%d=diffyear=%d\n",(int) anproj1, (int)yearp,(int)ageminpar,(int) anproj1-(int)ageminpar);*/ |
| } /* end yearp */ |
} /* end yearp */ |
| } /* end k */ |
} /* end k */ |
| |
|
|
Line 9696 int back_prevalence_limit(double *p, dou
|
Line 9906 int back_prevalence_limit(double *p, dou
|
| }else{ |
}else{ |
| /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */ |
/* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */ |
| bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k,nres); |
bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k,nres); |
| |
/* printf("TOTOT\n"); */ |
| |
/* exit(1); */ |
| } |
} |
| fprintf(ficresplb,"%.0f ",age ); |
fprintf(ficresplb,"%.0f ",age ); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
|
Line 9854 int hPijx(double *p, int bage, int fage)
|
Line 10066 int hPijx(double *p, int bage, int fage)
|
| |
|
| /* nhstepm=nhstepm*YEARM; aff par mois*/ |
/* nhstepm=nhstepm*YEARM; aff par mois*/ |
| |
|
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); /* We can't have it at an upper level because of nhstepm */ |
| |
/* and memory limitations if stepm is small */ |
| |
|
| /* oldm=oldms;savm=savms; */ |
/* oldm=oldms;savm=savms; */ |
| /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */ |
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */ |
| hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k); |
hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k); |
|
Line 11417 Please run with mle=-1 to get a correct
|
Line 11631 Please run with mle=-1 to get a correct
|
| This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
| Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
| }else{ |
}else{ |
| printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p); |
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)ageminpar); |
| } |
} |
| printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \ |
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \ |
| model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \ |
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \ |
|
Line 11479 Please run with mle=-1 to get a correct
|
Line 11693 Please run with mle=-1 to get a correct
|
| printf(" Error in movingaverage mobilav=%d\n",mobilav); |
printf(" Error in movingaverage mobilav=%d\n",mobilav); |
| } |
} |
| } |
} |
| /* /\* Prevalence for each covariates in probs[age][status][cov] *\/ */ |
/* else if(mobilavproj==-1){ /\* Forcing raw observed prevalences *\/ */ |
| /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */ |
/* for(i=1;i<=AGESUP;i++) */ |
| |
/* for(j=1;j<=nlstate;j++) */ |
| |
/* for(k=1;k<=ncovcombmax;k++) */ |
| |
/* mobaverages[i][j][k]=probs[i][j][k]; */ |
| |
/* /\* /\\* Prevalence for each covariates in probs[age][status][cov] *\\/ *\/ */ |
| |
/* /\* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); *\/ */ |
| |
/* } */ |
| else if (mobilavproj !=0) { |
else if (mobilavproj !=0) { |
| printf("Movingaveraging projected observed prevalence\n"); |
printf("Movingaveraging projected observed prevalence\n"); |
| fprintf(ficlog,"Movingaveraging projected observed prevalence\n"); |
fprintf(ficlog,"Movingaveraging projected observed prevalence\n"); |
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src