/* $Id$
$State$
$Log$
+ Revision 1.208 2015/11/17 14:31:57 brouard
+ Summary: temporary
+
Revision 1.207 2015/10/27 17:36:57 brouard
*** empty log message ***
/* If we start from prlim again, prlim tends to a constant matrix */
int i, ii,j,k;
- double min, max, maxmin, maxmax,sumnew=0.;
+ double *min, *max, *meandiff, maxmax,sumnew=0.;
/* double **matprod2(); */ /* test */
double **out, cov[NCOVMAX+1], **pmij();
double **newm;
- double agefin, delaymax=100. ; /* 100 Max number of years to converge */
+ double agefin, delaymax=200. ; /* 100 Max number of years to converge */
int ncvloop=0;
+ min=vector(1,nlstate);
+ max=vector(1,nlstate);
+ meandiff=vector(1,nlstate);
+
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);
savm=oldm;
oldm=newm;
- maxmax=0.;
- for(j=1;j<=nlstate;j++){
- min=1.;
- max=0.;
- for(i=1; i<=nlstate; i++) {
- sumnew=0;
- for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
+
+ for(j=1; j<=nlstate; j++){
+ max[j]=0.;
+ min[j]=1.;
+ }
+ for(i=1;i<=nlstate;i++){
+ sumnew=0;
+ for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
+ for(j=1; j<=nlstate; j++){
prlim[i][j]= newm[i][j]/(1-sumnew);
- 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); */
+ max[j]=FMAX(max[j],prlim[i][j]);
+ min[j]=FMIN(min[j],prlim[i][j]);
}
- maxmin=(max-min)/(max+min)*2;
- maxmax=FMAX(maxmax,maxmin);
- /* for(i=1; i<=nlstate; i++) { */
- /* sumnew=0.; */
- /* sumnew+=prlim[i][j]; */
- /* } */
- /* prmimj = sumnew/(float)nlstate; /\* Means of various prevalence limits. */
+ }
+
+ maxmax=0.;
+ for(j=1; j<=nlstate; j++){
+ meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */
+ maxmax=FMAX(maxmax,meandiff[j]);
+ /* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */
} /* j loop */
*ncvyear= (int)age- (int)agefin;
/* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
if(maxmax < ftolpl){
- /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
+ /* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
+ free_vector(min,1,nlstate);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
return prlim;
}
} /* age loop */
/* After some age loop it doesn't converge */
- printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g within %.0f years. \n\
+ printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
- /* printf(" age= %d newm\n",(int)age); */
- /* for(i=1; i<=nlstate+ndeath; i++) { */
- /* for(j=1;j<=nlstate+ndeath;j++){ */
- /* printf(" %lf", newm[i][j]); */
- /* } */
- /* printf("\n"); */
- /* } */
- /* printf("\n"); */
- /* printf("prlim\n"); */
- /* for(i=1; i<=nlstate; i++) { */
- /* sumnew=0; */
- /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */
- /* for(j=1;j<=nlstate;j++){ */
- /* prlim[i][j]= newm[i][j]/(1-sumnew); */
- /* printf(" %lf", prlim[i][j]); */
- /* } */
- /* printf("\n"); */
- /* } */
- /* printf("\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); */
+ free_vector(min,1,nlstate);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
+
return prlim; /* should not reach here */
}
}
/************ Variance ******************/
- void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
+ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
{
/* Variance of health expectancies */
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
/* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
nhstepm is the number of hstepm from age to agelim
nstepm is the number of stepm from age to agelim.
- Look at function hpijx to understand why (it is linked to memory size questions)
+ Look at function hpijx to understand why because of memory size limitations,
we decided (b) to get a life expectancy respecting the most precise curvature of the
survival function given by stepm (the optimization length). Unfortunately it
means that if the survival funtion is printed every two years of age and if
for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
+
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
if (popbased==1) {
if(mobilav ==0){
}
}
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */
for(j=1; j<= nlstate; j++){
for(h=0; h<=nhstepm; h++){
for(i=1, gp[h][j]=0.;i<=nlstate;i++)
gp[h][j] += prlim[i][i]*p3mat[i][j][h];
}
}
- /* This for computing probability of death (h=1 means
+ /* Next for computing probability of death (h=1 means
computed over hstepm matrices product = hstepm*stepm months)
as a weighted average of prlim.
*/
for(i=1; i<=npar; i++) /* Computes gradient x - delta */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear, ij);
+
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij);
if (popbased==1) {
if(mobilav ==0){
}
}
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
+
for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */
for(h=0; h<=nhstepm; h++){
for(i=1, gm[h][j]=0.;i<=nlstate;i++)
varppt[j][i]=doldmp[j][i];
/* end ppptj */
/* x centered again */
- hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyear,ij);
+
+ prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij);
if (popbased==1) {
if(mobilav ==0){
computed over hstepm (estepm) matrices product = hstepm*stepm months)
as a weighted average of prlim.
*/
+ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
for(j=nlstate+1;j<=nlstate+ndeath;j++){
for(i=1,gmp[j]=0.;i<= nlstate; i++)
gmp[j] += prlim[i][i]*p3mat[i][j][1];
} /* end varevsij */
/************ Variance of prevlim ******************/
- void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[])
+ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[])
{
/* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
if (stepm >= YEARM) hstepm=1;
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
gradg=matrix(1,npar,1,nlstate);
mgp=matrix(1,npar,1,nlstate);
mgm=matrix(1,npar,1,nlstate);
for(i=1; i<=npar; i++){ /* Computes gradient */
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Missing or not useful because 1 year */
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+ else
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
for(i=1;i<=nlstate;i++){
gp[i] = prlim[i][i];
mgp[theta][i] = prlim[i][i];
}
for(i=1; i<=npar; i++) /* Computes gradient */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+ else
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
for(i=1;i<=nlstate;i++){
gm[i] = prlim[i][i];
mgm[theta][i] = prlim[i][i];
}
for(i=1;i<=nlstate;i++)
gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];
+ /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */
} /* End theta */
trgradg =matrix(1,nlstate,1,npar);
for(j=1; j<=nlstate;j++)
for(theta=1; theta <=npar; theta++)
trgradg[j][theta]=gradg[theta][j];
- if((int)age==68 ||(int)age== 69 ){
- printf("\nmgm mgp %d ",(int)age);
- for(j=1; j<=nlstate;j++){
- printf("%d ",j);
- for(theta=1; theta <=npar; theta++)
- printf("%d %lf %lf",theta,mgm[theta][j],mgp[theta][j]);
- printf("\n ");
- }
- }
- if((int)age==68 ||(int)age== 69 ){
- printf("\n gradg %d ",(int)age);
- for(j=1; j<=nlstate;j++){
- printf("%d ",j);
- for(theta=1; theta <=npar; theta++)
- printf("%d %lf ",theta,gradg[theta][j]);
- printf("\n ");
- }
- }
+ /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
+ /* printf("\nmgm mgp %d ",(int)age); */
+ /* for(j=1; j<=nlstate;j++){ */
+ /* printf(" %d ",j); */
+ /* for(theta=1; theta <=npar; theta++) */
+ /* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */
+ /* printf("\n "); */
+ /* } */
+ /* } */
+ /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
+ /* printf("\n gradg %d ",(int)age); */
+ /* for(j=1; j<=nlstate;j++){ */
+ /* printf("%d ",j); */
+ /* for(theta=1; theta <=npar; theta++) */
+ /* printf("%d %lf ",theta,gradg[theta][j]); */
+ /* printf("\n "); */
+ /* } */
+ /* } */
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] =0.;
- if((int)age==68 ||(int)age== 69 ){
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81){
matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
}else{
}
- int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyear){
+ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
int i, j, k, i1 ;
/* double ftolpl = 1.e-10; */
for (age=agebase; age<=agelim; age++){
/* for (age=agebase; age<=agebase; age++){ */
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k);
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
fprintf(ficrespl,"%.0f ",age );
for(j=1;j<=cptcoveff;j++)
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
tot += prlim[i][i];
fprintf(ficrespl," %.5f", prlim[i][i]);
}
- fprintf(ficrespl," %.3f %d\n", tot, *ncvyear);
+ fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
} /* Age */
/* was end of cptcod */
} /* cptcov */
#endif
int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
int i,j, k, n=MAXN,iter=0,m,size=100, cptcod;
- int ncvyearnp=0;
- int *ncvyear=&ncvyearnp; /* Number of years needed for the period prevalence to converge */
+ int ncvyear=0; /* Number of years needed for the period prevalence to converge */
int jj, ll, li, lj, lk;
int numlinepar=0; /* Current linenumber of parameter file */
int num_filled;
if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
if (num_filled != 8) {
- printf("Not 8\n");
+ printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
+ printf("but line=%s\n",line);
}
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 *\/ */
- ftolpl=6.e-3; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
+ /*ftolpl=6.e-4; *//* 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 */
fflush(ficlog);
fflush(ficres);
-
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
+ while(fgets(line, MAXLINE, ficpar)) {
+ /* If line starts with a # it is a comment */
+ if (line[0] == '#') {
+ numlinepar++;
fputs(line,stdout);
fputs(line,ficparo);
- }
- ungetc(c,ficpar);
+ fputs(line,ficlog);
+ continue;
+ }else
+ break;
+ }
+
+ /* while((c=getc(ficpar))=='#' && c!= EOF){ */
+ /* ungetc(c,ficpar); */
+ /* fgets(line, MAXLINE, ficpar); */
+ /* fputs(line,stdout); */
+ /* fputs(line,ficparo); */
+ /* } */
+ /* ungetc(c,ficpar); */
estepm=0;
- fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
+ if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){
+
+ if (num_filled != 6) {
+ printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n");
+ printf("but line=%s\n",line);
+ goto end;
+ }
+ printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl);
+ }
+ /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
+ /*ftolpl=6.e-4;*/ /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
+
+ /* fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); */
if (estepm==0 || estepm < stepm) estepm=stepm;
if (fage <= 2) {
bage = ageminpar;
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */
prlim=matrix(1,nlstate,1,nlstate);
- prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, ncvyear);
+ prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, &ncvyear);
fclose(ficrespl);
#ifdef FREEEXIT2
cptcod= 0; /* To be deleted */
printf("varevsij %d \n",vpopbased);
fprintf(ficlog, "varevsij %d \n",vpopbased);
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
if(vpopbased==1)
fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
printf("Computing age specific period (stable) prevalences in each health state \n");
fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
for(age=bage; age <=fage ;age++){
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
if (vpopbased==1) {
if(mobilav ==0){
for(i=1; i<=nlstate;i++)
varpl=matrix(1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
- varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, strstart);
+ varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
/*}*/
}
git://henry.ined.fr / .git / commitdiff