--- imach/src/imach.c 2014/01/26 02:42:01 1.141
+++ imach/src/imach.c 2014/01/26 09:45:38 1.143
@@ -1,6 +1,17 @@
-/* $Id: imach.c,v 1.141 2014/01/26 02:42:01 brouard Exp $
+/* $Id: imach.c,v 1.143 2014/01/26 09:45:38 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ 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
* imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested...
@@ -440,11 +451,11 @@ extern int errno;
#define ODIRSEPARATOR '/'
#endif
-/* $Id: imach.c,v 1.141 2014/01/26 02:42:01 brouard Exp $ */
+/* $Id: imach.c,v 1.143 2014/01/26 09:45:38 brouard Exp $ */
/* $State: Exp $ */
-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 fullversion[]="$Revision: 1.141 $ $Date: 2014/01/26 02:42:01 $";
+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: 1.143 $ $Date: 2014/01/26 09:45:38 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -577,8 +588,8 @@ int **codtab; /**< codtab=imatrix(1,100,
int **Tvard, *Tprod, cptcovprod, *Tvaraff;
double *lsurv, *lpop, *tpop;
-double ftol=FTOL; /* Tolerance for computing Max Likelihood */
-double ftolhess; /* Tolerance for computing hessian */
+double ftol=FTOL; /**< Tolerance for computing Max Likelihood */
+double ftolhess; /**< Tolerance for computing hessian */
/**************** split *************************/
static int split( char *path, char *dirc, char *name, char *ext, char *finame )
@@ -1295,7 +1306,7 @@ double **prevalim(double **prlim, int nl
/*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 \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;
oldm=newm;
@@ -2198,8 +2209,8 @@ void freqsummary(char fileres[], int ia
first=1;
- for(k1=1; k1<=j;k1++){
- for(i1=1; i1<=ncodemax[k1];i1++){
+ for(k1=1; k1<=j;k1++){ /* Loop on covariates */
+ for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */
j1++;
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/
@@ -2207,10 +2218,10 @@ void freqsummary(char fileres[], int ia
for (jk=-5; jk<=nlstate+ndeath; jk++)
for(m=iagemin; m <= iagemax+3; m++)
freq[i][jk][m]=0;
-
- for (i=1; i<=nlstate; i++)
- for(m=iagemin; m <= iagemax+3; m++)
- prop[i][m]=0;
+
+ for (i=1; i<=nlstate; i++)
+ for(m=iagemin; m <= iagemax+3; m++)
+ prop[i][m]=0;
dateintsum=0;
k2cpt=0;
@@ -2248,6 +2259,9 @@ void freqsummary(char fileres[], int ia
fprintf(ficresp, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
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++)
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
@@ -3225,9 +3239,9 @@ 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 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 plot \"%s\" u 1:($3) not w l 1 ",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)) not w l 2 ",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 lt 3 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 3 ",subdirf(fileresprobmorprev));
fprintf(fichtm,"\n
File (multiple files are possible if covariates are present): %s\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
fprintf(fichtm,"\n
Probability is computed over estepm=%d months.
![](\"%s%s.png\")
\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
/* fprintf(fichtm,"\n
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
![](\"varmuptjgr%s%s.png\")
\n", stepm,YEARM,digitp,digit);
@@ -3822,17 +3836,17 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
if (i==cpt) 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 ++) {
if (i==cpt) 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 ++) {
if (i==cpt) 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*/
@@ -3855,14 +3869,14 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
if (j==i) 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);
for (j=1; j<= nlstate+1 ; j ++) {
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
}
- if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
- else fprintf(ficgp,"\" t\"\" w l 0,");
+ if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 1");
+ else fprintf(ficgp,"\" t\"\" w l lt 1,");
}
}
@@ -3949,7 +3963,7 @@ plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
ij=1;/* To be checked else nbcode[0][0] wrong */
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]]]);
ij++;
}
@@ -5491,12 +5505,32 @@ run imach with mle=-1 to get a correct t
m=pow(2,cptcoveff);
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(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */
- 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(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 ncodemax=2*/
+ 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++;
if (h>m)
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][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]]);
@@ -6076,14 +6110,14 @@ Interval (in months) between two waves:
agebase=ageminpar;
agelim=agemaxpar;
ftolpl=1.e-10;
- i1=cptcoveff;
+ i1=pow(2,cptcoveff);
if (cptcovn < 1){i1=1;}
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;
/* 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#******");
printf("\n#******");
fprintf(ficlog,"\n#******");
@@ -6104,9 +6138,9 @@ Interval (in months) between two waves:
for(i=1; i<=nlstate;i++)
fprintf(ficrespl," %.5f", prlim[i][i]);
fprintf(ficrespl,"\n");
- }
- }
- }
+ } /* Age */
+ /* was end of cptcod */
+ } /* cptcov */
fclose(ficrespl);
/*------------- h Pij x at various ages ------------*/