|
version 1.267, 2017/05/13 10:25:05
|
version 1.282, 2018/02/27 22:50:02
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
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 |
| |
|
| |
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 |
Revision 1.267 2017/05/13 10:25:05 brouard |
| Summary: temporary save for backprojection |
Summary: temporary save for backprojection |
| |
|
|
Line 989 typedef struct {
|
Line 1036 typedef struct {
|
| #define YEARM 12. /**< Number of months per year */ |
#define YEARM 12. /**< Number of months per year */ |
| /* #define AGESUP 130 */ |
/* #define AGESUP 130 */ |
| #define AGESUP 150 |
#define AGESUP 150 |
| |
#define AGEINF 0 |
| #define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */ |
#define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */ |
| #define AGEBASE 40 |
#define AGEBASE 40 |
| #define AGEOVERFLOW 1.e20 |
#define AGEOVERFLOW 1.e20 |
|
Line 1080 FILE *ficrescveij;
|
Line 1128 FILE *ficrescveij;
|
| char filerescve[FILENAMELENGTH]; |
char filerescve[FILENAMELENGTH]; |
| FILE *ficresvij; |
FILE *ficresvij; |
| char fileresv[FILENAMELENGTH]; |
char fileresv[FILENAMELENGTH]; |
| FILE *ficresvpl; |
|
| char fileresvpl[FILENAMELENGTH]; |
|
| char title[MAXLINE]; |
char title[MAXLINE]; |
| char model[MAXLINE]; /**< The model line */ |
char model[MAXLINE]; /**< The model line */ |
| char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH]; |
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH]; |
|
Line 1179 double *agedc;
|
Line 1226 double *agedc;
|
| double **covar; /**< covar[j,i], value of jth covariate for individual i, |
double **covar; /**< covar[j,i], value of jth covariate for individual i, |
| * covar=matrix(0,NCOVMAX,1,n); |
* covar=matrix(0,NCOVMAX,1,n); |
| * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ |
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ |
| double **coqvar; /* Fixed quantitative covariate iqv */ |
double **coqvar; /* Fixed quantitative covariate nqv */ |
| double ***cotvar; /* Time varying covariate itv */ |
double ***cotvar; /* Time varying covariate ntv */ |
| double ***cotqvar; /* Time varying quantitative covariate itqv */ |
double ***cotqvar; /* Time varying quantitative covariate itqv */ |
| double idx; |
double idx; |
| int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
|
Line 2483 void powell(double p[], double **xi, int
|
Line 2530 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) |
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 |
/**< Computes the prevalence limit in each live state at age x and for covariate combination ij |
| (and selected quantitative values in nres) |
* (and selected quantitative values in nres) |
| by left multiplying the unit |
* by left multiplying the unit |
| matrix by transitions matrix until convergence is reached with precision ftolpl */ |
* 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= 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 */ |
* Wx is row vector: population in state 1, population in state 2, population dead |
| /* or prevalence in state 1, prevalence in state 2, 0 */ |
* or prevalence in state 1, prevalence in state 2, 0 |
| /* newm is the matrix after multiplications, its rows are identical at a factor */ |
* newm is the matrix after multiplications, its rows are identical at a factor. |
| /* Initial matrix pimij */ |
* 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.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ |
| /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ |
/* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ |
| /* 0, 0 , 1} */ |
/* 0, 0 , 1} */ |
|
Line 2742 Earliest age to start was %d-%d=%d, ncvl
|
Line 2792 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, 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, 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 */ |
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); */ |
/* printf(" Backward prevalim age=%d agefin=%d \n", (int) age, (int) agefin); */ |
| /* for(i=1; i<=nlstate+ndeath; i++) { */ |
/* for(i=1; i<=nlstate+ndeath; i++) { */ |
| /* printf("%d newm= ",i); */ |
/* printf("%d newm= ",i); */ |
|
Line 2753 Earliest age to start was %d-%d=%d, ncvl
|
Line 2803 Earliest age to start was %d-%d=%d, ncvl
|
| /* for(j=1;j<=nlstate+ndeath;j++) { */ |
/* for(j=1;j<=nlstate+ndeath;j++) { */ |
| /* printf("%f ",oldm[i][j]); */ |
/* printf("%f ",oldm[i][j]); */ |
| /* } */ |
/* } */ |
| /* printf(" pmmij "); */ |
/* printf(" bmmij "); */ |
| /* for(j=1;j<=nlstate+ndeath;j++) { */ |
/* for(j=1;j<=nlstate+ndeath;j++) { */ |
| /* printf("%f ",pmmij[i][j]); */ |
/* printf("%f ",pmmij[i][j]); */ |
| /* } */ |
/* } */ |
|
Line 2781 Earliest age to start was %d-%d=%d, ncvl
|
Line 2831 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 */ |
meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */ |
| maxmax=FMAX(maxmax,meandiff[i]); |
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); */ |
/* 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); |
*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){ |
if(maxmax < ftolpl){ |
| /* printf("OK Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
/* 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); |
free_vector(min,1,nlstate); |
|
Line 2905 double **pmij(double **ps, double *cov,
|
Line 2955 double **pmij(double **ps, double *cov,
|
| double **out, **pmij(); |
double **out, **pmij(); |
| double sumnew=0.; |
double sumnew=0.; |
| double agefin; |
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 **dnewm, **dsavm, **doldm; |
| double **bbmij; |
double **bbmij; |
| |
|
|
Line 2914 double **pmij(double **ps, double *cov,
|
Line 2964 double **pmij(double **ps, double *cov,
|
| dsavm=ddsavms; |
dsavm=ddsavms; |
| |
|
| agefin=cov[2]; |
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 |
/* 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 */ |
the observed prevalence (with this covariate ij) at beginning of transition */ |
| /* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
/* 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 */ |
pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */ |
| /* outputs pmmij which is a stochastic matrix */ |
/* outputs pmmij which is a stochastic matrix in row */ |
| /* We do have the matrix Px in savm and we need pij */ |
|
| |
/* 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 */ |
| |
/* printf(" agefin=%d, ii=%d, ij=%d, prev=%f\n",(int)agefin,ii, ij, prevacurrent[(int)agefin][ii][ij]); */ |
| |
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++){ |
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++){ |
for (ii=1;ii<=nlstate;ii++){ |
| /* sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; */ |
/* 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 */ |
} /* 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){ */ |
/* 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){ */ |
/* }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 */ |
/* }else */ |
| doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
dsavm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
| } /*End ii */ |
} /*End ii */ |
| }else{ /* We put the identity matrix */ |
} /* 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 */ |
| for (ii=1;ii<=nlstate+ndeath;ii++){ |
|
| doldm[ii][j]=(ii==j ? 1. : 0.0); |
ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* Bug Valgrind */ |
| } /*End ii */ |
/* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ |
| /* 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)] */ |
|
| /* end bmij */ |
/* end bmij */ |
| return ps; /*pointer is unchanged */ |
return ps; /*pointer is unchanged */ |
| } |
} |
|
Line 3223 double ***hbxij(double ***po, int nhstep
|
Line 3284 double ***hbxij(double ***po, int nhstep
|
| newm=savm; |
newm=savm; |
| /* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
| cov[1]=1.; |
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 *\/ */ |
/* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ |
| cov[2]=agexact; |
cov[2]=agexact; |
| if(nagesqr==1) |
if(nagesqr==1) |
|
Line 3279 double ***hbxij(double ***po, int nhstep
|
Line 3340 double ***hbxij(double ***po, int nhstep
|
| for(i=1; i<=nlstate+ndeath; i++) |
for(i=1; i<=nlstate+ndeath; i++) |
| for(j=1;j<=nlstate+ndeath;j++) { |
for(j=1;j<=nlstate+ndeath;j++) { |
| po[i][j][h]=newm[i][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 */ |
} /* end h */ |
| /* printf("\n H=%d \n",h); */ |
/* printf("\n H=%d nhs=%d \n",h, nhstepm); */ |
| return po; |
return po; |
| } |
} |
| |
|
|
Line 3738 double funcone( double *x)
|
Line 3800 double funcone( double *x)
|
| s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
| s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
| /* if(s2==-1){ */ |
/* 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); *\/ */ |
/* /\* exit(1); *\/ */ |
| /* } */ |
/* } */ |
| bbh=(double)bh[mi][i]/(double)stepm; |
bbh=(double)bh[mi][i]/(double)stepm; |
|
Line 3771 double funcone( double *x)
|
Line 3833 double funcone( double *x)
|
| fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ |
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 ", \ |
%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, |
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++){ |
for(k=1,llt=0.,l=0.; k<=nlstate; k++){ |
| llt +=ll[k]*gipmx/gsw; |
llt +=ll[k]*gipmx/gsw; |
| fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); |
fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); |
|
Line 3824 void likelione(FILE *ficres,double p[],
|
Line 3886 void likelione(FILE *ficres,double p[],
|
| else if(mle >=1) |
else if(mle >=1) |
| fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); |
fprintf(fichtm,"\n<br>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: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk)); |
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: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk)); |
| |
fprintf(fichtm,"\n<br>Equation of the model: <b>model=1+age+%s</b><br>\n",model); |
| |
|
| for (k=1; k<= nlstate ; k++) { |
for (k=1; k<= nlstate ; k++) { |
| fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ |
fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ |
|
Line 4291 void pstamp(FILE *fichier)
|
Line 4353 void pstamp(FILE *fichier)
|
| fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); |
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 ********************/ |
/************ Frequencies ********************/ |
| void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ |
void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ |
|
Line 4367 void freqsummary(char fileres[], double
|
Line 4366 void freqsummary(char fileres[], double
|
| int mi; /* Effective wave */ |
int mi; /* Effective wave */ |
| int first; |
int first; |
| double ***freq; /* Frequencies */ |
double ***freq; /* Frequencies */ |
| double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */ |
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; |
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; |
| double **meanqt; |
double **meanqt; |
| double *pp, **prop, *posprop, *pospropt; |
double *pp, **prop, *posprop, *pospropt; |
|
Line 4806 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4805 Title=%s <br>Datafile=%s Firstpass=%d La
|
| y[iage]= log(freq[i][k][iage]/freq[i][i][iage]); |
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]); */ |
/* 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 */ |
linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */ |
| pstart[s1]=b; |
pstart[s1]=b; |
| pstart[s1-1]=a; |
pstart[s1-1]=a; |
|
Line 4909 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4910 Title=%s <br>Datafile=%s Firstpass=%d La
|
| /* End of freqsummary */ |
/* 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 ********************/ |
/************ 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) |
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) |
| { |
{ |
|
Line 5441 void concatwav(int wav[], int **dh, int
|
Line 5508 void concatwav(int wav[], int **dh, int
|
| /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. |
/* 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 |
nhstepm is the number of hstepm from age to agelim |
| nstepm is the number of stepm from age to agelin. |
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 */ |
and note for a fixed period like estepm months */ |
| /* We decided (b) to get a life expectancy respecting the most precise curvature of the |
/* 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 |
survival function given by stepm (the optimization length). Unfortunately it |
|
Line 5672 void concatwav(int wav[], int **dh, int
|
Line 5739 void concatwav(int wav[], int **dh, int
|
| /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ |
/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ |
| |
|
| } |
} |
| |
|
| |
/* Standard deviation of expectancies ij */ |
| fprintf(ficresstdeij,"%3.0f",age ); |
fprintf(ficresstdeij,"%3.0f",age ); |
| for(i=1; i<=nlstate;i++){ |
for(i=1; i<=nlstate;i++){ |
| eip=0.; |
eip=0.; |
|
Line 5687 void concatwav(int wav[], int **dh, int
|
Line 5755 void concatwav(int wav[], int **dh, int
|
| } |
} |
| fprintf(ficresstdeij,"\n"); |
fprintf(ficresstdeij,"\n"); |
| |
|
| |
/* Variance of expectancies ij */ |
| fprintf(ficrescveij,"%3.0f",age ); |
fprintf(ficrescveij,"%3.0f",age ); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| for(j=1; j<=nlstate;j++){ |
for(j=1; j<=nlstate;j++){ |
|
Line 5720 void concatwav(int wav[], int **dh, int
|
Line 5789 void concatwav(int wav[], int **dh, int
|
| /************ Variance ******************/ |
/************ 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) |
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 */ |
/** Variance of health expectancies |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ |
* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl); |
| /* double **newm;*/ |
* double **newm; |
| /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/ |
* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) |
| |
*/ |
| |
|
| /* int movingaverage(); */ |
/* int movingaverage(); */ |
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
|
Line 5731 void concatwav(int wav[], int **dh, int
|
Line 5801 void concatwav(int wav[], int **dh, int
|
| int i, j, nhstepm, hstepm, h, nstepm ; |
int i, j, nhstepm, hstepm, h, nstepm ; |
| int k; |
int k; |
| double *xp; |
double *xp; |
| double **gp, **gm; /* for var eij */ |
double **gp, **gm; /**< for var eij */ |
| double ***gradg, ***trgradg; /*for var eij */ |
double ***gradg, ***trgradg; /**< for var eij */ |
| double **gradgp, **trgradgp; /* for var p point j */ |
double **gradgp, **trgradgp; /**< for var p point j */ |
| double *gpp, *gmp; /* 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 **varppt; /**< for var p point j nlstate to nlstate+ndeath */ |
| double ***p3mat; |
double ***p3mat; |
| double age,agelim, hf; |
double age,agelim, hf; |
| /* double ***mobaverage; */ |
/* double ***mobaverage; */ |
|
Line 5796 void concatwav(int wav[], int **dh, int
|
Line 5866 void concatwav(int wav[], int **dh, int
|
| /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ |
/* fprintf(fichtm, "#Local time at start: %s", strstart);*/ |
| fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); |
fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); |
| fprintf(fichtm,"\n<br>%s <br>\n",digitp); |
fprintf(fichtm,"\n<br>%s <br>\n",digitp); |
| /* } */ |
|
| varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); |
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); |
| pstamp(ficresvij); |
pstamp(ficresvij); |
| fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); |
fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); |
|
Line 5851 void concatwav(int wav[], int **dh, int
|
Line 5921 void concatwav(int wav[], int **dh, int
|
| for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ |
for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ |
| xp[i] = x[i] + (i==theta ?delti[theta]:0); |
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); |
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 (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
|
Line 5863 void concatwav(int wav[], int **dh, int
|
Line 5936 void concatwav(int wav[], int **dh, int
|
| prlim[i][i]=mobaverage[(int)age][i][ij]; |
prlim[i][i]=mobaverage[(int)age][i][ij]; |
| } |
} |
| } |
} |
| |
/**< 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=1 to nhstepm */ |
*/ |
| |
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(j=1; j<= nlstate; j++){ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
| gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
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) |
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(j=nlstate+1;j<=nlstate+ndeath;j++){ |
| for(i=1,gpp[j]=0.; i<= nlstate; i++) |
for(i=1,gpp[j]=0.; i<= nlstate; i++) |
| gpp[j] += prlim[i][i]*p3mat[i][j][1]; |
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 */ |
for(i=1; i<=npar; i++) /* Computes gradient x - delta */ |
| xp[i] = x[i] - (i==theta ?delti[theta]:0); |
xp[i] = x[i] - (i==theta ?delti[theta]:0); |
|
Line 5912 void concatwav(int wav[], int **dh, int
|
Line 5990 void concatwav(int wav[], int **dh, int
|
| for(i=1,gmp[j]=0.; i<= nlstate; i++) |
for(i=1,gmp[j]=0.; i<= nlstate; i++) |
| gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
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(j=1; j<= nlstate; j++) /* vareij */ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; |
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; |
| } |
} |
| |
/**< Gradient of overall mortality p.3 (or p.j) |
| for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */ |
*/ |
| |
for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */ |
| gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; |
gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; |
| } |
} |
| |
|
| } /* End theta */ |
} /* End theta */ |
| |
|
| |
/* We got the gradient matrix for each theta and state j */ |
| trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ |
trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ |
| |
|
| for(h=0; h<=nhstepm; h++) /* veij */ |
for(h=0; h<=nhstepm; h++) /* veij */ |
|
Line 5935 void concatwav(int wav[], int **dh, int
|
Line 6017 void concatwav(int wav[], int **dh, int
|
| for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ |
for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ |
| for(theta=1; theta <=npar; theta++) |
for(theta=1; theta <=npar; theta++) |
| trgradgp[j][theta]=gradgp[theta][j]; |
trgradgp[j][theta]=gradgp[theta][j]; |
| |
/**< as well as its transposed matrix |
| |
*/ |
| |
|
| hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ |
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ |
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| for(j=1;j<=nlstate;j++) |
for(j=1;j<=nlstate;j++) |
| vareij[i][j][(int)age] =0.; |
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(h=0;h<=nhstepm;h++){ |
| for(k=0;k<=nhstepm;k++){ |
for(k=0;k<=nhstepm;k++){ |
| matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); |
matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); |
|
Line 5952 void concatwav(int wav[], int **dh, int
|
Line 6040 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(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); |
matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); |
| for(j=nlstate+1;j<=nlstate+ndeath;j++) |
for(j=nlstate+1;j<=nlstate+ndeath;j++) |
|
Line 6040 void concatwav(int wav[], int **dh, int
|
Line 6132 void concatwav(int wav[], int **dh, int
|
| } /* end varevsij */ |
} /* end varevsij */ |
| |
|
| /************ Variance of prevlim ******************/ |
/************ Variance of prevlim ******************/ |
| void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres) |
void varprevlim(char fileresvpl[], FILE *ficresvpl, double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres) |
| { |
{ |
| /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ |
/* 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 **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; |
int i, j, nhstepm, hstepm; |
| double *xp; |
double *xp; |
| double *gp, *gm; |
double *gp, *gm; |
|
Line 6064 void concatwav(int wav[], int **dh, int
|
Line 6156 void concatwav(int wav[], int **dh, int
|
| fprintf(ficresvpl,"\n"); |
fprintf(ficresvpl,"\n"); |
| |
|
| xp=vector(1,npar); |
xp=vector(1,npar); |
| dnewm=matrix(1,nlstate,1,npar); |
dnewmpar=matrix(1,nlstate,1,npar); |
| doldm=matrix(1,nlstate,1,nlstate); |
doldm=matrix(1,nlstate,1,nlstate); |
| |
|
| hstepm=1*YEARM; /* Every year of age */ |
hstepm=1*YEARM; /* Every year of age */ |
|
Line 6134 void concatwav(int wav[], int **dh, int
|
Line 6226 void concatwav(int wav[], int **dh, int
|
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| varpl[i][(int)age] =0.; |
varpl[i][(int)age] =0.; |
| if((int)age==79 ||(int)age== 80 ||(int)age== 81){ |
if((int)age==79 ||(int)age== 80 ||(int)age== 81){ |
| matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); |
matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov); |
| matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); |
matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg); |
| }else{ |
}else{ |
| matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); |
matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov); |
| matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); |
matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg); |
| } |
} |
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ |
varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ |
|
Line 6159 void concatwav(int wav[], int **dh, int
|
Line 6251 void concatwav(int wav[], int **dh, int
|
| |
|
| free_vector(xp,1,npar); |
free_vector(xp,1,npar); |
| free_matrix(doldm,1,nlstate,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 one-step probabilities ******************/ |
|
| void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) |
|
| { |
|
| int i, j=0, k1, l1, tj; |
|
| int k2, l2, j1, z1; |
|
| int k=0, l; |
|
| int first=1, first1, first2; |
|
| double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
|
| double **dnewm,**doldm; |
|
| double *xp; |
|
| double *gp, *gm; |
|
| double **gradg, **trgradg; |
|
| double **mu; |
|
| double age, cov[NCOVMAX+1]; |
|
| double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ |
|
| int theta; |
|
| char fileresprob[FILENAMELENGTH]; |
|
| char fileresprobcov[FILENAMELENGTH]; |
|
| char fileresprobcor[FILENAMELENGTH]; |
|
| double ***varpij; |
|
| |
|
| strcpy(fileresprob,"PROB_"); |
/************ Variance of backprevalence limit ******************/ |
| strcat(fileresprob,fileres); |
void varbrevlim(char fileresvbl[], FILE *ficresvbl, 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) |
| if((ficresprob=fopen(fileresprob,"w"))==NULL) { |
{ |
| printf("Problem with resultfile: %s\n", fileresprob); |
/* Variance of backward prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
| } |
|
| strcpy(fileresprobcov,"PROBCOV_"); |
double **dnewmpar,**doldm; |
| strcat(fileresprobcov,fileresu); |
int i, j, nhstepm, hstepm; |
| if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { |
double *xp; |
| printf("Problem with resultfile: %s\n", fileresprobcov); |
double *gp, *gm; |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); |
double **gradg, **trgradg; |
| } |
double **mgm, **mgp; |
| strcpy(fileresprobcor,"PROBCOR_"); |
double age,agelim; |
| strcat(fileresprobcor,fileresu); |
int theta; |
| if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { |
|
| printf("Problem with resultfile: %s\n", fileresprobcor); |
pstamp(ficresvbl); |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); |
fprintf(ficresvbl,"# Standard deviation of back (stable) prevalences \n"); |
| } |
fprintf(ficresvbl,"# Age "); |
| printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
if(nresult >=1) |
| fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
fprintf(ficresvbl," Result# "); |
| printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
for(i=1; i<=nlstate;i++) |
| fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
fprintf(ficresvbl," %1d-%1d",i,i); |
| printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
fprintf(ficresvbl,"\n"); |
| fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
|
| pstamp(ficresprob); |
xp=vector(1,npar); |
| fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); |
dnewmpar=matrix(1,nlstate,1,npar); |
| fprintf(ficresprob,"# Age"); |
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); |
| |
|
| |
} |
| |
|
| |
/************ Variance of one-step probabilities ******************/ |
| |
void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) |
| |
{ |
| |
int i, j=0, k1, l1, tj; |
| |
int k2, l2, j1, z1; |
| |
int k=0, l; |
| |
int first=1, first1, first2; |
| |
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
| |
double **dnewm,**doldm; |
| |
double *xp; |
| |
double *gp, *gm; |
| |
double **gradg, **trgradg; |
| |
double **mu; |
| |
double age, cov[NCOVMAX+1]; |
| |
double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ |
| |
int theta; |
| |
char fileresprob[FILENAMELENGTH]; |
| |
char fileresprobcov[FILENAMELENGTH]; |
| |
char fileresprobcor[FILENAMELENGTH]; |
| |
double ***varpij; |
| |
|
| |
strcpy(fileresprob,"PROB_"); |
| |
strcat(fileresprob,fileres); |
| |
if((ficresprob=fopen(fileresprob,"w"))==NULL) { |
| |
printf("Problem with resultfile: %s\n", fileresprob); |
| |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); |
| |
} |
| |
strcpy(fileresprobcov,"PROBCOV_"); |
| |
strcat(fileresprobcov,fileresu); |
| |
if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { |
| |
printf("Problem with resultfile: %s\n", fileresprobcov); |
| |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); |
| |
} |
| |
strcpy(fileresprobcor,"PROBCOR_"); |
| |
strcat(fileresprobcor,fileresu); |
| |
if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { |
| |
printf("Problem with resultfile: %s\n", fileresprobcor); |
| |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); |
| |
} |
| |
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
| |
fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
| |
printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
| |
fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
| |
printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
| |
fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
| |
pstamp(ficresprob); |
| |
fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); |
| |
fprintf(ficresprob,"# Age"); |
| pstamp(ficresprobcov); |
pstamp(ficresprobcov); |
| fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); |
fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); |
| fprintf(ficresprobcov,"# Age"); |
fprintf(ficresprobcov,"# Age"); |
|
Line 6436 To be simple, these graphs help to under
|
Line 6653 To be simple, these graphs help to under
|
| } |
} |
| |
|
| /* Eigen vectors */ |
/* 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=sqrt(1.-v11*v11); *//* error */ |
| v21=(lc1-v1)/cv12*v11; |
v21=(lc1-v1)/cv12*v11; |
| v12=-v21; |
v12=-v21; |
|
Line 6467 To be simple, these graphs help to under
|
Line 6689 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,"\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,"\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", \ |
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)), \ |
mu1,std,v11,sqrt(fabs(lc1)),v12,sqrt(fabs(lc2)), \ |
| mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */ |
mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */ |
| }else{ |
}else{ |
| first=0; |
first=0; |
| fprintf(fichtmcov," %d (%.3f),",(int) age, c12); |
fprintf(fichtmcov," %d (%.3f),",(int) age, c12); |
|
Line 6505 void printinghtml(char fileresu[], char
|
Line 6727 void printinghtml(char fileresu[], char
|
| int lastpass, int stepm, int weightopt, char model[],\ |
int lastpass, int stepm, int weightopt, char model[],\ |
| int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ |
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ |
| int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \ |
int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \ |
| double jprev1, double mprev1,double anprev1, double dateprev1, \ |
double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \ |
| double jprev2, double mprev2,double anprev2, double dateprev2){ |
double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){ |
| int jj1, k1, i1, cpt, k4, nres; |
int jj1, k1, i1, cpt, k4, nres; |
| |
|
| fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
|
Line 6660 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 6882 divided by h: <sub>h</sub>P<sub>ij</sub>
|
| if(prevfcast==1){ |
if(prevfcast==1){ |
| /* Projection of prevalence up to period (stable) prevalence in each health state */ |
/* Projection of prevalence up to period (stable) prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\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. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| <img src=\"%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); |
<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, 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,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \ |
| |
from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \ |
| |
account but can visually be appreciated). 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. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| |
<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); |
| } |
} |
| } |
} |
| |
|
|
Line 6767 true period expectancies (those weighted
|
Line 6999 true period expectancies (those weighted
|
| } |
} |
| |
|
| /******************* Gnuplot file **************/ |
/******************* 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 bage, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){ |
| |
|
| char dirfileres[132],optfileres[132]; |
char dirfileres[132],optfileres[132]; |
| char gplotcondition[132], gplotlabel[132]; |
char gplotcondition[132], gplotlabel[132]; |
|
Line 6776 true period expectancies (those weighted
|
Line 7008 true period expectancies (those weighted
|
| int ng=0; |
int ng=0; |
| int vpopbased; |
int vpopbased; |
| int ioffset; /* variable offset for columns */ |
int ioffset; /* variable offset for columns */ |
| |
int iyearc=1; /* variable column for year of projection */ |
| |
int iagec=1; /* variable column for age of projection */ |
| int nres=0; /* Index of resultline */ |
int nres=0; /* Index of resultline */ |
| int istart=1; /* For starting graphs in projections */ |
int istart=1; /* For starting graphs in projections */ |
| |
|
|
Line 6789 true period expectancies (those weighted
|
Line 7023 true period expectancies (those weighted
|
| /*#endif */ |
/*#endif */ |
| m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
| |
|
| |
/* diagram of the model */ |
| |
fprintf(ficgp,"\n#Diagram of the model \n"); |
| |
fprintf(ficgp,"\ndelta=0.03;delta2=0.07;unset arrow;\n"); |
| |
fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate); |
| |
fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i*10+j from cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.95*(cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) - cos(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta2:0)), -0.95*(sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) - sin(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate); |
| |
|
| |
fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate); |
| |
fprintf(ficgp,"\n#show arrow\nunset label\n"); |
| |
fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)), yoff+sin(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate); |
| |
fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0. font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate); |
| |
fprintf(ficgp,"\n#show label\nunset border;unset xtics; unset ytics;\n"); |
| |
fprintf(ficgp,"\n\nset ter svg size 640, 480;set out \"%s_.svg\" \n",subdirf2(optionfilefiname,"D_")); |
| |
fprintf(ficgp,"unset log y; plot [-1.2:1.2][yoff-1.2:1.2] 1/0 not; set out;reset;\n"); |
| |
|
| /* Contribution to likelihood */ |
/* Contribution to likelihood */ |
| /* Plot the probability implied in the likelihood */ |
/* Plot the probability implied in the likelihood */ |
| fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); |
fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); |
|
Line 6858 true period expectancies (those weighted
|
Line 7106 true period expectancies (those weighted
|
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
| fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); |
fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); |
| fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); |
/* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */ |
| |
fprintf(ficgp,"set title \"Alive state %d %s\" font \"Helvetica,12\"\n",cpt,gplotlabel); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); |
| /* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */ |
/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */ |
| /* k1-1 error should be nres-1*/ |
/* k1-1 error should be nres-1*/ |
|
Line 6880 true period expectancies (those weighted
|
Line 7129 true period expectancies (those weighted
|
| |
|
| fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_")); |
fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_")); |
| if(cptcoveff ==0){ |
if(cptcoveff ==0){ |
| fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3", 2+(cpt-1), cpt ); |
fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3", 2+3*(cpt-1), cpt ); |
| }else{ |
}else{ |
| kl=0; |
kl=0; |
| for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ |
for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ |
|
Line 6931 true period expectancies (those weighted
|
Line 7180 true period expectancies (those weighted
|
| } |
} |
| } /* end covariate */ |
} /* end covariate */ |
| } /* end if no 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 6 dt 3,\"%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 4,\"%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 4"); |
| |
} /* end if backprojcast */ |
| } /* end if backcast */ |
} /* end if backcast */ |
| fprintf(ficgp,"\nset out ;unset label;\n"); |
/* fprintf(ficgp,"\nset out ;unset label;\n"); */ |
| |
fprintf(ficgp,"\nset out ;unset title;\n"); |
| } /* nres */ |
} /* nres */ |
| } /* k1 */ |
} /* k1 */ |
| } /* cpt */ |
} /* cpt */ |
|
Line 7224 set ter svg size 640, 480\nunset log y\n
|
Line 7493 set ter svg size 640, 480\nunset log y\n
|
| for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
| if(m != 1 && TKresult[nres]!= k1) |
if(m != 1 && TKresult[nres]!= k1) |
| continue; |
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,"("); |
strcpy(gplotlabel,"("); |
| fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt); |
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 */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
|
Line 7248 set ter svg size 640, 480\nunset log y\n
|
Line 7517 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,"\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\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
| set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; /* Offset */ |
k=3; /* Offset */ |
| for (i=1; i<= nlstate ; i ++){ /* State of origin */ |
for (i=1; i<= nlstate ; i ++){ /* State of arrival */ |
| if(i==1) |
if(i==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_")); |
| else |
else |
|
Line 7265 set ter svg size 640, 480\nunset log y\n
|
Line 7534 set ter svg size 640, 480\nunset log y\n
|
| /* for (j=2; j<= nlstate ; j ++) */ |
/* for (j=2; j<= nlstate ; j ++) */ |
| /* fprintf(ficgp,"+$%d",k+l+j-1); */ |
/* fprintf(ficgp,"+$%d",k+l+j-1); */ |
| /* /\* 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 */ |
} /* nlstate */ |
| fprintf(ficgp,"\nset out; unset label;\n"); |
fprintf(ficgp,"\nset out; unset label;\n"); |
| } /* end cpt state*/ |
} /* end cpt state*/ |
|
Line 7330 set ter svg size 640, 480\nunset log y\n
|
Line 7599 set ter svg size 640, 480\nunset log y\n
|
| /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
| fprintf(ficgp," u %d:(", ioffset); |
fprintf(ficgp," u %d:(", ioffset); |
| if(i==nlstate+1){ |
if(i==nlstate+1){ |
| fprintf(ficgp," $%d/(1.-$%d)):5 t 'pw.%d' with line lc variable ", \ |
fprintf(ficgp," $%d/(1.-$%d)):1 t 'pw.%d' with line lc variable ", \ |
| 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,cpt ); |
| fprintf(ficgp,",\\\n '' "); |
fprintf(ficgp,",\\\n '' "); |
| fprintf(ficgp," u %d:(",ioffset); |
fprintf(ficgp," u %d:(",ioffset); |
| fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \ |
fprintf(ficgp," (($1-$2) == %d ) ? $%d/(1.-$%d) : 1/0):1 with labels center not ", \ |
| offyear, \ |
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 |
}else |
| fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ |
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 ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| }else{ /* more than 2 covariates */ |
}else{ /* more than 2 covariates */ |
| if(cptcoveff ==1){ |
ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/ |
| ioffset=4; /* Age is in 4 */ |
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| }else{ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| ioffset=6; /* Age is in 6 */ |
iyearc=ioffset-1; |
| /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
iagec=ioffset; |
| /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
|
| } |
|
| fprintf(ficgp," u %d:(",ioffset); |
fprintf(ficgp," u %d:(",ioffset); |
| kl=0; |
kl=0; |
| strcpy(gplotcondition,"("); |
strcpy(gplotcondition,"("); |
|
Line 7369 set ter svg size 640, 480\nunset log y\n
|
Line 7636 set ter svg size 640, 480\nunset log y\n
|
| /*6+1+(i-1)+(nlstate+1)*nlstate; 6+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*/ |
/* '' 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){ |
if(i==nlstate+1){ |
| fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):5 t 'p.%d' with line lc variable", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \ |
| 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,iyearc, cpt ); |
| fprintf(ficgp,",\\\n '' "); |
fprintf(ficgp,",\\\n '' "); |
| fprintf(ficgp," u %d:(",ioffset); |
fprintf(ficgp," u %d:(",iagec); |
| fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ |
fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \ |
| offyear, \ |
iyearc, iagec, 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, iyearc ); |
| /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/ |
/* '' 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{ |
}else{ |
| fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ |
|
Line 7388 set ter svg size 640, 480\nunset log y\n
|
Line 7655 set ter svg size 640, 480\nunset log y\n
|
| } /* end covariate */ |
} /* end covariate */ |
| } /* End if prevfcast */ |
} /* 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)):1 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," (($1-$2) == %d ) ? $%d : 1/0):1 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 */ |
| |
ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/ |
| |
/*# 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 */ |
| |
iyearc=ioffset-1; |
| |
iagec=ioffset; |
| |
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 <cptcoveff && cptcoveff>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):%d t 'bw%d' with line lc variable", gplotcondition, \ |
| |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1),iyearc,cpt ); |
| |
fprintf(ficgp,",\\\n '' "); |
| |
fprintf(ficgp," u %d:(",iagec); |
| |
/* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */ |
| |
fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d : 1/0):%d with labels center not ", gplotcondition, \ |
| |
iyearc,iagec,offbyear, \ |
| |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), iyearc ); |
| |
/* '' 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 */ |
/* 9eme writing MLE parameters */ |
| fprintf(ficgp,"\n##############\n#9eme MLE estimated parameters\n#############\n"); |
fprintf(ficgp,"\n##############\n#9eme MLE estimated parameters\n#############\n"); |
|
Line 7433 set ter svg size 640, 480\nunset log y\n
|
Line 7815 set ter svg size 640, 480\nunset log y\n
|
| continue; |
continue; |
| fprintf(ficgp,"\n\n# Combination of dummy k1=%d which is ",k1); |
fprintf(ficgp,"\n\n# Combination of dummy k1=%d which is ",k1); |
| strcpy(gplotlabel,"("); |
strcpy(gplotlabel,"("); |
| sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1); |
/*sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);*/ |
| for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ |
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 */ |
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,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
|
Line 7450 set ter svg size 640, 480\nunset log y\n
|
Line 7832 set ter svg size 640, 480\nunset log y\n
|
| strcpy(gplotlabel+strlen(gplotlabel),")"); |
strcpy(gplotlabel+strlen(gplotlabel),")"); |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres); |
| fprintf(ficgp,"\nset label \"%s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel); |
fprintf(ficgp,"\nset key outside "); |
| |
/* fprintf(ficgp,"\nset label \"%s\" at graph 1.2,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel); */ |
| |
fprintf(ficgp,"\nset title \"%s\" font \"Helvetica,12\"\n",gplotlabel); |
| fprintf(ficgp,"\nset ter svg size 640, 480 "); |
fprintf(ficgp,"\nset ter svg size 640, 480 "); |
| if (ng==1){ |
if (ng==1){ |
| fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ |
fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ |
|
Line 7494 set ter svg size 640, 480\nunset log y\n
|
Line 7878 set ter svg size 640, 480\nunset log y\n
|
| /* for(j=3; j <=ncovmodel-nagesqr; j++) { */ |
/* for(j=3; j <=ncovmodel-nagesqr; j++) { */ |
| for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */ |
for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */ |
| /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ |
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ |
| if(j==Tage[ij]) { /* Product by age */ |
if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ |
| if(ij <=cptcovage) { /* 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(DummyV[j]==0){ |
if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ |
| fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; |
if(DummyV[j]==0){ |
| }else{ /* quantitative */ |
fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; |
| fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */ |
}else{ /* quantitative */ |
| /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ |
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]]); |
|
| } |
} |
| }else{ /* Vn*Vm Vn is quanti */ |
ij++; |
| 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 */ |
}else if(cptcovprod >0){ |
| fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); |
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 */ |
} else{ /* simple covariate */ |
| /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */ |
/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */ |
| if(Dummy[j]==0){ |
if(Dummy[j]==0){ |
|
Line 7551 set ter svg size 640, 480\nunset log y\n
|
Line 7939 set ter svg size 640, 480\nunset log y\n
|
| |
|
| ij=1; |
ij=1; |
| for(j=3; j <=ncovmodel-nagesqr; j++){ |
for(j=3; j <=ncovmodel-nagesqr; j++){ |
| if((j-2)==Tage[ij]) { /* Bug valgrind */ |
if(cptcovage >0){ |
| if(ij <=cptcovage) { /* Bug valgrind */ |
if((j-2)==Tage[ij]) { /* Bug valgrind */ |
| fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]); |
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,Tvar[j-2])]); */ |
fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]); |
| ij++; |
/* 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 */ |
}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,")"); |
fprintf(ficgp,")"); |
| } |
} |
| fprintf(ficgp,")"); |
fprintf(ficgp,")"); |
| if(ng ==2) |
if(ng ==2) |
| fprintf(ficgp," t \"p%d%d\" ", k2,k); |
fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); |
| else /* ng= 3 */ |
else /* ng= 3 */ |
| fprintf(ficgp," t \"i%d%d\" ", k2,k); |
fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); |
| }else{ /* end ng <> 1 */ |
}else{ /* end ng <> 1 */ |
| if( k !=k2) /* logit p11 is hard to draw */ |
if( k !=k2) /* logit p11 is hard to draw */ |
| fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k); |
fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); |
| } |
} |
| if ((k+k2)!= (nlstate*2+ndeath) && ng != 1) |
if ((k+k2)!= (nlstate*2+ndeath) && ng != 1) |
| fprintf(ficgp,","); |
fprintf(ficgp,","); |
|
Line 7579 set ter svg size 640, 480\nunset log y\n
|
Line 7968 set ter svg size 640, 480\nunset log y\n
|
| i=i+ncovmodel; |
i=i+ncovmodel; |
| } /* end k */ |
} /* end k */ |
| } /* end k2 */ |
} /* end k2 */ |
| fprintf(ficgp,"\n set out; unset label;\n"); |
/* fprintf(ficgp,"\n set out; unset label;set key default;\n"); */ |
| |
fprintf(ficgp,"\n set out; unset title;set key default;\n"); |
| } /* end k1 */ |
} /* end k1 */ |
| } /* end ng */ |
} /* end ng */ |
| /* avoid: */ |
/* avoid: */ |
|
Line 7603 set ter svg size 640, 480\nunset log y\n
|
Line 7993 set ter svg size 640, 480\nunset log y\n
|
| double *agemingoodr, *agemaxgoodr; |
double *agemingoodr, *agemaxgoodr; |
| |
|
| |
|
| /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ |
/* modcovmax=2*cptcoveff; Max number of modalities. We suppose */ |
| /* a covariate has 2 modalities, should be equal to ncovcombmax *\/ */ |
/* a covariate has 2 modalities, should be equal to ncovcombmax */ |
| |
|
| sumnewp = vector(1,ncovcombmax); |
sumnewp = vector(1,ncovcombmax); |
| sumnewm = vector(1,ncovcombmax); |
sumnewm = vector(1,ncovcombmax); |
|
Line 7724 set ter svg size 640, 480\nunset log y\n
|
Line 8114 set ter svg size 640, 480\nunset log y\n
|
| sumr+=probs[(int)age][i][cptcod]; |
sumr+=probs[(int)age][i][cptcod]; |
| } |
} |
| if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
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 */ |
} /* end bad */ |
| /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
| if(fabs(sumr - 1.) > 1.e-3) { /* bad */ |
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 */ |
} /* end bad */ |
| }/* age */ |
}/* age */ |
| |
|
|
Line 7764 set ter svg size 640, 480\nunset log y\n
|
Line 8154 set ter svg size 640, 480\nunset log y\n
|
| sumr+=mobaverage[(int)age][i][cptcod]; |
sumr+=mobaverage[(int)age][i][cptcod]; |
| } |
} |
| if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
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 */ |
} /* end bad */ |
| /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
| if(fabs(sumr - 1.) > 1.e-3) { /* bad */ |
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 */ |
} /* end bad */ |
| }/* age */ |
}/* age */ |
| |
|
|
Line 7817 set ter svg size 640, 480\nunset log y\n
|
Line 8207 set ter svg size 640, 480\nunset log y\n
|
| |
|
| |
|
| /************** Forecasting ******************/ |
/************** 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 ***prev, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){ |
| /* proj1, year, month, day of starting projection |
/* proj1, year, month, day of starting projection |
| agemin, agemax range of age |
agemin, agemax range of age |
| dateprev1 dateprev2 range of dates during which prevalence is computed |
dateprev1 dateprev2 range of dates during which prevalence is computed |
|
Line 7856 void prevforecast(char fileres[], double
|
Line 8246 void prevforecast(char fileres[], double
|
| if(estepm < stepm){ |
if(estepm < stepm){ |
| printf ("Problem %d lower than %d\n",estepm, stepm); |
printf ("Problem %d lower than %d\n",estepm, stepm); |
| } |
} |
| else hstepm=estepm; |
else{ |
| |
hstepm=estepm; |
| |
} |
| |
if(estepm > stepm){ /* Yes every two year */ |
| |
stepsize=2; |
| |
} |
| |
|
| hstepm=hstepm/stepm; |
hstepm=hstepm/stepm; |
| yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and |
yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and |
|
Line 7901 void prevforecast(char fileres[], double
|
Line 8296 void prevforecast(char fileres[], double
|
| for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { |
for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { |
| fprintf(ficresf,"\n"); |
fprintf(ficresf,"\n"); |
| fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); |
fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); |
| for (agec=fage; agec>=(ageminpar-1); agec--){ |
/* for (agec=fage; agec>=(ageminpar-1); agec--){ */ |
| |
for (agec=fage; agec>=(bage); agec--){ |
| nhstepm=(int) rint((agelim-agec)*YEARM/stepm); |
nhstepm=(int) rint((agelim-agec)*YEARM/stepm); |
| nhstepm = nhstepm/hstepm; |
nhstepm = nhstepm/hstepm; |
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| |
/* We compute pii at age agec over nhstepm);*/ |
| hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, 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 */ |
| for (h=0; h<=nhstepm; h++){ |
for (h=0; h<=nhstepm; h++){ |
| if (h*hstepm/YEARM*stepm ==yearp) { |
if (h*hstepm/YEARM*stepm ==yearp) { |
| fprintf(ficresf,"\n"); |
break; |
| 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); |
fprintf(ficresf,"\n"); |
| } |
for(j=1;j<=cptcoveff;j++) |
| for(j=1; j<=nlstate+ndeath;j++) { |
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| ppij=0.; |
fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); |
| for(i=1; i<=nlstate;i++) { |
|
| /* if (mobilav>=1) */ |
for(j=1; j<=nlstate+ndeath;j++) { |
| ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k]; |
ppij=0.; |
| /* else { */ /* even if mobilav==-1 we use mobaverage */ |
for(i=1; i<=nlstate;i++) { |
| /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ |
if (mobilav>=1) |
| /* } */ |
ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k]; |
| if (h*hstepm/YEARM*stepm== yearp) { |
else { /* even if mobilav==-1 we use mobaverage, probs may not sums to 1 */ |
| fprintf(ficresf," %.3f", p3mat[i][j][h]); |
ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; |
| } |
|
| } /* end i */ |
|
| if (h*hstepm/YEARM*stepm==yearp) { |
|
| fprintf(ficresf," %.3f", ppij); |
|
| } |
} |
| }/* end j */ |
fprintf(ficresf," %.3f", p3mat[i][j][h]); |
| } /* end h */ |
} /* end i */ |
| |
fprintf(ficresf," %.3f", ppij); |
| |
}/* end j */ |
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| } /* end agec */ |
} /* end agec */ |
| /* diffyear=(int) anproj1+yearp-ageminpar-1; */ |
/* diffyear=(int) anproj1+yearp-ageminpar-1; */ |
|
Line 7945 void prevforecast(char fileres[], double
|
Line 8340 void prevforecast(char fileres[], double
|
| |
|
| } |
} |
| |
|
| /* /\************** Back Forecasting ******************\/ */ |
/************** Back Forecasting ******************/ |
| 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){ |
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 |
/* back1, year, month, day of starting backection |
| agemin, agemax range of age |
agemin, agemax range of age |
| dateprev1 dateprev2 range of dates during which prevalence is computed |
dateprev1 dateprev2 range of dates during which prevalence is computed |
| anback2 year of en of backection (same day and month as back1). |
anback2 year of end of backprojection (same day and month as back1). |
| |
prevacurrent and prev are prevalences. |
| */ |
*/ |
| 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 agec; /* generic age */ |
| double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; |
double agelim, ppij, ppi, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; |
| double *popeffectif,*popcount; |
double *popeffectif,*popcount; |
| double ***p3mat; |
double ***p3mat; |
| /* double ***mobaverage; */ |
/* double ***mobaverage; */ |
| char fileresfb[FILENAMELENGTH]; |
char fileresfb[FILENAMELENGTH]; |
| |
|
| agelim=AGESUP; |
agelim=AGEINF; |
| /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people |
/* 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). |
in each health status at the date of interview (if between dateprev1 and dateprev2). |
| We still use firstpass and lastpass as another selection. |
We still use firstpass and lastpass as another selection. |
|
Line 7989 void prevbackforecast(char fileres[], do
|
Line 8385 void prevbackforecast(char fileres[], do
|
| if(estepm < stepm){ |
if(estepm < stepm){ |
| printf ("Problem %d lower than %d\n",estepm, stepm); |
printf ("Problem %d lower than %d\n",estepm, stepm); |
| } |
} |
| else hstepm=estepm; |
else{ |
| |
hstepm=estepm; |
| |
} |
| |
if(estepm >= stepm){ /* Yes every two year */ |
| |
stepsize=2; |
| |
} |
| |
|
| hstepm=hstepm/stepm; |
hstepm=hstepm/stepm; |
| yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and |
yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and |
|
Line 8006 void prevbackforecast(char fileres[], do
|
Line 8407 void prevbackforecast(char fileres[], do
|
| if (cptcovn < 1){i1=1;} |
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); |
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"); |
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(nres=1; nres <= nresult; nres++) /* For each resultline */ |
| for(k=1; k<=i1;k++){ |
for(k=1; k<=i1;k++){ |
| if(i1 != 1 && TKresult[nres]!= k) |
if(i1 != 1 && TKresult[nres]!= k) |
|
Line 8021 void prevbackforecast(char fileres[], do
|
Line 8419 void prevbackforecast(char fileres[], do
|
| printf("\nCombination (%d) projection ignored because no cases \n",k); |
printf("\nCombination (%d) projection ignored because no cases \n",k); |
| continue; |
continue; |
| } |
} |
| fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#"); |
fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#"); |
| for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
| fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| } |
} |
|
Line 8031 void prevbackforecast(char fileres[], do
|
Line 8429 void prevbackforecast(char fileres[], do
|
| fprintf(ficresfb," yearbproj age"); |
fprintf(ficresfb," yearbproj age"); |
| for(j=1; j<=nlstate+ndeath;j++){ |
for(j=1; j<=nlstate+ndeath;j++){ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| fprintf(ficresfb," p%d%d",i,j); |
fprintf(ficresfb," b%d%d",i,j); |
| fprintf(ficresfb," p.%d",j); |
fprintf(ficresfb," b.%d",j); |
| } |
} |
| for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) { |
for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) { |
| /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */ |
/* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */ |
| fprintf(ficresfb,"\n"); |
fprintf(ficresfb,"\n"); |
| fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); |
fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); |
| /* for (agec=fage; agec>=(ageminpar-1); agec--){ */ |
/* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */ |
| /* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */ |
/* for (agec=bage; agec<=agemax-1; agec++){ /\* testing *\/ */ |
| for (agec=fage; agec>=fage-20; agec--){ /* testing up to 10 */ |
for (agec=bage; agec<=fage; agec++){ /* testing */ |
| nhstepm=(int) rint((agelim-agec)*YEARM/stepm); |
/* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/ |
| |
nhstepm=(int) (agec-agelim) *YEARM/stepm;/* nhstepm=(int) rint((agec-agelim)*YEARM/stepm);*/ |
| nhstepm = nhstepm/hstepm; |
nhstepm = nhstepm/hstepm; |
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| |
/* computes hbxij at age agec over 1 to nhstepm */ |
| |
/* printf("####prevbackforecast debug agec=%.2f nhstepm=%d\n",agec, nhstepm);fflush(stdout); */ |
| hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres); |
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++){ |
for (h=0; h<=nhstepm; h++){ |
| if (h*hstepm/YEARM*stepm ==yearp) { |
if (h*hstepm/YEARM*stepm ==-yearp) { |
| fprintf(ficresfb,"\n"); |
break; |
| 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); |
fprintf(ficresfb,"\n"); |
| } |
for(j=1;j<=cptcoveff;j++) |
| for(j=1; j<=nlstate+ndeath;j++) { |
fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| ppij=0.; |
fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec-h*hstepm/YEARM*stepm); |
| for(i=1; i<=nlstate;i++) { |
for(i=1; i<=nlstate+ndeath;i++) { |
| /* if (mobilav==1) */ |
ppij=0.;ppi=0.; |
| ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k]; |
for(j=1; j<=nlstate;j++) { |
| |
/* if (mobilav==1) */ |
| |
ppij=ppij+p3mat[i][j][h]*prevacurrent[(int)agec][j][k]; |
| |
ppi=ppi+prevacurrent[(int)agec][j][k]; |
| |
/* ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k]; */ |
| |
/* ppi=ppi+mobaverage[(int)agec][j][k]; */ |
| /* else { */ |
/* else { */ |
| /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ |
/* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ |
| /* } */ |
/* } */ |
| if (h*hstepm/YEARM*stepm== yearp) { |
fprintf(ficresfb," %.3f", p3mat[i][j][h]); |
| fprintf(ficresfb," %.3f", p3mat[i][j][h]); |
} /* end j */ |
| } |
if(ppi <0.99){ |
| } /* end i */ |
printf("Error in prevbackforecast, prevalence doesn't sum to 1 for state %d: %3f\n",i, ppi); |
| if (h*hstepm/YEARM*stepm==yearp) { |
fprintf(ficlog,"Error in prevbackforecast, prevalence doesn't sum to 1 for state %d: %3f\n",i, ppi); |
| fprintf(ficresfb," %.3f", ppij); |
} |
| } |
fprintf(ficresfb," %.3f", ppij); |
| }/* end j */ |
}/* end j */ |
| } /* end h */ |
|
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| } /* end agec */ |
} /* end agec */ |
| } /* end yearp */ |
} /* end yearp */ |
|
Line 8084 void prevbackforecast(char fileres[], do
|
Line 8491 void prevbackforecast(char fileres[], do
|
| |
|
| } |
} |
| |
|
| |
/* Variance of prevalence limit: varprlim */ |
| |
void varprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **prlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){ |
| |
/*------- Variance of period (stable) prevalence------*/ |
| |
|
| |
char fileresvpl[FILENAMELENGTH]; |
| |
FILE *ficresvpl; |
| |
double **oldm, **savm; |
| |
double **varpl; /* Variances of prevalence limits by age */ |
| |
int i1, k, nres, j ; |
| |
|
| |
strcpy(fileresvpl,"VPL_"); |
| |
strcat(fileresvpl,fileresu); |
| |
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { |
| |
printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); |
| |
exit(0); |
| |
} |
| |
printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); |
| |
fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);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(ficresvpl,"\n#****** "); |
| |
printf("\n#****** "); |
| |
fprintf(ficlog,"\n#****** "); |
| |
for(j=1;j<=cptcoveff;j++) { |
| |
fprintf(ficresvpl,"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(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); |
| |
fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); |
| |
} |
| |
fprintf(ficresvpl,"******\n"); |
| |
printf("******\n"); |
| |
fprintf(ficlog,"******\n"); |
| |
|
| |
varpl=matrix(1,nlstate,(int) bage, (int) fage); |
| |
oldm=oldms;savm=savms; |
| |
varprevlim(fileresvpl, ficresvpl, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyearp, k, strstart, nres); |
| |
free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
| |
/*}*/ |
| |
} |
| |
|
| |
fclose(ficresvpl); |
| |
printf("done variance-covariance of period prevalence\n");fflush(stdout); |
| |
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); |
| |
|
| |
} |
| |
/* Variance of back prevalence: varbprlim */ |
| |
void varbprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **bprlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){ |
| |
/*------- Variance of back (stable) prevalence------*/ |
| |
|
| |
char fileresvbl[FILENAMELENGTH]; |
| |
FILE *ficresvbl; |
| |
|
| |
double **oldm, **savm; |
| |
double **varbpl; /* Variances of back prevalence limits by age */ |
| |
int i1, k, nres, j ; |
| |
|
| |
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); |
| |
|
| |
|
| |
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(fileresvbl, ficresvbl, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, ncvyearp, 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); |
| |
|
| |
} /* End of varbprlim */ |
| |
|
| /************** Forecasting *****not tested NB*************/ |
/************** 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){ */ |
/* 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){ */ |
| |
|
|
Line 10140 int main(int argc, char *argv[])
|
Line 10664 int main(int argc, char *argv[])
|
| int vpopbased=0; |
int vpopbased=0; |
| int nres=0; |
int nres=0; |
| int endishere=0; |
int endishere=0; |
| |
int noffset=0; |
| |
int ncurrv=0; /* Temporary variable */ |
| |
|
| char ca[32], cb[32]; |
char ca[32], cb[32]; |
| /* FILE *fichtm; *//* Html File */ |
/* FILE *fichtm; *//* Html File */ |
| /* FILE *ficgp;*/ /*Gnuplot File */ |
/* FILE *ficgp;*/ /*Gnuplot File */ |
|
Line 10192 int main(int argc, char *argv[])
|
Line 10718 int main(int argc, char *argv[])
|
| double *delti; /* Scale */ |
double *delti; /* Scale */ |
| double ***eij, ***vareij; |
double ***eij, ***vareij; |
| double **varpl; /* Variances of prevalence limits by age */ |
double **varpl; /* Variances of prevalence limits by age */ |
| |
|
| double *epj, vepp; |
double *epj, vepp; |
| |
|
| double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; |
double dateprev1, dateprev2; |
| double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000; |
double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0; |
| |
double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0; |
| |
|
| double **ximort; |
double **ximort; |
| char *alph[]={"a","a","b","c","d","e"}, str[4]="1234"; |
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234"; |
|
Line 10273 int main(int argc, char *argv[])
|
Line 10801 int main(int argc, char *argv[])
|
| if(pathr[0] == '\0') break; /* Dirty */ |
if(pathr[0] == '\0') break; /* Dirty */ |
| } |
} |
| } |
} |
| |
else if (argc<=2){ |
| |
strcpy(pathtot,argv[1]); |
| |
} |
| else{ |
else{ |
| strcpy(pathtot,argv[1]); |
strcpy(pathtot,argv[1]); |
| |
strcpy(z,argv[2]); |
| |
printf("\nargv[2]=%s z=%c\n",argv[2],z[0]); |
| } |
} |
| /*if(getcwd(pathcd, MAXLINE)!= NULL)printf ("Error pathcd\n");*/ |
/*if(getcwd(pathcd, MAXLINE)!= NULL)printf ("Error pathcd\n");*/ |
| /*cygwin_split_path(pathtot,path,optionfile); |
/*cygwin_split_path(pathtot,path,optionfile); |
|
Line 10352 int main(int argc, char *argv[])
|
Line 10885 int main(int argc, char *argv[])
|
| exit(70); |
exit(70); |
| } |
} |
| |
|
| |
|
| |
|
| strcpy(filereso,"o"); |
strcpy(filereso,"o"); |
| strcat(filereso,fileresu); |
strcat(filereso,fileresu); |
| if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */ |
if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */ |
|
Line 10362 int main(int argc, char *argv[])
|
Line 10893 int main(int argc, char *argv[])
|
| fflush(ficlog); |
fflush(ficlog); |
| goto end; |
goto end; |
| } |
} |
| |
/*-------- Rewriting parameter file ----------*/ |
| |
strcpy(rfileres,"r"); /* "Rparameterfile */ |
| |
strcat(rfileres,optionfilefiname); /* Parameter file first name */ |
| |
strcat(rfileres,"."); /* */ |
| |
strcat(rfileres,optionfilext); /* Other files have txt extension */ |
| |
if((ficres =fopen(rfileres,"w"))==NULL) { |
| |
printf("Problem writing new parameter file: %s\n", rfileres);goto end; |
| |
fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end; |
| |
fflush(ficlog); |
| |
goto end; |
| |
} |
| |
fprintf(ficres,"#IMaCh %s\n",version); |
| |
|
| |
|
| /* Reads comments: lines beginning with '#' */ |
/* Reads comments: lines beginning with '#' */ |
| numlinepar=0; |
numlinepar=0; |
| |
/* Is it a BOM UTF-8 Windows file? */ |
| /* First parameter line */ |
/* First parameter line */ |
| while(fgets(line, MAXLINE, ficpar)) { |
while(fgets(line, MAXLINE, ficpar)) { |
| |
noffset=0; |
| |
if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */ |
| |
{ |
| |
noffset=noffset+3; |
| |
printf("# File is an UTF8 Bom.\n"); // 0xBF |
| |
} |
| |
else if( line[0] == (char)0xFE && line[1] == (char)0xFF) |
| |
{ |
| |
noffset=noffset+2; |
| |
printf("# File is an UTF16BE BOM file\n"); |
| |
} |
| |
else if( line[0] == 0 && line[1] == 0) |
| |
{ |
| |
if( line[2] == (char)0xFE && line[3] == (char)0xFF){ |
| |
noffset=noffset+4; |
| |
printf("# File is an UTF16BE BOM file\n"); |
| |
} |
| |
} else{ |
| |
;/*printf(" Not a BOM file\n");*/ |
| |
} |
| |
|
| /* If line starts with a # it is a comment */ |
/* If line starts with a # it is a comment */ |
| if (line[0] == '#') { |
if (line[noffset] == '#') { |
| numlinepar++; |
numlinepar++; |
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| |
fputs(line,ficres); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| continue; |
continue; |
| }else |
}else |
|
Line 10393 int main(int argc, char *argv[])
|
Line 10959 int main(int argc, char *argv[])
|
| numlinepar++; |
numlinepar++; |
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| |
fputs(line,ficres); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| continue; |
continue; |
| }else |
}else |
|
Line 10415 int main(int argc, char *argv[])
|
Line 10982 int main(int argc, char *argv[])
|
| numlinepar++; |
numlinepar++; |
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| |
fputs(line,ficres); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| continue; |
continue; |
| }else |
}else |
| break; |
break; |
| } |
} |
| if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ |
if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ |
| if (num_filled == 0){ |
if (num_filled != 1){ |
| printf("ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line); |
printf("ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line); |
| fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line); |
fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line); |
| model[0]='\0'; |
|
| goto end; |
|
| } else if (num_filled != 1){ |
|
| printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); |
|
| fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); |
|
| model[0]='\0'; |
model[0]='\0'; |
| goto end; |
goto end; |
| } |
} |
|
Line 10450 int main(int argc, char *argv[])
|
Line 11013 int main(int argc, char *argv[])
|
| fflush(ficlog); |
fflush(ficlog); |
| /* if(model[0]=='#'|| model[0]== '\0'){ */ |
/* if(model[0]=='#'|| model[0]== '\0'){ */ |
| if(model[0]=='#'){ |
if(model[0]=='#'){ |
| printf("Error in 'model' line: model should start with 'model=1+age+' and end with '.' \n \ |
printf("Error in 'model' line: model should start with 'model=1+age+' and end without space \n \ |
| 'model=1+age+.' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age.' or \n \ |
'model=1+age+' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age' or \n \ |
| 'model=1+age+V1+V2.' or 'model=1+age+V1+V2+V1*V2.' etc. \n"); \ |
'model=1+age+V1+V2' or 'model=1+age+V1+V2+V1*V2' etc. \n"); \ |
| if(mle != -1){ |
if(mle != -1){ |
| printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter file.\n"); |
printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter vectors and subdiagonal covariance matrix.\n"); |
| exit(1); |
exit(1); |
| } |
} |
| } |
} |
|
Line 10475 int main(int argc, char *argv[])
|
Line 11038 int main(int argc, char *argv[])
|
| |
|
| |
|
| covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ |
covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ |
| coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */ |
if(nqv>=1)coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */ |
| cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/ |
if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */ |
| 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*/ |
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+v3+v2*v4+v5*age makes cptcovn = 5 |
| v1+v2*age+v2*v3 makes cptcovn = 3 |
v1+v2*age+v2*v3 makes cptcovn = 3 |
|
Line 10676 Please run with mle=-1 to get a correct
|
Line 11239 Please run with mle=-1 to get a correct
|
| |
|
| fflush(ficlog); |
fflush(ficlog); |
| |
|
| /*-------- Rewriting parameter file ----------*/ |
|
| strcpy(rfileres,"r"); /* "Rparameterfile */ |
|
| strcat(rfileres,optionfilefiname); /* Parameter file first name*/ |
|
| strcat(rfileres,"."); /* */ |
|
| strcat(rfileres,optionfilext); /* Other files have txt extension */ |
|
| if((ficres =fopen(rfileres,"w"))==NULL) { |
|
| printf("Problem writing new parameter file: %s\n", rfileres);goto end; |
|
| fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end; |
|
| } |
|
| fprintf(ficres,"#%s\n",version); |
|
| } /* End of mle != -3 */ |
} /* End of mle != -3 */ |
| |
|
| /* Main data |
/* Main data |
|
Line 11023 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 11576 Title=%s <br>Datafile=%s Firstpass=%d La
|
| firstpass, lastpass, stepm, weightopt, model); |
firstpass, lastpass, stepm, weightopt, model); |
| |
|
| fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
| fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
fprintf(fichtm,"<h4>Parameter line 2</h4><ul><li>Tolerance for the convergence of the likelihood: ftol=%f \n<li>Interval for the elementary matrix (in month): stepm=%d",\ |
| |
ftol, stepm); |
| |
fprintf(fichtm,"\n<li>Number of fixed dummy covariates: ncovcol=%d ", ncovcol); |
| |
ncurrv=1; |
| |
for(i=ncurrv; i <=ncovcol; i++) fprintf(fichtm,"V%d ", i); |
| |
fprintf(fichtm,"\n<li> Number of fixed quantitative variables: nqv=%d ", nqv); |
| |
ncurrv=i; |
| |
for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i); |
| |
fprintf(fichtm,"\n<li> Number of time varying (wave varying) covariates: ntv=%d ", ntv); |
| |
ncurrv=i; |
| |
for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i); |
| |
fprintf(fichtm,"\n<li>Number of quantitative time varying covariates: nqtv=%d ", nqtv); |
| |
ncurrv=i; |
| |
for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i); |
| |
fprintf(fichtm,"\n<li>Weights column \n<br>Number of alive states: nlstate=%d <br>Number of death states (not really implemented): ndeath=%d \n<li>Number of waves: maxwav=%d \n<li>Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n<li>Does the weight column be taken into account (1), or not (0): weight=%d</ul>\n", \ |
| |
nlstate, ndeath, maxwav, mle, weightopt); |
| |
|
| |
fprintf(fichtm,"<h4> Diagram of states <a href=\"%s_.svg\">%s_.svg</a></h4> \n\ |
| |
<img src=\"%s_.svg\">", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_")); |
| |
|
| |
|
| |
fprintf(fichtm,"\n<h4>Some descriptive statistics </h4>\n<br>Total number of observations=%d <br>\n\ |
| Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\ |
Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\ |
| Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\ |
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\ |
| imx,agemin,agemax,jmin,jmax,jmean); |
imx,agemin,agemax,jmin,jmax,jmean); |
| pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
| oldms= 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 */ |
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
| savms= 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 */ |
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ |
| |
|
| /* For Powell, parameters are in a vector p[] starting at p[1] |
/* 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] */ |
so we point p on param[1][1] so that p[1] maps on param[1][1][1] */ |
|
Line 11586 Please run with mle=-1 to get a correct
|
Line 12160 Please run with mle=-1 to get a correct
|
| fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
| fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
| /* day and month of proj2 are not used but only year anproj2.*/ |
/* day and month of proj2 are not used but only year anproj2.*/ |
| |
dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.; |
| |
dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.; |
| |
|
| } |
} |
| break; |
break; |
| case 12: |
case 12: |
|
Line 11601 Please run with mle=-1 to get a correct
|
Line 12178 Please run with mle=-1 to get a correct
|
| fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| /* day and month of proj2 are not used but only year anproj2.*/ |
/* day and month of proj2 are not used but only year anproj2.*/ |
| |
dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.; |
| |
dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.; |
| } |
} |
| break; |
break; |
| case 13: |
case 13: |
|
Line 11651 Please run with mle=-1 to get a correct
|
Line 12230 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\ |
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); |
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
| }else{ |
}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); */ |
| |
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, backcast, pathc,p, (int)anproj1-bage, (int)anback1-fage); |
| } |
} |
| printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \ |
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \ |
| model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \ |
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \ |
| jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2); |
jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2); |
| |
|
| /*------------ free_vector -------------*/ |
/*------------ free_vector -------------*/ |
| /* chdir(path); */ |
/* chdir(path); */ |
|
Line 11691 Please run with mle=-1 to get a correct
|
Line 12271 Please run with mle=-1 to get a correct
|
| k=1; |
k=1; |
| varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart); |
varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart); |
| |
|
| /* Prevalence for each covariates in probs[age][status][cov] */ |
/* Prevalence for each covariate combination in probs[age][status][cov] */ |
| probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
probs= ma3x(AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
| for(i=1;i<=AGESUP;i++) |
for(i=AGEINF;i<=AGESUP;i++) |
| for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */ |
for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */ |
| for(k=1;k<=ncovcombmax;k++) |
for(k=1;k<=ncovcombmax;k++) |
| probs[i][j][k]=0.; |
probs[i][j][k]=0.; |
| prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, |
| |
ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
| if (mobilav!=0 ||mobilavproj !=0 ) { |
if (mobilav!=0 ||mobilavproj !=0 ) { |
| mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
mobaverages= ma3x(AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
| for(i=1;i<=AGESUP;i++) |
for(i=AGEINF;i<=AGESUP;i++) |
| for(j=1;j<=nlstate;j++) |
for(j=1;j<=nlstate+ndeath;j++) |
| for(k=1;k<=ncovcombmax;k++) |
for(k=1;k<=ncovcombmax;k++) |
| mobaverages[i][j][k]=0.; |
mobaverages[i][j][k]=0.; |
| mobaverage=mobaverages; |
mobaverage=mobaverages; |
|
Line 11712 Please run with mle=-1 to get a correct
|
Line 12293 Please run with mle=-1 to get a correct
|
| fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); |
fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); |
| printf(" Error in movingaverage mobilav=%d\n",mobilav); |
printf(" Error in movingaverage mobilav=%d\n",mobilav); |
| } |
} |
| } |
} else if (mobilavproj !=0) { |
| /* else if(mobilavproj==-1){ /\* Forcing raw observed prevalences *\/ */ |
|
| /* for(i=1;i<=AGESUP;i++) */ |
|
| /* for(j=1;j<=nlstate;j++) */ |
|
| /* for(k=1;k<=ncovcombmax;k++) */ |
|
| /* mobaverages[i][j][k]=probs[i][j][k]; */ |
|
| /* /\* /\\* Prevalence for each covariates in probs[age][status][cov] *\\/ *\/ */ |
|
| /* /\* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); *\/ */ |
|
| /* } */ |
|
| else if (mobilavproj !=0) { |
|
| printf("Movingaveraging projected observed prevalence\n"); |
printf("Movingaveraging projected observed prevalence\n"); |
| fprintf(ficlog,"Movingaveraging projected observed prevalence\n"); |
fprintf(ficlog,"Movingaveraging projected observed prevalence\n"); |
| if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){ |
if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){ |
| fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj); |
fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj); |
| printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj); |
printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj); |
| } |
} |
| |
}else{ |
| |
printf("Internal error moving average\n"); |
| |
fflush(stdout); |
| |
exit(1); |
| } |
} |
| }/* end if moving average */ |
}/* end if moving average */ |
| |
|
| /*---------- Forecasting ------------------*/ |
/*---------- Forecasting ------------------*/ |
| /*if((stepm == 1) && (strcmp(model,".")==0)){*/ |
|
| if(prevfcast==1){ |
if(prevfcast==1){ |
| /* if(stepm ==1){*/ |
/* if(stepm ==1){*/ |
| prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); |
prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); |
| } |
} |
| |
|
| |
/* Backcasting */ |
| if(backcast==1){ |
if(backcast==1){ |
| ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
| ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
|
Line 11745 Please run with mle=-1 to get a correct
|
Line 12322 Please run with mle=-1 to get a correct
|
| /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Back Prevalence limit (period or stable prevalence) --------------*/ |
| |
|
| bprlim=matrix(1,nlstate,1,nlstate); |
bprlim=matrix(1,nlstate,1,nlstate); |
| |
|
| back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); |
back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); |
| fclose(ficresplb); |
fclose(ficresplb); |
| |
|
| hBijx(p, bage, fage, mobaverage); |
hBijx(p, bage, fage, mobaverage); |
| fclose(ficrespijb); |
fclose(ficrespijb); |
| free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
|
| |
|
| prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, |
prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, |
| bage, fage, firstpass, lastpass, anback2, p, cptcoveff); |
mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); |
| |
varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); |
| |
|
| |
|
| |
free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
| free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath); |
free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath); |
| free_matrix(ddsavms, 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(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath); |
| } |
} /* end Backcasting */ |
| |
|
| |
|
| /* ------ Other prevalence ratios------------ */ |
/* ------ Other prevalence ratios------------ */ |
| |
|
|
Line 11810 Please run with mle=-1 to get a correct
|
Line 12391 Please run with mle=-1 to get a correct
|
| fclose(ficreseij); |
fclose(ficreseij); |
| printf("done evsij\n");fflush(stdout); |
printf("done evsij\n");fflush(stdout); |
| fprintf(ficlog,"done evsij\n");fflush(ficlog); |
fprintf(ficlog,"done evsij\n");fflush(ficlog); |
| |
|
| |
|
| /*---------- State-specific expectancies and variances ------------*/ |
/*---------- State-specific expectancies and variances ------------*/ |
| |
|
| |
|
| strcpy(filerest,"T_"); |
strcpy(filerest,"T_"); |
| strcat(filerest,fileresu); |
strcat(filerest,fileresu); |
| if((ficrest=fopen(filerest,"w"))==NULL) { |
if((ficrest=fopen(filerest,"w"))==NULL) { |
|
Line 11822 Please run with mle=-1 to get a correct
|
Line 12403 Please run with mle=-1 to get a correct
|
| } |
} |
| printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout); |
printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout); |
| fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog); |
fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog); |
| |
|
| |
|
| strcpy(fileresstde,"STDE_"); |
strcpy(fileresstde,"STDE_"); |
| strcat(fileresstde,fileresu); |
strcat(fileresstde,fileresu); |
| if((ficresstdeij=fopen(fileresstde,"w"))==NULL) { |
if((ficresstdeij=fopen(fileresstde,"w"))==NULL) { |
|
Line 11851 Please run with mle=-1 to get a correct
|
Line 12430 Please run with mle=-1 to get a correct
|
| printf(" Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout); |
printf(" Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout); |
| fprintf(ficlog," Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog); |
fprintf(ficlog," Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog); |
| |
|
| /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
|
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
|
| |
|
| i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */ |
i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */ |
| if (cptcovn < 1){i1=1;} |
if (cptcovn < 1){i1=1;} |
| |
|
|
Line 11915 Please run with mle=-1 to get a correct
|
Line 12491 Please run with mle=-1 to get a correct
|
| vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
| pstamp(ficrest); |
pstamp(ficrest); |
| |
|
| |
epj=vector(1,nlstate+1); |
| for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| oldm=oldms;savm=savms; /* ZZ Segmentation fault */ |
oldm=oldms;savm=savms; /* ZZ Segmentation fault */ |
| cptcod= 0; /* To be deleted */ |
cptcod= 0; /* To be deleted */ |
|
Line 11931 Please run with mle=-1 to get a correct
|
Line 12507 Please run with mle=-1 to get a correct
|
| for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
| fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
| /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ |
/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ |
| epj=vector(1,nlstate+1); |
|
| printf("Computing age specific period (stable) prevalences in each health state \n"); |
printf("Computing age specific period (stable) prevalences in each health state \n"); |
| fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); |
fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); |
| for(age=bage; age <=fage ;age++){ |
for(age=bage; age <=fage ;age++){ |
|
Line 11969 Please run with mle=-1 to get a correct
|
Line 12544 Please run with mle=-1 to get a correct
|
| fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
| } |
} |
| } /* End vpopbased */ |
} /* End vpopbased */ |
| |
free_vector(epj,1,nlstate+1); |
| free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| free_vector(epj,1,nlstate+1); |
|
| printf("done selection\n");fflush(stdout); |
printf("done selection\n");fflush(stdout); |
| fprintf(ficlog,"done selection\n");fflush(ficlog); |
fprintf(ficlog,"done selection\n");fflush(ficlog); |
| |
|
| /*}*/ |
|
| } /* End k selection */ |
} /* End k selection */ |
| |
|
| printf("done State-specific expectancies\n");fflush(stdout); |
printf("done State-specific expectancies\n");fflush(stdout); |
| fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); |
fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); |
| |
|
| /*------- Variance of period (stable) prevalence------*/ |
/* variance-covariance of period prevalence*/ |
| |
varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); |
| strcpy(fileresvpl,"VPL_"); |
|
| strcat(fileresvpl,fileresu); |
|
| if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { |
|
| printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); |
|
| exit(0); |
|
| } |
|
| printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); |
|
| fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);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(ficresvpl,"\n#****** "); |
|
| printf("\n#****** "); |
|
| fprintf(ficlog,"\n#****** "); |
|
| for(j=1;j<=cptcoveff;j++) { |
|
| fprintf(ficresvpl,"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(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); |
|
| fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); |
|
| } |
|
| fprintf(ficresvpl,"******\n"); |
|
| printf("******\n"); |
|
| fprintf(ficlog,"******\n"); |
|
| |
|
| varpl=matrix(1,nlstate,(int) bage, (int) fage); |
|
| oldm=oldms;savm=savms; |
|
| varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart, nres); |
|
| free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
|
| /*}*/ |
|
| } |
|
| |
|
| 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_vector(weight,1,n); |
| free_imatrix(Tvard,1,NCOVMAX,1,2); |
free_imatrix(Tvard,1,NCOVMAX,1,2); |
|
Line 12046 Please run with mle=-1 to get a correct
|
Line 12575 Please run with mle=-1 to get a correct
|
| |
|
| /*---------- End : free ----------------*/ |
/*---------- End : free ----------------*/ |
| if (mobilav!=0 ||mobilavproj !=0) |
if (mobilav!=0 ||mobilavproj !=0) |
| free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */ |
free_ma3x(mobaverages,AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */ |
| free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
free_ma3x(probs,AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
| free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
| free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
| } /* mle==-3 arrives here for freeing */ |
} /* mle==-3 arrives here for freeing */ |
|
Line 12055 Please run with mle=-1 to get a correct
|
Line 12584 Please run with mle=-1 to get a correct
|
| free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(newms, 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_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); |
if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n); |
| free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n); |
if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); |
| free_matrix(coqvar,1,maxwav,1,n); |
if(nqv>=1)free_matrix(coqvar,1,nqv,1,n); |
| free_matrix(covar,0,NCOVMAX,1,n); |
free_matrix(covar,0,NCOVMAX,1,n); |
| free_matrix(matcov,1,npar,1,npar); |
free_matrix(matcov,1,npar,1,npar); |
| free_matrix(hess,1,npar,1,npar); |
free_matrix(hess,1,npar,1,npar); |
| /*free_vector(delti,1,npar);*/ |
/*free_vector(delti,1,npar);*/ |
| free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
| free_matrix(agev,1,maxwav,1,imx); |
free_matrix(agev,1,maxwav,1,imx); |
| |
free_ma3x(paramstart,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
| free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
| |
|
| free_ivector(ncodemax,1,NCOVMAX); |
free_ivector(ncodemax,1,NCOVMAX); |
|
Line 12140 Please run with mle=-1 to get a correct
|
Line 12670 Please run with mle=-1 to get a correct
|
| fclose(ficlog); |
fclose(ficlog); |
| /*------ End -----------*/ |
/*------ End -----------*/ |
| |
|
| |
|
| |
/* Executes gnuplot */ |
| |
|
| printf("Before Current directory %s!\n",pathcd); |
printf("Before Current directory %s!\n",pathcd); |
| #ifdef WIN32 |
#ifdef WIN32 |
|
Line 12208 end:
|
Line 12740 end:
|
| printf("\nType q for exiting: "); fflush(stdout); |
printf("\nType q for exiting: "); fflush(stdout); |
| scanf("%s",z); |
scanf("%s",z); |
| } |
} |
| |
exit(0); |
| } |
} |
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src