--- imach/src/imach.c 2017/05/13 07:26:12 1.266
+++ imach/src/imach.c 2017/05/18 20:09:32 1.268
@@ -1,6 +1,12 @@
-/* $Id: imach.c,v 1.266 2017/05/13 07:26:12 brouard Exp $
+/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.268 2017/05/18 20:09:32 brouard
+ Summary: backprojection and confidence intervals of backprevalence
+
+ Revision 1.267 2017/05/13 10:25:05 brouard
+ Summary: temporary save for backprojection
+
Revision 1.266 2017/05/13 07:26:12 brouard
Summary: Version 0.99r13 (improvements and bugs fixed)
@@ -818,7 +824,7 @@ Back prevalence and projections:
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
- - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+ - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k, nres);
Computes the transition matrix starting at age 'age' over
'nhstepm*hstepm*stepm' months (i.e. until
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
@@ -986,6 +992,7 @@ typedef struct {
#define YEARM 12. /**< Number of months per year */
/* #define AGESUP 130 */
#define AGESUP 150
+#define AGEINF 0
#define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */
#define AGEBASE 40
#define AGEOVERFLOW 1.e20
@@ -1000,12 +1007,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.266 2017/05/13 07:26:12 brouard Exp $ */
+/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 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.266 $ $Date: 2017/05/13 07:26:12 $";
+char fullversion[]="$Revision: 1.268 $ $Date: 2017/05/18 20:09:32 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1079,6 +1086,8 @@ FILE *ficresvij;
char fileresv[FILENAMELENGTH];
FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];
+FILE *ficresvbl;
+char fileresvbl[FILENAMELENGTH];
char title[MAXLINE];
char model[MAXLINE]; /**< The model line */
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH];
@@ -1176,8 +1185,8 @@ double *agedc;
double **covar; /**< covar[j,i], value of jth covariate for individual i,
* covar=matrix(0,NCOVMAX,1,n);
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
-double **coqvar; /* Fixed quantitative covariate iqv */
-double ***cotvar; /* Time varying covariate itv */
+double **coqvar; /* Fixed quantitative covariate nqv */
+double ***cotvar; /* Time varying covariate ntv */
double ***cotqvar; /* Time varying quantitative covariate itqv */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
@@ -2739,7 +2748,7 @@ Earliest age to start was %d-%d=%d, ncvl
/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, ageminpar, agemaxpar, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */
/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij)); /* Bug Valgrind */
- /* if((int)age == 70){ */
+ /* if((int)age == 86 || (int)age == 87){ */
/* printf(" Backward prevalim age=%d agefin=%d \n", (int) age, (int) agefin); */
/* for(i=1; i<=nlstate+ndeath; i++) { */
/* printf("%d newm= ",i); */
@@ -2750,7 +2759,7 @@ Earliest age to start was %d-%d=%d, ncvl
/* for(j=1;j<=nlstate+ndeath;j++) { */
/* printf("%f ",oldm[i][j]); */
/* } */
- /* printf(" pmmij "); */
+ /* printf(" bmmij "); */
/* for(j=1;j<=nlstate+ndeath;j++) { */
/* printf("%f ",pmmij[i][j]); */
/* } */
@@ -2778,9 +2787,9 @@ Earliest age to start was %d-%d=%d, ncvl
meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */
maxmax=FMAX(maxmax,meandiff[i]);
/* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */
- } /* j loop */
+ } /* i loop */
*ncvyear= -( (int)age- (int)agefin);
- /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear);*/
+ /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
if(maxmax < ftolpl){
/* printf("OK Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
free_vector(min,1,nlstate);
@@ -2902,7 +2911,7 @@ double **pmij(double **ps, double *cov,
double **out, **pmij();
double sumnew=0.;
double agefin;
-
+ double k3=0.; /* constant of the w_x diagonal matrixe (in order for B to sum to 1 even for death state) */
double **dnewm, **dsavm, **doldm;
double **bbmij;
@@ -2911,55 +2920,65 @@ double **pmij(double **ps, double *cov,
dsavm=ddsavms;
agefin=cov[2];
+ /* Bx = Diag(w_x) P_x Diag(Sum_i w^i_x p^ij_x */
/* bmij *//* age is cov[2], ij is included in cov, but we need for
the observed prevalence (with this covariate ij) at beginning of transition */
/* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+
+ /* P_x */
pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */
- /* outputs pmmij which is a stochastic matrix */
- /* We do have the matrix Px in savm and we need pij */
+ /* outputs pmmij which is a stochastic matrix in row */
+
+ /* Diag(w_x) */
+ /* Problem with prevacurrent which can be zero */
+ sumnew=0.;
+ /* for (ii=1;ii<=nlstate+ndeath;ii++){ */
+ for (ii=1;ii<=nlstate;ii++){ /* Only on live states */
+ sumnew+=prevacurrent[(int)agefin][ii][ij];
+ }
+ if(sumnew >0.01){ /* At least some value in the prevalence */
+ for (ii=1;ii<=nlstate+ndeath;ii++){
+ for (j=1;j<=nlstate+ndeath;j++)
+ doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);
+ }
+ }else{
+ for (ii=1;ii<=nlstate+ndeath;ii++){
+ for (j=1;j<=nlstate+ndeath;j++)
+ doldm[ii][j]=(ii==j ? 1./nlstate : 0.0);
+ }
+ /* if(sumnew <0.9){ */
+ /* printf("Problem internal bmij B: sum on i wi <0.9: j=%d, sum_i wi=%lf,agefin=%d\n",j,sumnew, (int)agefin); */
+ /* } */
+ }
+ k3=0.0; /* We put the last diagonal to 0 */
+ for (ii=nlstate+1;ii<=nlstate+ndeath;ii++){
+ doldm[ii][ii]= k3;
+ }
+ /* End doldm, At the end doldm is diag[(w_i)] */
+
+ /* left Product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm) */
+ bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* Bug Valgrind */
+
+ /* Diag(Sum_i w^i_x p^ij_x */
+ /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */
for (j=1;j<=nlstate+ndeath;j++){
- sumnew=0.; /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */
+ sumnew=0.;
for (ii=1;ii<=nlstate;ii++){
/* sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; */
- sumnew+=pmmij[ii][j]*prevacurrent[(int)agefin][ii][ij]; /* Yes prevalence at beginning of transition */
+ sumnew+=pmmij[ii][j]*doldm[ii][ii]; /* Yes prevalence at beginning of transition */
} /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */
- if(sumnew >= 1.e-10){
- for (ii=1;ii<=nlstate+ndeath;ii++){
+ for (ii=1;ii<=nlstate+ndeath;ii++){
/* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */
- /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+ /* dsavm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
/* }else if(agefin >= agemaxpar+stepm/YEARM){ */
- /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+ /* dsavm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
/* }else */
- doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);
- } /*End ii */
- }else{ /* We put the identity matrix */
- for (ii=1;ii<=nlstate+ndeath;ii++){
- doldm[ii][j]=(ii==j ? 1. : 0.0);
- } /*End ii */
- /* printf("Problem internal bmij A: sum_i w_i*p_ij=N.j/N.. <1.e-10 i=%d, j=%d, sumnew=%lf,agefin=%d\n",ii,j,sumnew, (int)agefin); */
- }
- } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] or identity*/
- /* left Product of this diag matrix by dsavm=Px (dnewm=dsavm*doldm) */
- /* bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /\* Bug Valgrind *\/ */
- bbmij=matprod2(dnewm, pmmij,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */
- /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */
- /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
- /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
- /* left Product of this matrix by diag matrix of prevalences (savm) */
- for (j=1;j<=nlstate+ndeath;j++){
- sumnew=0.;
- for (ii=1;ii<=nlstate+ndeath;ii++){
- sumnew+=prevacurrent[(int)agefin][ii][ij];
- dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0);
- }
- /* if(sumnew <0.9){ */
- /* printf("Problem internal bmij B: sum on i wi <0.9: j=%d, sum_i wi=%lf,agefin=%d\n",j,sumnew, (int)agefin); */
- /* } */
- } /* End j, At the end dsavm is diag[(w_i)] */
- /* What if dsavm doesn't sum ii to 1? */
- /* ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /\* Bug Valgrind *\/ */
- ps=matprod2(ps, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */
- /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
+ dsavm[ii][j]=(ii==j ? 1./sumnew : 0.0);
+ } /*End ii */
+ } /* End j, At the end dsavm is diag[1/(w_1p1i+w_2 p2i)] for ALL states even if the sum is only for live states */
+
+ ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* Bug Valgrind */
+ /* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
/* end bmij */
return ps; /*pointer is unchanged */
}
@@ -3174,18 +3193,18 @@ double ***hpxij(double ***po, int nhstep
}
for(i=1; i<=nlstate+ndeath; i++)
for(j=1;j<=nlstate+ndeath;j++) {
- po[i][j][h]=newm[i][j];
- /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+ po[i][j][h]=newm[i][j];
+ /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
}
/*printf("h=%d ",h);*/
} /* end h */
- /* printf("\n H=%d \n",h); */
+ /* printf("\n H=%d \n",h); */
return po;
}
/************* Higher Back Matrix Product ***************/
/* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */
-double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij )
+double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij, int nres )
{
/* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over
'nhstepm*hstepm*stepm' months (i.e. until
@@ -3230,18 +3249,26 @@ double ***hbxij(double ***po, int nhstep
/* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
/* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
-
}
- for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
- /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
- for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
+ for (k=1; k<=nsq;k++) { /* For single varying covariates only */
+ /* Here comes the value of quantitative after renumbering k with single quantitative covariates */
+ cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
+ /* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
+ }
+ for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 */
+ if(Dummy[Tvar[Tage[k]]]){
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ } else{
+ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
+ }
+ /* printf("hBxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
+ }
+ for (k=1; k<=cptcovprod;k++){ /* Useless because included in cptcovn */
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
- /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
-
+ }
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
+
/* Careful transposed matrix */
/* age is in cov[2], prevacurrent at beginning of transition. */
/* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
@@ -3268,11 +3295,12 @@ double ***hbxij(double ***po, int nhstep
for(i=1; i<=nlstate+ndeath; i++)
for(j=1;j<=nlstate+ndeath;j++) {
po[i][j][h]=newm[i][j];
- /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+ /* if(h==nhstepm) */
+ /* printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]); */
}
- /*printf("h=%d ",h);*/
+ /* printf("h=%d %.1f ",h, agexact); */
} /* end h */
- /* printf("\n H=%d \n",h); */
+ /* printf("\n H=%d nhs=%d \n",h, nhstepm); */
return po;
}
@@ -3727,7 +3755,7 @@ double funcone( double *x)
s1=s[mw[mi][i]][i];
s2=s[mw[mi+1][i]][i];
/* if(s2==-1){ */
- /* printf(" s1=%d, s2=%d i=%d \n", s1, s2, i); */
+ /* printf(" ERROR s1=%d, s2=%d i=%d \n", s1, s2, i); */
/* /\* exit(1); *\/ */
/* } */
bbh=(double)bh[mi][i]/(double)stepm;
@@ -3760,7 +3788,7 @@ double funcone( double *x)
fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
%11.6f %11.6f %11.6f ", \
num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
- 2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
+ 2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
for(k=1,llt=0.,l=0.; k<=nlstate; k++){
llt +=ll[k]*gipmx/gsw;
fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
@@ -4280,70 +4308,7 @@ void pstamp(FILE *fichier)
fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
}
-int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) {
- /* y=a+bx regression */
- double sumx = 0.0; /* sum of x */
- double sumx2 = 0.0; /* sum of x**2 */
- double sumxy = 0.0; /* sum of x * y */
- double sumy = 0.0; /* sum of y */
- double sumy2 = 0.0; /* sum of y**2 */
- double sume2; /* sum of square or residuals */
- double yhat;
-
- double denom=0;
- int i;
- int ne=*no;
-
- for ( i=ifi, ne=0;i<=ila;i++) {
- if(!isfinite(x[i]) || !isfinite(y[i])){
- /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
- continue;
- }
- ne=ne+1;
- sumx += x[i];
- sumx2 += x[i]*x[i];
- sumxy += x[i] * y[i];
- sumy += y[i];
- sumy2 += y[i]*y[i];
- denom = (ne * sumx2 - sumx*sumx);
- /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
- }
-
- denom = (ne * sumx2 - sumx*sumx);
- if (denom == 0) {
- // vertical, slope m is infinity
- *b = INFINITY;
- *a = 0;
- if (r) *r = 0;
- return 1;
- }
-
- *b = (ne * sumxy - sumx * sumy) / denom;
- *a = (sumy * sumx2 - sumx * sumxy) / denom;
- if (r!=NULL) {
- *r = (sumxy - sumx * sumy / ne) / /* compute correlation coeff */
- sqrt((sumx2 - sumx*sumx/ne) *
- (sumy2 - sumy*sumy/ne));
- }
- *no=ne;
- for ( i=ifi, ne=0;i<=ila;i++) {
- if(!isfinite(x[i]) || !isfinite(y[i])){
- /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
- continue;
- }
- ne=ne+1;
- yhat = y[i] - *a -*b* x[i];
- sume2 += yhat * yhat ;
-
- denom = (ne * sumx2 - sumx*sumx);
- /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
- }
- *sb = sqrt(sume2/(ne-2)/(sumx2 - sumx * sumx /ne));
- *sa= *sb * sqrt(sumx2/ne);
-
- return 0;
-}
/************ Frequencies ********************/
void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
@@ -4356,8 +4321,8 @@ void freqsummary(char fileres[], double
int mi; /* Effective wave */
int first;
double ***freq; /* Frequencies */
- double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */
- int no;
+ 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 **meanqt;
double *pp, **prop, *posprop, *pospropt;
@@ -4795,6 +4760,8 @@ Title=%s
Datafile=%s Firstpass=%d La
y[iage]= log(freq[i][k][iage]/freq[i][i][iage]);
/* printf("i=%d, k=%d, s1=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,s1,j1,jj, iage, y[iage]); */
}
+ /* Some are not finite, but linreg will ignore these ages */
+ no=0;
linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */
pstart[s1]=b;
pstart[s1-1]=a;
@@ -4898,6 +4865,72 @@ Title=%s
Datafile=%s Firstpass=%d La
/* End of freqsummary */
}
+/* Simple linear regression */
+int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) {
+
+ /* y=a+bx regression */
+ double sumx = 0.0; /* sum of x */
+ double sumx2 = 0.0; /* sum of x**2 */
+ double sumxy = 0.0; /* sum of x * y */
+ double sumy = 0.0; /* sum of y */
+ double sumy2 = 0.0; /* sum of y**2 */
+ double sume2 = 0.0; /* sum of square or residuals */
+ double yhat;
+
+ double denom=0;
+ int i;
+ int ne=*no;
+
+ for ( i=ifi, ne=0;i<=ila;i++) {
+ if(!isfinite(x[i]) || !isfinite(y[i])){
+ /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+ continue;
+ }
+ ne=ne+1;
+ sumx += x[i];
+ sumx2 += x[i]*x[i];
+ sumxy += x[i] * y[i];
+ sumy += y[i];
+ sumy2 += y[i]*y[i];
+ denom = (ne * sumx2 - sumx*sumx);
+ /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+ }
+
+ denom = (ne * sumx2 - sumx*sumx);
+ if (denom == 0) {
+ // vertical, slope m is infinity
+ *b = INFINITY;
+ *a = 0;
+ if (r) *r = 0;
+ return 1;
+ }
+
+ *b = (ne * sumxy - sumx * sumy) / denom;
+ *a = (sumy * sumx2 - sumx * sumxy) / denom;
+ if (r!=NULL) {
+ *r = (sumxy - sumx * sumy / ne) / /* compute correlation coeff */
+ sqrt((sumx2 - sumx*sumx/ne) *
+ (sumy2 - sumy*sumy/ne));
+ }
+ *no=ne;
+ for ( i=ifi, ne=0;i<=ila;i++) {
+ if(!isfinite(x[i]) || !isfinite(y[i])){
+ /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+ continue;
+ }
+ ne=ne+1;
+ yhat = y[i] - *a -*b* x[i];
+ sume2 += yhat * yhat ;
+
+ denom = (ne * sumx2 - sumx*sumx);
+ /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+ }
+ *sb = sqrt(sume2/(double)(ne-2)/(sumx2 - sumx * sumx /(double)ne));
+ *sa= *sb * sqrt(sumx2/ne);
+
+ return 0;
+}
+
/************ Prevalence ********************/
void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
{
@@ -6034,7 +6067,7 @@ void concatwav(int wav[], int **dh, int
/* 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);*/
- double **dnewm,**doldm;
+ double **dnewmpar,**doldm;
int i, j, nhstepm, hstepm;
double *xp;
double *gp, *gm;
@@ -6053,7 +6086,7 @@ void concatwav(int wav[], int **dh, int
fprintf(ficresvpl,"\n");
xp=vector(1,npar);
- dnewm=matrix(1,nlstate,1,npar);
+ dnewmpar=matrix(1,nlstate,1,npar);
doldm=matrix(1,nlstate,1,nlstate);
hstepm=1*YEARM; /* Every year of age */
@@ -6123,11 +6156,11 @@ void concatwav(int wav[], int **dh, int
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] =0.;
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);
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg);
}else{
- matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
- matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg);
}
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
@@ -6148,7 +6181,132 @@ void concatwav(int wav[], int **dh, int
free_vector(xp,1,npar);
free_matrix(doldm,1,nlstate,1,npar);
- free_matrix(dnewm,1,nlstate,1,nlstate);
+ free_matrix(dnewmpar,1,nlstate,1,nlstate);
+
+}
+
+
+/************ Variance of backprevalence limit ******************/
+ void varbrevlim(char fileres[], double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres)
+{
+ /* Variance of backward 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);*/
+
+ double **dnewmpar,**doldm;
+ int i, j, nhstepm, hstepm;
+ double *xp;
+ double *gp, *gm;
+ double **gradg, **trgradg;
+ double **mgm, **mgp;
+ double age,agelim;
+ int theta;
+
+ pstamp(ficresvbl);
+ fprintf(ficresvbl,"# Standard deviation of back (stable) prevalences \n");
+ fprintf(ficresvbl,"# Age ");
+ if(nresult >=1)
+ fprintf(ficresvbl," Result# ");
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresvbl," %1d-%1d",i,i);
+ fprintf(ficresvbl,"\n");
+
+ xp=vector(1,npar);
+ dnewmpar=matrix(1,nlstate,1,npar);
+ doldm=matrix(1,nlstate,1,nlstate);
+
+ hstepm=1*YEARM; /* Every year of age */
+ hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */
+ agelim = AGEINF;
+ for (age=fage; age>=bage; age --){ /* If stepm=6 months */
+ nhstepm=(int) rint((age-agelim)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ if (stepm >= YEARM) hstepm=1;
+ nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+ gradg=matrix(1,npar,1,nlstate);
+ mgp=matrix(1,npar,1,nlstate);
+ mgm=matrix(1,npar,1,nlstate);
+ gp=vector(1,nlstate);
+ gm=vector(1,nlstate);
+
+ for(theta=1; theta <=npar; theta++){
+ for(i=1; i<=npar; i++){ /* Computes gradient */
+ xp[i] = x[i] + (i==theta ?delti[theta]:0);
+ }
+ if(mobilavproj > 0 )
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ else
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ for(i=1;i<=nlstate;i++){
+ gp[i] = bprlim[i][i];
+ mgp[theta][i] = bprlim[i][i];
+ }
+ for(i=1; i<=npar; i++) /* Computes gradient */
+ xp[i] = x[i] - (i==theta ?delti[theta]:0);
+ if(mobilavproj > 0 )
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ else
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ for(i=1;i<=nlstate;i++){
+ gm[i] = bprlim[i][i];
+ mgm[theta][i] = bprlim[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==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++)
+ varbpl[i][(int)age] =0.;
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81){
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg);
+ }else{
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg);
+ }
+ for(i=1;i<=nlstate;i++)
+ varbpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
+
+ fprintf(ficresvbl,"%.0f ",age );
+ if(nresult >=1)
+ fprintf(ficresvbl,"%d ",nres );
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresvbl," %.5f (%.5f)",bprlim[i][i],sqrt(varbpl[i][(int)age]));
+ fprintf(ficresvbl,"\n");
+ free_vector(gp,1,nlstate);
+ free_vector(gm,1,nlstate);
+ free_matrix(mgm,1,npar,1,nlstate);
+ free_matrix(mgp,1,npar,1,nlstate);
+ free_matrix(gradg,1,npar,1,nlstate);
+ free_matrix(trgradg,1,nlstate,1,npar);
+ } /* End age */
+
+ free_vector(xp,1,npar);
+ free_matrix(doldm,1,nlstate,1,npar);
+ free_matrix(dnewmpar,1,nlstate,1,nlstate);
}
@@ -6649,10 +6807,17 @@ divided by h: hPij
if(prevfcast==1){
/* Projection of prevalence up to period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\
+ fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
\
![](\"%s_%d-%d-%d.svg\")
", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
}
}
+ if(backcast==1){
+ /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg
\
+", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ }
+ }
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d-%d-%d.svg
\
@@ -6756,7 +6921,7 @@ true period expectancies (those weighted
}
/******************* Gnuplot file **************/
- void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){
char dirfileres[132],optfileres[132];
char gplotcondition[132], gplotlabel[132];
@@ -6920,6 +7085,25 @@ true period expectancies (those weighted
}
} /* end covariate */
} /* end if no covariate */
+ if(backcast == 1){
+ fprintf(ficgp,", \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+ /* k1-1 error should be nres-1*/
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+ 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 lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"\" w l lt 1");
+ } /* end if backprojcast */
} /* end if backcast */
fprintf(ficgp,"\nset out ;unset label;\n");
} /* nres */
@@ -7213,7 +7397,7 @@ set ter svg size 640, 480\nunset log y\n
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
if(m != 1 && TKresult[nres]!= k1)
continue;
- for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life ending state */
+ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */
strcpy(gplotlabel,"(");
fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt);
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
@@ -7237,11 +7421,11 @@ set ter svg size 640, 480\nunset log y\n
}
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
- fprintf(ficgp,"set label \"Ending alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
+ fprintf(ficgp,"set label \"Origin alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
k=3; /* Offset */
- for (i=1; i<= nlstate ; i ++){ /* State of origin */
+ for (i=1; i<= nlstate ; i ++){ /* State of arrival */
if(i==1)
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
else
@@ -7254,7 +7438,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (j=2; j<= nlstate ; j ++) */
/* fprintf(ficgp,"+$%d",k+l+j-1); */
/* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
- fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
+ fprintf(ficgp,") t \"bprev(%d,%d)\" w l",cpt,i);
} /* nlstate */
fprintf(ficgp,"\nset out; unset label;\n");
} /* end cpt state*/
@@ -7325,7 +7509,7 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficgp," u %d:(",ioffset);
fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \
offyear, \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate );
}else
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
@@ -7364,7 +7548,7 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficgp," u %d:(",ioffset);
fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \
offyear, \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate );
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
}else{
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
@@ -7377,6 +7561,123 @@ set ter svg size 640, 480\nunset log y\n
} /* end covariate */
} /* End if prevfcast */
+ if(backcast==1){
+ /* Back projection from cross-sectional to stable (mixed) for each covariate */
+
+ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+ for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+ if(m != 1 && TKresult[nres]!= k1)
+ continue;
+ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+ strcpy(gplotlabel,"(");
+ fprintf(ficgp,"\n#\n#\n#Back projection of prevalence to stable (mixed) back prevalence: 'BPROJ_' files, covariatecombination#=%d originstate=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
+ }
+ for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+ fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ }
+ strcpy(gplotlabel+strlen(gplotlabel),")");
+ fprintf(ficgp,"\n#\n");
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
+
+ fprintf(ficgp,"# hbijx=backprobability over h years, hb.jx is weighted by observed prev at destination state\n ");
+ fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ fprintf(ficgp,"set label \"Origin alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
+
+ /* for (i=1; i<= nlstate+1 ; i ++){ /\* nlstate +1 p11 p21 p.1 *\/ */
+ istart=nlstate+1; /* Could be one if by state, but nlstate+1 is w.i projection only */
+ /*istart=1;*/ /* Could be one if by state, but nlstate+1 is w.i projection only */
+ for (i=istart; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ if(i==istart){
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"FB_"));
+ }else{
+ fprintf(ficgp,",\\\n '' ");
+ }
+ if(cptcoveff ==0){ /* No covariate */
+ ioffset=2; /* Age is in 2 */
+ /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ /*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ fprintf(ficgp," u %d:(", ioffset);
+ if(i==nlstate+1){
+ fprintf(ficgp," $%d/(1.-$%d)):5 t 'bw%d' with line lc variable ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ fprintf(ficgp,",\\\n '' ");
+ fprintf(ficgp," u %d:(",ioffset);
+ fprintf(ficgp," (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", \
+ offbyear, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
+ }else
+ fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i );
+ }else{ /* more than 2 covariates */
+ if(cptcoveff ==1){
+ ioffset=4; /* Age is in 4 */
+ }else{
+ ioffset=6; /* Age is in 6 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ }
+ fprintf(ficgp," u %d:(",ioffset);
+ kl=0;
+ strcpy(gplotcondition,"(");
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */
+ kl++;
+ sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
+ kl++;
+ if(k 1)
+ sprintf(gplotcondition+strlen(gplotcondition)," && ");
+ }
+ strcpy(gplotcondition+strlen(gplotcondition),")");
+ /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+ /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
+ /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
+ /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
+ if(i==nlstate+1){
+ fprintf(ficgp,"%s ? $%d : 1/0):5 t 'bw%d' with line lc variable", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt );
+ fprintf(ficgp,",\\\n '' ");
+ fprintf(ficgp," u %d:(",ioffset);
+ /* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */
+ fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", gplotcondition, \
+ offbyear, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
+/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
+ }else{
+ /* fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ */
+ fprintf(ficgp,"%s ? $%d : 1/0) t 'b%d%d' with line ", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), cpt,i );
+ }
+ } /* end if covariate */
+ } /* nlstate */
+ fprintf(ficgp,"\nset out; unset label;\n");
+ } /* end cpt state*/
+ } /* end covariate */
+ } /* End if backcast */
+
/* 9eme writing MLE parameters */
fprintf(ficgp,"\n##############\n#9eme MLE estimated parameters\n#############\n");
@@ -7483,36 +7784,40 @@ set ter svg size 640, 480\nunset log y\n
/* for(j=3; j <=ncovmodel-nagesqr; j++) { */
for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
- if(j==Tage[ij]) { /* Product by age */
- if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
- if(DummyV[j]==0){
- fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
- }else{ /* quantitative */
- fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
- /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
- }
- ij++;
- }
- }else if(j==Tprod[ijp]) { /* */
- /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
- if(ijp <=cptcovprod) { /* Product */
- if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
- if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
- /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
- fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
- }else{ /* Vn is dummy and Vm is quanti */
- /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
- fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
+ if(j==Tage[ij]) { /* Product by age To be looked at!!*/
+ if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
+ if(DummyV[j]==0){
+ fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
+ }else{ /* quantitative */
+ fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
+ /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
}
- }else{ /* Vn*Vm Vn is quanti */
- if(DummyV[Tvard[ijp][2]]==0){
- fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
- }else{ /* Both quanti */
- fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ ij++;
+ }
+ }
+ }else if(cptcovprod >0){
+ if(j==Tprod[ijp]) { /* */
+ /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+ if(ijp <=cptcovprod) { /* Product */
+ if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
+ if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
+ /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
+ fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
+ }else{ /* Vn is dummy and Vm is quanti */
+ /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+ fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ }
+ }else{ /* Vn*Vm Vn is quanti */
+ if(DummyV[Tvard[ijp][2]]==0){
+ fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
+ }else{ /* Both quanti */
+ fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ }
}
+ ijp++;
}
- ijp++;
- }
+ } /* end Tprod */
} else{ /* simple covariate */
/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */
if(Dummy[j]==0){
@@ -7540,15 +7845,16 @@ set ter svg size 640, 480\nunset log y\n
ij=1;
for(j=3; j <=ncovmodel-nagesqr; j++){
- if((j-2)==Tage[ij]) { /* Bug valgrind */
- if(ij <=cptcovage) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);
- /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
- ij++;
- }
- }
- else
- fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */
+ if(cptcovage >0){
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ if(ij <=cptcovage) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
+ ij++;
+ }
+ }
+ }else
+ fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */
}
fprintf(ficgp,")");
}
@@ -7713,11 +8019,11 @@ set ter svg size 640, 480\nunset log y\n
sumr+=probs[(int)age][i][cptcod];
}
if(fabs(sum - 1.) > 1.e-3) { /* bad */
- printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, bage);
+ printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
} /* end bad */
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
if(fabs(sumr - 1.) > 1.e-3) { /* bad */
- printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, bage);
+ printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
} /* end bad */
}/* age */
@@ -7753,11 +8059,11 @@ set ter svg size 640, 480\nunset log y\n
sumr+=mobaverage[(int)age][i][cptcod];
}
if(fabs(sum - 1.) > 1.e-3) { /* bad */
- printf("Moving average B1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you decrease fage=%d?\n",cptcod, sum, (int) age, fage);
+ printf("Moving average B1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you decrease fage=%d?\n",cptcod, sum, (int) age, (int)fage);
} /* end bad */
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
if(fabs(sumr - 1.) > 1.e-3) { /* bad */
- printf("Moving average B2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase fage=%d\n",cptcod,sumr, (int)age, fage);
+ printf("Moving average B2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase fage=%d\n",cptcod,sumr, (int)age, (int)fage);
} /* end bad */
}/* age */
@@ -7806,13 +8112,13 @@ set ter svg size 640, 480\nunset log y\n
/************** Forecasting ******************/
- void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
+void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
/* proj1, year, month, day of starting projection
agemin, agemax range of age
dateprev1 dateprev2 range of dates during which prevalence is computed
anproj2 year of en of projection (same day and month as proj1).
*/
- int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */
double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double *popeffectif,*popcount;
@@ -7895,32 +8201,31 @@ set ter svg size 640, 480\nunset log y\n
nhstepm = nhstepm/hstepm;
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
+ /* We compute pii at age agec over nhstepm);*/
hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k,nres);
-
+ /* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */
for (h=0; h<=nhstepm; h++){
if (h*hstepm/YEARM*stepm ==yearp) {
- fprintf(ficresf,"\n");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
- }
- for(j=1; j<=nlstate+ndeath;j++) {
- ppij=0.;
- for(i=1; i<=nlstate;i++) {
- /* if (mobilav>=1) */
- ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
- /* else { */ /* even if mobilav==-1 we use mobaverage */
- /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
- /* } */
- if (h*hstepm/YEARM*stepm== yearp) {
- fprintf(ficresf," %.3f", p3mat[i][j][h]);
- }
- } /* end i */
- if (h*hstepm/YEARM*stepm==yearp) {
- fprintf(ficresf," %.3f", ppij);
- }
- }/* end j */
- } /* end h */
+ break;
+ }
+ }
+ fprintf(ficresf,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
+
+ for(j=1; j<=nlstate+ndeath;j++) {
+ ppij=0.;
+ for(i=1; i<=nlstate;i++) {
+ /* if (mobilav>=1) */
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
+ /* else { */ /* even if mobilav==-1 we use mobaverage */
+ /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
+ /* } */
+ fprintf(ficresf," %.3f", p3mat[i][j][h]);
+ } /* end i */
+ fprintf(ficresf," %.3f", ppij);
+ }/* end j */
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
} /* end agec */
/* diffyear=(int) anproj1+yearp-ageminpar-1; */
@@ -7935,134 +8240,151 @@ set ter svg size 640, 480\nunset log y\n
}
/* /\************** Back Forecasting ******************\/ */
-/* void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ */
-/* /\* back1, year, month, day of starting backection */
-/* agemin, agemax range of age */
-/* dateprev1 dateprev2 range of dates during which prevalence is computed */
-/* anback2 year of en of backection (same day and month as back1). */
-/* *\/ */
-/* int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; */
-/* double agec; /\* generic age *\/ */
-/* double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; */
-/* double *popeffectif,*popcount; */
-/* double ***p3mat; */
-/* /\* double ***mobaverage; *\/ */
-/* char fileresfb[FILENAMELENGTH]; */
-
-/* agelim=AGESUP; */
-/* /\* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people */
-/* in each health status at the date of interview (if between dateprev1 and dateprev2). */
-/* We still use firstpass and lastpass as another selection. */
-/* *\/ */
-/* /\* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ *\/ */
-/* /\* firstpass, lastpass, stepm, weightopt, model); *\/ */
-/* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
-
-/* strcpy(fileresfb,"FB_"); */
-/* strcat(fileresfb,fileresu); */
-/* if((ficresfb=fopen(fileresfb,"w"))==NULL) { */
-/* printf("Problem with back forecast resultfile: %s\n", fileresfb); */
-/* fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb); */
-/* } */
-/* printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
-/* fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
-
-/* if (cptcoveff==0) ncodemax[cptcoveff]=1; */
-
-/* /\* if (mobilav!=0) { *\/ */
-/* /\* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
-/* /\* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ *\/ */
-/* /\* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); *\/ */
-/* /\* printf(" Error in movingaverage mobilav=%d\n",mobilav); *\/ */
-/* /\* } *\/ */
-/* /\* } *\/ */
-
-/* stepsize=(int) (stepm+YEARM-1)/YEARM; */
-/* if (stepm<=12) stepsize=1; */
-/* if(estepm < stepm){ */
-/* printf ("Problem %d lower than %d\n",estepm, stepm); */
-/* } */
-/* else hstepm=estepm; */
-
-/* hstepm=hstepm/stepm; */
-/* yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp and */
-/* fractional in yp1 *\/ */
-/* anprojmean=yp; */
-/* yp2=modf((yp1*12),&yp); */
-/* mprojmean=yp; */
-/* yp1=modf((yp2*30.5),&yp); */
-/* jprojmean=yp; */
-/* if(jprojmean==0) jprojmean=1; */
-/* if(mprojmean==0) jprojmean=1; */
-
-/* i1=cptcoveff; */
-/* if (cptcovn < 1){i1=1;} */
+void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
+ /* back1, year, month, day of starting backection
+ agemin, agemax range of age
+ dateprev1 dateprev2 range of dates during which prevalence is computed
+ anback2 year of en of backection (same day and month as back1).
+ */
+ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ double agec; /* generic age */
+ double agelim, ppij, ppi, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double *popeffectif,*popcount;
+ double ***p3mat;
+ /* double ***mobaverage; */
+ char fileresfb[FILENAMELENGTH];
+
+ agelim=AGEINF;
+ /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+ in each health status at the date of interview (if between dateprev1 and dateprev2).
+ We still use firstpass and lastpass as another selection.
+ */
+ /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
+ /* firstpass, lastpass, stepm, weightopt, model); */
+
+ /*Do we need to compute prevalence again?*/
+
+ /* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
-/* fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); */
+ strcpy(fileresfb,"FB_");
+ strcat(fileresfb,fileresu);
+ if((ficresfb=fopen(fileresfb,"w"))==NULL) {
+ printf("Problem with back forecast resultfile: %s\n", fileresfb);
+ fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb);
+ }
+ printf("\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
+ fprintf(ficlog,"\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
-/* fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); */
-
-/* /\* if (h==(int)(YEARM*yearp)){ *\/ */
-/* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
-/* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
-/* k=k+1; */
-/* fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#"); */
-/* for(j=1;j<=cptcoveff;j++) { */
-/* fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
-/* } */
-/* fprintf(ficresfb," yearbproj age"); */
-/* for(j=1; j<=nlstate+ndeath;j++){ */
-/* for(i=1; i<=nlstate;i++) */
-/* fprintf(ficresfb," p%d%d",i,j); */
-/* fprintf(ficresfb," p.%d",j); */
-/* } */
-/* for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) { */
-/* /\* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { *\/ */
-/* fprintf(ficresfb,"\n"); */
-/* fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */
-/* for (agec=fage; agec>=(ageminpar-1); agec--){ */
-/* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
-/* nhstepm = nhstepm/hstepm; */
-/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
-/* oldm=oldms;savm=savms; */
-/* hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm,oldm,savm, dnewm, doldm, dsavm, k); */
-/* for (h=0; h<=nhstepm; h++){ */
-/* if (h*hstepm/YEARM*stepm ==yearp) { */
-/* fprintf(ficresfb,"\n"); */
-/* for(j=1;j<=cptcoveff;j++) */
-/* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
-/* fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm); */
-/* } */
-/* for(j=1; j<=nlstate+ndeath;j++) { */
-/* ppij=0.; */
-/* for(i=1; i<=nlstate;i++) { */
-/* if (mobilav==1) */
-/* ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; */
-/* else { */
-/* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; */
-/* } */
-/* if (h*hstepm/YEARM*stepm== yearp) { */
-/* fprintf(ficresfb," %.3f", p3mat[i][j][h]); */
-/* } */
-/* } /\* end i *\/ */
-/* if (h*hstepm/YEARM*stepm==yearp) { */
-/* fprintf(ficresfb," %.3f", ppij); */
-/* } */
-/* }/\* end j *\/ */
-/* } /\* end h *\/ */
-/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
-/* } /\* end agec *\/ */
-/* } /\* end yearp *\/ */
-/* } /\* end cptcod *\/ */
-/* } /\* end cptcov *\/ */
-
-/* /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
-
-/* fclose(ficresfb); */
-/* printf("End of Computing Back forecasting \n"); */
-/* fprintf(ficlog,"End of Computing Back forecasting\n"); */
+ if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ if (stepm<=12) stepsize=1;
+ if(estepm < stepm){
+ printf ("Problem %d lower than %d\n",estepm, stepm);
+ }
+ else hstepm=estepm;
+
+ hstepm=hstepm/stepm;
+ yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
+ fractional in yp1 */
+ anprojmean=yp;
+ yp2=modf((yp1*12),&yp);
+ mprojmean=yp;
+ yp1=modf((yp2*30.5),&yp);
+ jprojmean=yp;
+ if(jprojmean==0) jprojmean=1;
+ if(mprojmean==0) jprojmean=1;
+
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
+
+ fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+ printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+
+ fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
+
+ /* if (h==(int)(YEARM*yearp)){ */
+ /* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
+ /* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
+ /* k=k+1; */
+ for(nres=1; nres <= nresult; nres++) /* For each resultline */
+ for(k=1; k<=i1;k++){
+ if(i1 != 1 && TKresult[nres]!= k)
+ continue;
+ if(invalidvarcomb[k]){
+ printf("\nCombination (%d) projection ignored because no cases \n",k);
+ continue;
+ }
+ fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+ fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ }
+ fprintf(ficresfb," yearbproj age");
+ for(j=1; j<=nlstate+ndeath;j++){
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresfb," b%d%d",i,j);
+ fprintf(ficresfb," b.%d",j);
+ }
+ for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {
+ /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */
+ fprintf(ficresfb,"\n");
+ fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
+ printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
+ /* for (agec=fage; agec>=(ageminpar-1); agec--){ */
+ /* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
+ for (agec=bage; agec<=agemax-1; agec++){ /* testing */
+ /* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/
+ nhstepm=(int) rint((agec-agelim)*YEARM/stepm);
+ nhstepm = nhstepm/hstepm;
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ /* computes hbxij at age agec over 1 to nhstepm */
+ hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
+ /* hpxij(p3mat,nhstepm,agec,hstepm,p, nlstate,stepm,oldm,savm, k,nres); */
+ /* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */
+ /* printf(" agec=%.2f\n",agec);fflush(stdout); */
+ for (h=0; h<=nhstepm; h++){
+ if (h*hstepm/YEARM*stepm ==-yearp) {
+ break;
+ }
+ }
+ fprintf(ficresfb,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec-h*hstepm/YEARM*stepm);
+ for(i=1; i<=nlstate+ndeath;i++) {
+ ppij=0.;ppi=0.;
+ for(j=1; j<=nlstate;j++) {
+ /* if (mobilav==1) */
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k];
+ ppi=ppi+mobaverage[(int)agec][j][k];
+ /* else { */
+ /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
+ /* } */
+ fprintf(ficresfb," %.3f", p3mat[i][j][h]);
+ } /* end j */
+ if(ppi <0.99){
+ printf("Error in prevbackforecast, prevalence doesn't sum to 1 for state %d: %3f\n",i, ppi);
+ fprintf(ficlog,"Error in prevbackforecast, prevalence doesn't sum to 1 for state %d: %3f\n",i, ppi);
+ }
+ fprintf(ficresfb," %.3f", ppij);
+ }/* end j */
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ } /* end agec */
+ } /* end yearp */
+ } /* end k */
+
+ /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+
+ fclose(ficresfb);
+ printf("End of Computing Back forecasting \n");
+ fprintf(ficlog,"End of Computing Back forecasting\n");
-/* } */
+}
/************** Forecasting *****not tested NB*************/
/* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */
@@ -10071,7 +10393,7 @@ int hPijx(double *p, int bage, int fage)
/* oldm=oldms;savm=savms; */
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */
- hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
+ hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
for(i=1; i<=nlstate;i++)
@@ -10172,6 +10494,7 @@ int main(int argc, char *argv[])
double *delti; /* Scale */
double ***eij, ***vareij;
double **varpl; /* Variances of prevalence limits by age */
+ double **varbpl; /* Variances of back prevalence limits by age */
double *epj, vepp;
double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
@@ -10455,9 +10778,9 @@ int main(int argc, char *argv[])
covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */
- coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */
- cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/
- cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */
+ if(nqv>=1)coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */
+ if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */
+ if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
/* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
v1+v2*age+v2*v3 makes cptcovn = 3
@@ -11008,10 +11331,10 @@ Youngest age at first (selected) pass %.
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\
imx,agemin,agemax,jmin,jmax,jmean);
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
+ oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
+ newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
+ savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
+ oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
/* For Powell, parameters are in a vector p[] starting at p[1]
so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
@@ -11631,7 +11954,7 @@ Please run with mle=-1 to get a correct
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
}else{
- printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)ageminpar);
+ printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1);
}
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
@@ -11681,7 +12004,7 @@ Please run with mle=-1 to get a correct
if (mobilav!=0 ||mobilavproj !=0 ) {
mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
for(i=1;i<=AGESUP;i++)
- for(j=1;j<=nlstate;j++)
+ for(j=1;j<=nlstate+ndeath;j++)
for(k=1;k<=ncovcombmax;k++)
mobaverages[i][j][k]=0.;
mobaverage=mobaverages;
@@ -11730,15 +12053,67 @@ Please run with mle=-1 to get a correct
hBijx(p, bage, fage, mobaverage);
fclose(ficrespijb);
- free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
+ /* free_matrix(bprlim,1,nlstate,1,nlstate); /\*here or after loop ? *\/ */
+
+ prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
+ bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
+
+ /*------- Variance of back (stable) prevalence------*/
+
+ strcpy(fileresvbl,"VBL_");
+ strcat(fileresvbl,fileresu);
+ if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {
+ printf("Problem with variance of back (stable) prevalence resultfile: %s\n", fileresvbl);
+ exit(0);
+ }
+ printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);
+ fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);
+
+ /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
+
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
+
+ for(nres=1; nres <= nresult; nres++) /* For each resultline */
+ for(k=1; k<=i1;k++){
+ if(i1 != 1 && TKresult[nres]!= k)
+ continue;
+ fprintf(ficresvbl,"\n#****** ");
+ printf("\n#****** ");
+ fprintf(ficlog,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+ printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+ fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+ fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+ }
+ fprintf(ficresvbl,"******\n");
+ printf("******\n");
+ fprintf(ficlog,"******\n");
+
+ varbpl=matrix(1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+
+ varbrevlim(fileres, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, &ncvyear, k, strstart, nres);
+ free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);
+ /*}*/
+ }
+
+ fclose(ficresvbl);
+ printf("done variance-covariance of back prevalence\n");fflush(stdout);
+ fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);
- /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
- bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
}
-
+ free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
+
/* ------ Other prevalence ratios------------ */
@@ -12009,6 +12384,7 @@ Please run with mle=-1 to get a correct
fclose(ficresvpl);
printf("done variance-covariance of period prevalence\n");fflush(stdout);
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
+
free_vector(weight,1,n);
free_imatrix(Tvard,1,NCOVMAX,1,2);
@@ -12035,9 +12411,9 @@ Please run with mle=-1 to get a correct
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
- free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
- free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n);
- free_matrix(coqvar,1,maxwav,1,n);
+ if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n);
+ if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
+ if(nqv>=1)free_matrix(coqvar,1,nqv,1,n);
free_matrix(covar,0,NCOVMAX,1,n);
free_matrix(matcov,1,npar,1,npar);
free_matrix(hess,1,npar,1,npar);