version 1.141, 2014/01/26 02:42:01
|
version 1.143, 2014/01/26 09:45:38
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.143 2014/01/26 09:45:38 brouard |
|
Summary: Version 0.98nR (to be improved, but gives same optimization results as 0.98k. Nice, promising |
|
|
|
* imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... |
|
(Module): Version 0.98nR Running ok, but output format still only works for three covariates. |
|
|
|
Revision 1.142 2014/01/26 03:57:36 brouard |
|
Summary: gnuplot changed plot w l 1 has to be changed to plot w l lt 2 |
|
|
|
* imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... |
|
|
Revision 1.141 2014/01/26 02:42:01 brouard |
Revision 1.141 2014/01/26 02:42:01 brouard |
* imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... |
* imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... |
|
|
Line 443 extern int errno;
|
Line 454 extern int errno;
|
/* $Id$ */ |
/* $Id$ */ |
/* $State$ */ |
/* $State$ */ |
|
|
char version[]="Imach version 0.98n, January 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Sicentific Research 25293121)"; |
char version[]="Imach version 0.98nR2, January 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)"; |
char fullversion[]="$Revision$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
char strstart[80]; |
char strstart[80]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
Line 577 int **codtab; /**< codtab=imatrix(1,100,
|
Line 588 int **codtab; /**< codtab=imatrix(1,100,
|
int **Tvard, *Tprod, cptcovprod, *Tvaraff; |
int **Tvard, *Tprod, cptcovprod, *Tvaraff; |
double *lsurv, *lpop, *tpop; |
double *lsurv, *lpop, *tpop; |
|
|
double ftol=FTOL; /* Tolerance for computing Max Likelihood */ |
double ftol=FTOL; /**< Tolerance for computing Max Likelihood */ |
double ftolhess; /* Tolerance for computing hessian */ |
double ftolhess; /**< Tolerance for computing hessian */ |
|
|
/**************** split *************************/ |
/**************** split *************************/ |
static int split( char *path, char *dirc, char *name, char *ext, char *finame ) |
static int split( char *path, char *dirc, char *name, char *ext, char *finame ) |
Line 1295 double **prevalim(double **prlim, int nl
|
Line 1306 double **prevalim(double **prlim, int nl
|
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ |
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ |
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ |
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ |
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ |
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ |
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); |
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */ |
|
|
savm=oldm; |
savm=oldm; |
oldm=newm; |
oldm=newm; |
Line 2198 void freqsummary(char fileres[], int ia
|
Line 2209 void freqsummary(char fileres[], int ia
|
|
|
first=1; |
first=1; |
|
|
for(k1=1; k1<=j;k1++){ |
for(k1=1; k1<=j;k1++){ /* Loop on covariates */ |
for(i1=1; i1<=ncodemax[k1];i1++){ |
for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */ |
j1++; |
j1++; |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
scanf("%d", i);*/ |
scanf("%d", i);*/ |
Line 2207 void freqsummary(char fileres[], int ia
|
Line 2218 void freqsummary(char fileres[], int ia
|
for (jk=-5; jk<=nlstate+ndeath; jk++) |
for (jk=-5; jk<=nlstate+ndeath; jk++) |
for(m=iagemin; m <= iagemax+3; m++) |
for(m=iagemin; m <= iagemax+3; m++) |
freq[i][jk][m]=0; |
freq[i][jk][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++) |
prop[i][m]=0; |
prop[i][m]=0; |
|
|
dateintsum=0; |
dateintsum=0; |
k2cpt=0; |
k2cpt=0; |
Line 2248 void freqsummary(char fileres[], int ia
|
Line 2259 void freqsummary(char fileres[], int ia
|
fprintf(ficresp, "\n#********** Variable "); |
fprintf(ficresp, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
fprintf(ficresp, "**********\n#"); |
fprintf(ficresp, "**********\n#"); |
|
fprintf(ficlog, "\n#********** Variable "); |
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
|
fprintf(ficlog, "**********\n#"); |
} |
} |
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); |
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); |
Line 3225 void varevsij(char optionfilefiname[], d
|
Line 3239 void varevsij(char optionfilefiname[], d
|
/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l 1 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 2 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l lt 3 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l 2 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 3 ",subdirf(fileresprobmorprev)); |
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); |
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); |
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); |
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); |
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit); |
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit); |
Line 3822 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 3836 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
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\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 1,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
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\"95\%% CI\" w l 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"95\%% CI\" w l lt 2,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
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 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1)); |
fprintf(ficgp,"\" t\"\" w l lt 2,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 3",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1)); |
} |
} |
} |
} |
/*2 eme*/ |
/*2 eme*/ |
Line 3855 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 3869 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); |
if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
} |
} |
fprintf(ficgp,"\" t\"\" w l 0,"); |
fprintf(ficgp,"\" t\"\" w l lt 1,"); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
for (j=1; j<= nlstate+1 ; j ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); |
if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
} |
} |
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0"); |
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 1"); |
else fprintf(ficgp,"\" t\"\" w l 0,"); |
else fprintf(ficgp,"\" t\"\" w l lt 1,"); |
} |
} |
} |
} |
|
|
Line 3949 plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
|
Line 3963 plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
|
fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
ij=1;/* To be checked else nbcode[0][0] wrong */ |
ij=1;/* To be checked else nbcode[0][0] wrong */ |
for(j=3; j <=ncovmodel; j++) { |
for(j=3; j <=ncovmodel; j++) { |
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { |
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */ |
fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); |
fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); |
ij++; |
ij++; |
} |
} |
Line 5491 run imach with mle=-1 to get a correct t
|
Line 5505 run imach with mle=-1 to get a correct t
|
m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
|
|
for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ |
for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ |
for(i=1; i <=(m/pow(2,k));i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ |
for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ |
for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */ |
for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/ |
for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ |
for(cpt=1; cpt <=pow(2,k-1); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ |
h++; |
h++; |
if (h>m) |
if (h>m) |
h=1; |
h=1; |
|
/**< codtab(h,k) k |
|
* h 1 2 3 4 |
|
*______________________________ |
|
* 1 i=1 1 i=1 1 i=1 1 i=1 1 |
|
* 2 2 1 1 1 |
|
* 3 i=2 1 2 1 1 |
|
* 4 2 2 1 1 |
|
* 5 i=3 1 i=2 1 2 1 |
|
* 6 2 1 2 1 |
|
* 7 i=4 1 2 2 1 |
|
* 8 2 2 2 1 |
|
* 9 i=5 1 i=3 1 i=2 1 1 |
|
* 10 2 1 1 1 |
|
* 11 i=6 1 2 1 1 |
|
* 12 2 2 1 1 |
|
* 13 i=7 1 i=4 1 2 1 |
|
* 14 2 1 2 1 |
|
* 15 i=8 1 2 2 1 |
|
* 16 2 2 2 1 |
|
*/ |
codtab[h][k]=j; |
codtab[h][k]=j; |
codtab[h][Tvar[k]]=j; |
codtab[h][Tvar[k]]=j; |
printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); |
printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); |
Line 6076 Interval (in months) between two waves:
|
Line 6110 Interval (in months) between two waves:
|
agebase=ageminpar; |
agebase=ageminpar; |
agelim=agemaxpar; |
agelim=agemaxpar; |
ftolpl=1.e-10; |
ftolpl=1.e-10; |
i1=cptcoveff; |
i1=pow(2,cptcoveff); |
if (cptcovn < 1){i1=1;} |
if (cptcovn < 1){i1=1;} |
|
|
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
//for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
k=k+1; |
k=k+1; |
/* to clean */ |
/* to clean */ |
printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode); |
//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]); |
fprintf(ficrespl,"\n#******"); |
fprintf(ficrespl,"\n#******"); |
printf("\n#******"); |
printf("\n#******"); |
fprintf(ficlog,"\n#******"); |
fprintf(ficlog,"\n#******"); |
Line 6104 Interval (in months) between two waves:
|
Line 6138 Interval (in months) between two waves:
|
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
fprintf(ficrespl," %.5f", prlim[i][i]); |
fprintf(ficrespl," %.5f", prlim[i][i]); |
fprintf(ficrespl,"\n"); |
fprintf(ficrespl,"\n"); |
} |
} /* Age */ |
} |
/* was end of cptcod */ |
} |
} /* cptcov */ |
fclose(ficrespl); |
fclose(ficrespl); |
|
|
/*------------- h Pij x at various ages ------------*/ |
/*------------- h Pij x at various ages ------------*/ |