/* $Id$
$State$
$Log$
+ Revision 1.191 2015/07/14 10:00:33 brouard
+ Summary: Some fixes
+
Revision 1.190 2015/05/05 08:51:13 brouard
Summary: Adding digits in output parameters (7 digits instead of 6)
/* #define DEBUG */
/* #define DEBUGBRENT */
#define POWELL /* Instead of NLOPT */
+#define POWELLF1F3 /* Skip test */
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
/* $Id$ */
/* $State$ */
-char version[]="Imach version 0.98q2, April 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char version[]="Imach version 0.98q3, July 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int m,nb;
long *num;
-int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage,*cens;
+int firstpass=0, lastpass=4,*cod, *Tage,*cens;
+int *ncodemax; /* ncodemax[j]= Number of modalities of the j th
+ covariate for which somebody answered excluding
+ undefined. Usually 2: 0 and 1. */
+int *ncodemaxwundef; /* ncodemax[j]= Number of modalities of the j th
+ covariate for which somebody answered including
+ undefined. Usually 3: -1, 0 and 1. */
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
double **pmmij, ***probs;
double *ageexmed,*agecens;
xicom[j]=xi[j];
}
- axs=0.0;
- xxss=1; /* 1 and using scale */
+ /* axs=0.0; */
+ /* xxss=1; /\* 1 and using scale *\/ */
xxs=1;
- do{
+ /* do{ */
ax=0.;
xx= xxs;
mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */
/* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */
/* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */
/* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/
- if (fx != fx){
- xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */
- printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx);
- }
- }while(fx != fx);
+ /* if (fx != fx){ */
+ /* xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */
+ /* printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); */
+ /* } */
+ /* }while(fx != fx); */
#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);
printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
- for (i=1;i<=n;i++) {
+ for (i=1;i<=n;i++) {
printf(" %d %.12f",i, p[i]);
fprintf(ficlog," %d %.12lf",i, p[i]);
fprintf(ficrespow," %.12lf", p[i]);
free_vector(ptt,1,n);
free_vector(pt,1,n);
return;
- }
+ } /* enough precision */
if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");
for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
ptt[j]=2.0*p[j]-pt[j];
pt[j]=p[j];
}
fptt=(*func)(ptt); /* f_3 */
+#ifdef POWELLF1F3
+#else
if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
+#endif
/* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */
/* From x1 (P0) distance of x2 is at h and x3 is 2h */
/* Let f"(x2) be the 2nd derivative equal everywhere. */
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("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,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
- }
+ printf("directest= %.12lf, 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,"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 */
#endif
#ifdef DEBUGLINMIN
printf("\n");
fprintf(ficlog,"\n");
#endif
- } /* end of t negative */
+ } /* end of t or directest negative */
+#ifdef POWELLF1F3
+#else
} /* end if (fptt < fp) */
- }
+#endif
+ } /* loop iteration */
}
/**** Prevalence limit (stable or period prevalence) ****************/
cptcoveff=0;
- for (k=-1; k < maxncov; k++) Ndum[k]=0;
for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
/* Loop on covariates without age and products */
for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */
+ for (k=-1; k < maxncov; k++) Ndum[k]=0;
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the
modality of this covariate Vj*/
ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
/* getting the maximum value of the modality of the covariate
(should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
female is 1, then modmaxcovj=1.*/
- } /* end for loop on individuals */
+ } /* end for loop on individuals i */
printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
+ fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
cptcode=modmaxcovj;
/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
/*for (i=0; i<=cptcode; i++) {*/
- for (i=modmincovj; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */
- printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], i, Ndum[i]);
- if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */
- ncodemax[j]++; /* ncodemax[j]= Number of non-null modalities of the j th covariate. */
+ for (k=modmincovj; k<=modmaxcovj; k++) { /* k=-1 ? 0 and 1*//* For each value k of the modality of model-cov j */
+ printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
+ fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
+ if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */
+ if( k != -1){
+ ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th
+ covariate for which somebody answered excluding
+ undefined. Usually 2: 0 and 1. */
+ }
+ ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
+ covariate for which somebody answered including
+ undefined. Usually 3: -1, 0 and 1. */
}
/* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for
historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
nbcode[Tvar[j]][2]=1;
nbcode[Tvar[j]][3]=2;
*/
- ij=1; /* ij is similar to i but can jumps over null modalities */
- for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 */
- for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */
- /*recode from 0 */
- if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */
- nbcode[Tvar[j]][ij]=k; /* stores the modality k in an array nbcode.
- k is a modality. If we have model=V1+V1*sex
- then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
- ij++;
+ ij=0; /* ij is similar to i but can jumps over null modalities */
+ for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 to 1*/
+ if (Ndum[i] == 0) { /* If at least one individual responded to this modality k */
+ break;
}
- if (ij > ncodemax[j]) break;
- } /* end of loop on */
- } /* end of loop on modality */
+ ij++;
+ nbcode[Tvar[j]][ij]=i; /* stores the original modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
+ cptcode = ij; /* New max modality for covar j */
+ } /* end of loop on modality i=-1 to 1 or more */
+
+ /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
+ /* /\*recode from 0 *\/ */
+ /* k is a modality. If we have model=V1+V1*sex */
+ /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
+ /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */
+ /* } */
+ /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */
+ /* if (ij > ncodemax[j]) { */
+ /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
+ /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
+ /* break; */
+ /* } */
+ /* } /\* end of loop on modality k *\/ */
} /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/
for (k=-1; k< maxncov; k++) Ndum[k]=0;
Ndum[ij]++; /* Might be supersed V1 + V1*age */
}
- ij=1;
+ ij=0;
for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
if((Ndum[i]!=0) && (i<=ncovcol)){
+ ij++;
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
Tvaraff[ij]=i; /*For printing (unclear) */
- ij++;
- }else
- Tvaraff[ij]=0;
+ }else{
+ /* Tvaraff[ij]=0; */
+ }
}
- ij--;
+ /* ij--; */
cptcoveff=ij; /*Number of total covariates*/
}
jj1=0;
for(k1=1; k1<=m;k1++){
- for(i1=1; i1<=ncodemax[k1];i1++){
+ /* for(i1=1; i1<=ncodemax[k1];i1++){ */
jj1++;
if (cptcovn > 0) {
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
- for (cpt=1; cpt<=cptcoveff;cpt++)
+ for (cpt=1; cpt<=cptcoveff;cpt++){
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
+ printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);fflush(stdout);
+ }
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}
/* Pij */
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
<img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
}
- } /* end i1 */
+ /* } /\* end i1 *\/ */
}/* End k1 */
fprintf(fichtm,"</ul>");
jj1=0;
for(k1=1; k1<=m;k1++){
- for(i1=1; i1<=ncodemax[k1];i1++){
+ /* for(i1=1; i1<=ncodemax[k1];i1++){ */
jj1++;
if (cptcovn > 0) {
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
drawn in addition to the population based expectancies computed using\
observed and cahotic prevalences: %s%d.png<br>\
<img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
- } /* end i1 */
+ /* } /\* end i1 *\/ */
}/* End k1 */
fprintf(fichtm,"</ul>");
fflush(fichtm);
if (strlen(model) >1){ /* If there is at least 1 covariate */
j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0;
if (strstr(model,"AGE") !=0){
- printf("Error. AGE must be in lower case 'age' model=1+age+%s ",model);
- fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s ",model);fflush(ficlog);
+ printf("Error. AGE must be in lower case 'age' model=1+age+%s. ",model);
+ fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s. ",model);fflush(ficlog);
return 1;
}
if (strstr(model,"v") !=0){
printf(" strpt=%s, model=%s\n",strpt, model);
if(strpt != model){
printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
- 'model=1+age+age*age+V1' or 'model=1+age+age*age+V1+V1*age', please swap as well as \n \
+ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model);
fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
- 'model=1+age+age*age+V1' or 'model=1+age+age*age+V1+V1*age', please swap as well as \n \
+ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model); fflush(ficlog);
return 1;
}
}
else if(mle==-3) { /* Main Wizard */
prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
- printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
- fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
+ printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
+ fprintf(ficlog," You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
matcov=matrix(1,npar,1,npar);
}
s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */
tab=ivector(1,NCOVMAX);
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
+ ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
/* Reads data from file datafile */
if (readdata(datafile, firstobs, lastobs, &imx)==1)
}
/*--------- results files --------------*/
- fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
+ fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
free_ivector(ncodemax,1,NCOVMAX);
+ free_ivector(ncodemaxwundef,1,NCOVMAX);
free_ivector(Tvar,1,NCOVMAX);
free_ivector(Tprod,1,NCOVMAX);
free_ivector(Tvaraff,1,NCOVMAX);
git://henry.ined.fr / .git / commitdiff