/* $Id$
$State$
$Log$
+ Revision 1.201 2015/09/15 17:34:58 brouard
+ Summary: 0.98r0
+
+ - Some new graphs like suvival functions
+ - Some bugs fixed like model=1+age+V2.
+
Revision 1.200 2015/09/09 16:53:55 brouard
Summary: Big bug thanks to Flavia
/* #define DEBUG */
/* #define DEBUGBRENT */
+#define DEBUGLINMIN
#define POWELL /* Instead of NLOPT */
#define POWELLF1F3 /* Skip test */
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
int outcmd=0;
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
-char fileresu[FILENAMELENGTH]; /* Without r in front */
+char fileresu[FILENAMELENGTH]; /* fileres without r in front */
char filelog[FILENAMELENGTH]; /* Log file */
char filerest[FILENAMELENGTH];
char fileregp[FILENAMELENGTH];
nrfunc=func;
for (j=1;j<=n;j++) {
pcom[j]=p[j];
- xicom[j]=xi[j];
+ xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */
}
/* axs=0.0; */
#ifdef DEBUGLINMIN
printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb);
+ fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb);
#endif
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/
/* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */
#endif
#ifdef DEBUGLINMIN
printf("linmin end ");
+ fprintf(ficlog,"linmin end ");
#endif
for (j=1;j<=n;j++) {
/* printf(" before xi[%d]=%12.8f", j,xi[j]); */
/* printf("\n"); */
#ifdef DEBUGLINMIN
printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
+ fprintf(ficlog,"Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
for (j=1;j<=n;j++) {
- printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]);
- if(j % ncovmodel == 0)
+ printf(" xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
+ fprintf(ficlog," xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
+ if(j % ncovmodel == 0){
printf("\n");
+ fprintf(ficlog,"\n");
+ }
}
#endif
free_vector(xicom,1,n);
xits=vector(1,n);
*fret=(*func)(p);
for (j=1;j<=n;j++) pt[j]=p[j];
- rcurr_time = time(NULL);
+ rcurr_time = time(NULL);
for (*iter=1;;++(*iter)) {
fp=(*fret); /* From former iteration or initial value */
ibig=0;
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */
t= t- del*SQR(fp-fptt);
#endif
- directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */
+ directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */
#ifdef DEBUG
printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
if (t < 0.0) { /* Then we use it for new direction */
#else
if (directest*t < 0.0) { /* Contradiction between both tests */
- printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
+ printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
- fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
+ fprintf(ficlog,"directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
}
if (directest < 0.0) { /* Then we use it for new direction */
#ifdef DEBUGLINMIN
printf("Before linmin in direction P%d-P0\n",n);
for (j=1;j<=n;j++) {
- printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0)
+ printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
printf("\n");
+ fprintf(ficlog,"\n");
+ }
}
#endif
linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
#ifdef DEBUGLINMIN
for (j=1;j<=n;j++) {
printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0)
+ fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
printf("\n");
+ fprintf(ficlog,"\n");
+ }
}
#endif
for (j=1;j<=n;j++) {
/* double **matprod2(); */ /* test */
double **out, cov[NCOVMAX+1], **pmij();
double **newm;
- double agefin, delaymax=50 ; /* Max number of years to converge */
+ double agefin, delaymax=100 ; /* Max number of years to converge */
+ long int ncvyear=0, ncvloop=0;
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
cov[1]=1.;
/* Even if hstepm = 1, at least one multiplication by the unit matrix */
+ /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
+ ncvloop++;
newm=savm;
/* Covariates have to be included here again */
cov[2]=agefin;
sumnew=0;
for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
prlim[i][j]= newm[i][j]/(1-sumnew);
- /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/
max=FMAX(max,prlim[i][j]);
min=FMIN(min,prlim[i][j]);
+ /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */
}
maxmin=max-min;
maxmax=FMAX(maxmax,maxmin);
} /* j loop */
if(maxmax < ftolpl){
+ /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age-(int)agefin); */
return prlim;
}
} /* age loop */
+ printf("Warning: the stable prevalence did not converge with the required precision ftolpl=6*10^5*ftol=%g. \n\
+Earliest age to start was %d-%d=%d, ncvloop=%ld, ncvyear=%d\n\
+Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin);
return prlim; /* should not reach here */
}
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
if(globpr){
- fprintf(ficresilk,"%9ld %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
+ fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f\
%11.6f %11.6f %11.6f ", \
- num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
+ num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
for(k=1,llt=0.,l=0.; k<=nlstate; k++){
llt +=ll[k]*gipmx/gsw;
if(*globpri !=0){ /* Just counts and sums, no printings */
strcpy(fileresilk,"ILK_");
- strcat(fileresilk,fileres);
+ strcat(fileresilk,fileresu);
if((ficresilk=fopen(fileresilk,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresilk);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
}
fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
- fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight 2wlli out sav ");
+ fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight 2wlli out sav ");
/* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
for(k=1; k<=nlstate; k++)
fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
*fretone=(*funcone)(p);
if(*globpri !=0){
fclose(ficresilk);
- fprintf(fichtm,"\n<br>File of contributions to the likelihood (if mle=1): <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
- fflush(fichtm);
+ fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with initial parameters and mle >= 1. You should at least run with mle >= 1 and starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
+ fprintf(fichtm,"<br>- The first 3 individuals are drawn with lines. The function drawn is -2Log(L) in log scale: <a href=\"%s.png\">%s.png</a><br> \
+<img src=\"%s.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ fflush(fichtm);
}
return;
}
/*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
strcat(fileresprobmorprev,digit); /* Tvar to be done */
strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
- strcat(fileresprobmorprev,fileres);
+ strcat(fileresprobmorprev,fileresu);
if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprobmorprev);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
}
strcpy(fileresprobcov,"PROBCOV_");
- strcat(fileresprobcov,fileres);
+ strcat(fileresprobcov,fileresu);
if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprobcov);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
}
strcpy(fileresprobcor,"PROBCOR_");
- strcat(fileresprobcor,fileres);
+ strcat(fileresprobcor,fileresu);
if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprobcor);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
<img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Pij */
- fprintf(fichtm,"<br>\n- Pij or
Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+ fprintf(fichtm,"<br>\n- Pij or
conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Quasi-incidences */
fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
/*#endif */
m=pow(2,cptcoveff);
+ /* Contribution to likelihood */
+ /* Plot the probability implied in the likelihood */
+ fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
+ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
+ /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
+ fprintf(ficgp,"\nset ter png size 640, 480");
+/* good for mle=4 plot by number of matrix products.
+ replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */
+/* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */
+ /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
+ fprintf(ficgp,"\nset out \"%s.png\";",subdirf2(optionfilefiname,"ILK_"));
+ fprintf(ficgp,"\nplot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk));
+ fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk));
+ fprintf(ficgp,"\nset out\n");
+ /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
+
strcpy(dirfileres,optionfilefiname);
strcpy(optfileres,"vpl");
/* 1eme*/
} /* end covariate */
/* Survival functions (period) from state i in state j by final state j */
- for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
+ for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */
k=3;
fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
} /* end cpt state*/
} /* end covariate */
- /* CV preval stable (period) */
- for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
+ /* CV preval stable (period) for each covariate */
+ for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
k=3;
fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
}
-int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){
+ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl){
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
int i, j, k, i1 ;
- double ftolpl = 1.e-10;
+ /* double ftolpl = 1.e-10; */
double age, agebase, agelim;
- strcpy(filerespl,"PL_");
- strcat(filerespl,fileresu);
- if((ficrespl=fopen(filerespl,"w"))==NULL) {
- printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
- fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
- }
- printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
- fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
- pstamp(ficrespl);
- fprintf(ficrespl,"# Period (stable) prevalence \n");
- fprintf(ficrespl,"#Age ");
- for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
- fprintf(ficrespl,"\n");
+ strcpy(filerespl,"PL_");
+ strcat(filerespl,fileresu);
+ if((ficrespl=fopen(filerespl,"w"))==NULL) {
+ printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+ fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+ }
+ printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+ fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+ pstamp(ficrespl);
+ fprintf(ficrespl,"# Period (stable) prevalence \n");
+ fprintf(ficrespl,"#Age ");
+ for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
+ fprintf(ficrespl,"\n");
/* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
}
printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
}
-
+ ftolpl=6*ftol*1.e5; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
/* Third parameter line */
while(fgets(line, MAXLINE, ficpar)) {
/* If line starts with a # it is a comment */
* 15 i=8 1 2 2 2
* 16 2 2 2 2
*/
- for(h=1; h <=100 ;h++){
- /* printf("h=%2d ", h); */
- /* for(k=1; k <=10; k++){ */
- /* printf("k=%d %d ",k,codtabm(h,k)); */
- /* codtab[h][k]=codtabm(h,k); */
- /* } */
- /* printf("\n"); */
- }
+ /* /\* for(h=1; h <=100 ;h++){ *\/ */
+ /* /\* printf("h=%2d ", h); *\/ */
+ /* /\* for(k=1; k <=10; k++){ *\/ */
+ /* /\* printf("k=%d %d ",k,codtabm(h,k)); *\/ */
+ /* /\* codtab[h][k]=codtabm(h,k); *\/ */
+ /* /\* } *\/ */
+ /* /\* printf("\n"); *\/ */
+ /* } */
/* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
/* 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*\/ */
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */
prlim=matrix(1,nlstate,1,nlstate);
- prevalence_limit(p, prlim, ageminpar, agemaxpar);
+ prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl);
fclose(ficrespl);
#ifdef FREEEXIT2
git://henry.ined.fr / .git / commitdiff