--- imach/src/imach.c 2017/05/23 08:39:25 1.269 +++ imach/src/imach.c 2018/05/01 17:57:25 1.287 @@ -1,6 +1,62 @@ -/* $Id: imach.c,v 1.269 2017/05/23 08:39:25 brouard Exp $ +/* $Id: imach.c,v 1.287 2018/05/01 17:57:25 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.287 2018/05/01 17:57:25 brouard + Summary: Bug fixed by providing frequencies only for non missing covariates + + Revision 1.286 2018/04/27 14:27:04 brouard + Summary: some minor bugs + + Revision 1.285 2018/04/21 21:02:16 brouard + Summary: Some bugs fixed, valgrind tested + + Revision 1.284 2018/04/20 05:22:13 brouard + Summary: Computing mean and stdeviation of fixed quantitative variables + + Revision 1.283 2018/04/19 14:49:16 brouard + Summary: Some minor bugs fixed + + Revision 1.282 2018/02/27 22:50:02 brouard + *** empty log message *** + + Revision 1.281 2018/02/27 19:25:23 brouard + Summary: Adding second argument for quitting + + Revision 1.280 2018/02/21 07:58:13 brouard + Summary: 0.99r15 + + New Makefile with recent VirtualBox 5.26. Bug in sqrt negatve in imach.c + + Revision 1.279 2017/07/20 13:35:01 brouard + Summary: temporary working + + Revision 1.278 2017/07/19 14:09:02 brouard + Summary: Bug for mobil_average=0 and prevforecast fixed(?) + + Revision 1.277 2017/07/17 08:53:49 brouard + Summary: BOM files can be read now + + Revision 1.276 2017/06/30 15:48:31 brouard + Summary: Graphs improvements + + Revision 1.275 2017/06/30 13:39:33 brouard + Summary: Saito's color + + Revision 1.274 2017/06/29 09:47:08 brouard + Summary: Version 0.99r14 + + Revision 1.273 2017/06/27 11:06:02 brouard + Summary: More documentation on projections + + Revision 1.272 2017/06/27 10:22:40 brouard + Summary: Color of backprojection changed from 6 to 5(yellow) + + Revision 1.271 2017/06/27 10:17:50 brouard + Summary: Some bug with rint + + Revision 1.270 2017/05/24 05:45:29 brouard + *** empty log message *** + Revision 1.269 2017/05/23 08:39:25 brouard Summary: Code into subroutine, cleanings @@ -1010,12 +1066,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.269 2017/05/23 08:39:25 brouard Exp $ */ +/* $Id: imach.c,v 1.287 2018/05/01 17:57:25 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; -char fullversion[]="$Revision: 1.269 $ $Date: 2017/05/23 08:39:25 $"; +char copyright[]="April 2018,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; +char fullversion[]="$Revision: 1.287 $ $Date: 2018/05/01 17:57:25 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -2489,15 +2545,18 @@ void powell(double p[], double **xi, int double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij, int nres) { - /* Computes the prevalence limit in each live state at age x and for covariate combination ij - (and selected quantitative values in nres) - by left multiplying the unit - matrix by transitions matrix until convergence is reached with precision ftolpl */ - /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ - /* Wx is row vector: population in state 1, population in state 2, population dead */ - /* or prevalence in state 1, prevalence in state 2, 0 */ - /* newm is the matrix after multiplications, its rows are identical at a factor */ - /* Initial matrix pimij */ + /**< Computes the prevalence limit in each live state at age x and for covariate combination ij + * (and selected quantitative values in nres) + * by left multiplying the unit + * matrix by transitions matrix until convergence is reached with precision ftolpl + * Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I + * Wx is row vector: population in state 1, population in state 2, population dead + * or prevalence in state 1, prevalence in state 2, 0 + * newm is the matrix after multiplications, its rows are identical at a factor. + * Inputs are the parameter, age, a tolerance for the prevalence limit ftolpl. + * Output is prlim. + * Initial matrix pimij + */ /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ /* 0, 0 , 1} */ @@ -3240,7 +3299,7 @@ double ***hbxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - agexact=age-((h-1)*hstepm + (d))*stepm/YEARM; /* age just before transition, d or d-1? */ + agexact=age-( (h-1)*hstepm + (d) )*stepm/YEARM; /* age just before transition, d or d-1? */ /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ cov[2]=agexact; if(nagesqr==1) @@ -3842,7 +3901,7 @@ void likelione(FILE *ficres,double p[], else if(mle >=1) fprintf(fichtm,"\n
File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); - + fprintf(fichtm,"\n
Equation of the model: model=1+age+%s
\n",model); for (k=1; k<= nlstate ; k++) { fprintf(fichtm,"
- Probability p%dj by origin %d and destination j. Dot's sizes are related to corresponding weight: %s-p%dj.png
\ @@ -4324,7 +4383,7 @@ void freqsummary(char fileres[], double double ***freq; /* Frequencies */ double *x, *y, a=0.,b=0.,r=1., sa=0., sb=0.; /* for regression, y=b+m*x and r is the correlation coefficient */ int no=0, linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb); - double *meanq; + double *meanq, *stdq, *idq; double **meanqt; double *pp, **prop, *posprop, *pospropt; double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; @@ -4337,6 +4396,8 @@ void freqsummary(char fileres[], double pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ + stdq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ + idq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ meanqt=matrix(1,lastpass,1,nqtveff); strcpy(fileresp,"P_"); strcat(fileresp,fileresu); @@ -4445,13 +4506,16 @@ Title=%s
Datafile=%s Firstpass=%d La posprop[i]=0; pospropt[i]=0; } - /* for (z1=1; z1<= nqfveff; z1++) { */ - /* meanq[z1]+=0.; */ + for (z1=1; z1<= nqfveff; z1++) { /* zeroing for each combination j1 as well as for the total */ + idq[z1]=0.; + meanq[z1]=0.; + stdq[z1]=0.; + } + /* for (z1=1; z1<= nqtveff; z1++) { */ /* for(m=1;m<=lastpass;m++){ */ - /* meanqt[m][z1]=0.; */ - /* } */ - /* } */ - + /* meanqt[m][z1]=0.; */ + /* } */ + /* } */ /* dateintsum=0; */ /* k2cpt=0; */ @@ -4461,9 +4525,6 @@ Title=%s
Datafile=%s Firstpass=%d La if(j !=0){ if(anyvaryingduminmodel==0){ /* If All fixed covariates */ if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - /* for (z1=1; z1<= nqfveff; z1++) { */ - /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ - /* } */ for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ /* if(Tvaraff[z1] ==-20){ */ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ @@ -4484,7 +4545,7 @@ Title=%s
Datafile=%s Firstpass=%d La }/* end j==0 */ if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ /* for(m=firstpass; m<=lastpass; m++){ */ - for(mi=1; miDatafile=%s Firstpass=%d La }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ } /* end j==0 */ /* bool =0 we keep that guy which corresponds to the combination of dummy values */ - if(bool==1){ + if(bool==1){ /*Selected */ /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] and mw[mi+1][iind]. dh depends on stepm. */ agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ @@ -4522,6 +4583,11 @@ Title=%s
Datafile=%s Firstpass=%d La if(s[m][iind]==-1) printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.)); freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ + for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */ + idq[z1]=idq[z1]+weight[iind]; + meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /* Computes mean of quantitative with selected filter */ + stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */ + } /* if((int)agev[m][iind] == 55) */ /* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ @@ -4537,6 +4603,11 @@ Title=%s
Datafile=%s Firstpass=%d La bool=1; }/* end bool 2 */ } /* end m */ + /* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */ + /* idq[z1]=idq[z1]+weight[iind]; */ + /* meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /\* Computes mean of quantitative with selected filter *\/ */ + /* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/ /\* Computes mean of quantitative with selected filter *\/ */ + /* } */ } /* end bool */ } /* end iind = 1 to imx */ /* prop[s][age] is feeded for any initial and valid live state as well as @@ -4574,6 +4645,27 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtmfr, "**********\n"); fprintf(ficlog, "**********\n"); } + /* + Printing means of quantitative variables if any + */ + for (z1=1; z1<= nqfveff; z1++) { + fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]); + fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]); + if(weightopt==1){ + printf(" Weighted mean and standard deviation of"); + fprintf(ficlog," Weighted mean and standard deviation of"); + fprintf(ficresphtmfr," Weighted mean and standard deviation of"); + } + printf(" fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + fprintf(ficlog," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + fprintf(ficresphtmfr," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)

\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + } + /* for (z1=1; z1<= nqtveff; z1++) { */ + /* for(m=1;m<=lastpass;m++){ */ + /* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f

\n", z1, m, meanqt[m][z1]); */ + /* } */ + /* } */ + fprintf(ficresphtm,""); if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ fprintf(ficresp, " Age"); @@ -4808,7 +4900,7 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficlog,"\n"); } } - } + } /* end of state i */ printf("#Freqsummary\n"); fprintf(ficlog,"\n"); for(s1=-1; s1 <=nlstate+ndeath; s1++){ @@ -4854,7 +4946,9 @@ Title=%s
Datafile=%s Firstpass=%d La fclose(ficresp); fclose(ficresphtm); fclose(ficresphtmfr); + free_vector(idq,1,nqfveff); free_vector(meanq,1,nqfveff); + free_vector(stdq,1,nqfveff); free_matrix(meanqt,1,lastpass,1,nqtveff); free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); @@ -5270,9 +5364,11 @@ void concatwav(int wav[], int **dh, int /* *cptcov=0; */ for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ + for (k=1; k <= maxncov; k++) + for(j=1; j<=2; j++) + nbcode[k][j]=0; /* Valgrind */ /* Loop on covariates without age and products and no quantitative variable */ - /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */ for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ for (j=-1; (j < maxncov); j++) Ndum[j]=0; if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ @@ -5290,7 +5386,11 @@ void concatwav(int wav[], int **dh, int modmaxcovj=ij; else if (ij < modmincovj) modmincovj=ij; - if ((ij < -1) && (ij > NCOVMAX)){ + if (ij <0 || ij >1 ){ + printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); + fprintf(ficlog,"Information, currently IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); + } + if ((ij < -1) || (ij > NCOVMAX)){ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); exit(1); }else @@ -5336,12 +5436,18 @@ void concatwav(int wav[], int **dh, int /* nbcode[Tvar[j]][3]=2; */ /* To be continued (not working yet). */ ij=0; /* ij is similar to i but can jump over null modalities */ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + + /* for (i=modmincovj; i<=modmaxcovj; i++) { */ /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + /* Skipping the case of missing values by reducing nbcode to 0 and 1 and not -1, 0, 1 */ + /* model=V1+V2+V3, if V2=-1, 0 or 1, then nbcode[2][1]=0 and nbcode[2][2]=1 instead of + * nbcode[2][1]=-1, nbcode[2][2]=0 and nbcode[2][3]=1 */ + /*, could be restored in the future */ + for (i=0; i<=1; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ if (Ndum[i] == 0) { /* If nobody responded to this modality k */ break; } ij++; - nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/ + nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1 . Could be -1*/ cptcode = ij; /* New max modality for covar j */ } /* end of loop on modality i=-1 to 1 or more */ break; @@ -5357,21 +5463,7 @@ void concatwav(int wav[], int **dh, int break; } /* end switch */ } /* end dummy test */ - - /* 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*/ + } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; /* Look at fixed dummy (single or product) covariates to check empty modalities */ @@ -5464,7 +5556,7 @@ void concatwav(int wav[], int **dh, int /* 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 agelin. - Look at hpijx to understand the reason of that which relies in memory size + Look at hpijx to understand the reason which relies in memory size consideration and note for a fixed period like estepm months */ /* 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 @@ -5745,10 +5837,11 @@ void concatwav(int wav[], int **dh, int /************ 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 *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres) { - /* Variance of health expectancies */ - /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ - /* double **newm;*/ - /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/ + /** Variance of health expectancies + * double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl); + * double **newm; + * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) + */ /* int movingaverage(); */ double **dnewm,**doldm; @@ -5756,11 +5849,11 @@ void concatwav(int wav[], int **dh, int int i, j, nhstepm, hstepm, h, nstepm ; int k; double *xp; - double **gp, **gm; /* for var eij */ - double ***gradg, ***trgradg; /*for var eij */ - double **gradgp, **trgradgp; /* for var p point j */ - double *gpp, *gmp; /* for var p point j */ - double **varppt; /* for var p point j nlstate to nlstate+ndeath */ + double **gp, **gm; /**< for var eij */ + double ***gradg, ***trgradg; /**< for var eij */ + double **gradgp, **trgradgp; /**< for var p point j */ + double *gpp, *gmp; /**< for var p point j */ + double **varppt; /**< for var p point j nlstate to nlstate+ndeath */ double ***p3mat; double age,agelim, hf; /* double ***mobaverage; */ @@ -5821,7 +5914,7 @@ void concatwav(int wav[], int **dh, int /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ fprintf(fichtm,"\n

  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); fprintf(fichtm,"\n
    %s
    \n",digitp); - /* } */ + varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); pstamp(ficresvij); fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); @@ -5876,9 +5969,12 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - + /**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and + * returns into prlim . + */ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres); - + + /* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */ if (popbased==1) { if(mobilav ==0){ for(i=1; i<=nlstate;i++) @@ -5888,23 +5984,28 @@ void concatwav(int wav[], int **dh, int prlim[i][i]=mobaverage[(int)age][i][ij]; } } - - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=1 to nhstepm */ + /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}_x\f$ at horizon h. + */ + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */ + /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}_x\f$, which are the probability + * at horizon h in state j including mortality. + */ 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]; } } - /* Next for computing probability of death (h=1 means + /* Next for computing shifted+ probability of death (h=1 means computed over hstepm matrices product = hstepm*stepm months) - as a weighted average of prlim. + as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 . */ for(j=nlstate+1;j<=nlstate+ndeath;j++){ for(i=1,gpp[j]=0.; i<= nlstate; i++) gpp[j] += prlim[i][i]*p3mat[i][j][1]; - } - /* end probability of death */ + } + + /* Again with minus shift */ for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); @@ -5937,19 +6038,23 @@ void concatwav(int wav[], int **dh, int for(i=1,gmp[j]=0.; i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; } - /* end probability of death */ - + /* end shifting computations */ + + /**< Computing gradient matrix at horizon h + */ for(j=1; j<= nlstate; j++) /* vareij */ for(h=0; h<=nhstepm; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; } - - for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */ + /**< Gradient of overall mortality p.3 (or p.j) + */ + for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; } } /* End theta */ - + + /* We got the gradient matrix for each theta and state j */ trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ for(h=0; h<=nhstepm; h++) /* veij */ @@ -5960,13 +6065,19 @@ void concatwav(int wav[], int **dh, int for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ for(theta=1; theta <=npar; theta++) trgradgp[j][theta]=gradgp[theta][j]; - + /**< as well as its transposed matrix + */ hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1;i<=nlstate;i++) for(j=1;j<=nlstate;j++) vareij[i][j][(int)age] =0.; - + + /* Computing trgradg by matcov by gradg at age and summing over h + * and k (nhstepm) formula 15 of article + * Lievre-Brouard-Heathcote + */ + for(h=0;h<=nhstepm;h++){ for(k=0;k<=nhstepm;k++){ matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); @@ -5977,7 +6088,11 @@ void concatwav(int wav[], int **dh, int } } - /* pptj */ + /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of + * p.j overall mortality formula 49 but computed directly because + * we compute the grad (wix pijx) instead of grad (pijx),even if + * wix is independent of theta. + */ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); for(j=nlstate+1;j<=nlstate+ndeath;j++) @@ -6586,7 +6701,12 @@ To be simple, these graphs help to under } /* Eigen vectors */ - v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); + if(1+(v1-lc1)*(v1-lc1)/cv12/cv12 <1.e-5){ + printf(" Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12); + fprintf(ficlog," Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12); + v11=(1./sqrt(fabs(1+(v1-lc1)*(v1-lc1)/cv12/cv12))); + }else + v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); /*v21=sqrt(1.-v11*v11); *//* error */ v21=(lc1-v1)/cv12*v11; v12=-v21; @@ -6617,8 +6737,8 @@ To be simple, these graphs help to under fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ - mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \ - mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */ + mu1,std,v11,sqrt(fabs(lc1)),v12,sqrt(fabs(lc2)), \ + mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */ }else{ first=0; fprintf(fichtmcov," %d (%.3f),",(int) age, c12); @@ -6655,8 +6775,8 @@ void printinghtml(char fileresu[], char int lastpass, int stepm, int weightopt, char model[],\ int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \ - double jprev1, double mprev1,double anprev1, double dateprev1, \ - double jprev2, double mprev2,double anprev2, double dateprev2){ + double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \ + double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){ int jj1, k1, i1, cpt, k4, nres; fprintf(fichtm,"