|
version 1.268, 2017/05/18 20:09:32
|
version 1.283, 2018/04/19 14:49:16
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.283 2018/04/19 14:49:16 brouard |
| |
Summary: Some minor bugs fixed |
| |
|
| |
Revision 1.282 2018/02/27 22:50:02 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.281 2018/02/27 19:25:23 brouard |
| |
Summary: Adding second argument for quitting |
| |
|
| |
Revision 1.280 2018/02/21 07:58:13 brouard |
| |
Summary: 0.99r15 |
| |
|
| |
New Makefile with recent VirtualBox 5.26. Bug in sqrt negatve in imach.c |
| |
|
| |
Revision 1.279 2017/07/20 13:35:01 brouard |
| |
Summary: temporary working |
| |
|
| |
Revision 1.278 2017/07/19 14:09:02 brouard |
| |
Summary: Bug for mobil_average=0 and prevforecast fixed(?) |
| |
|
| |
Revision 1.277 2017/07/17 08:53:49 brouard |
| |
Summary: BOM files can be read now |
| |
|
| |
Revision 1.276 2017/06/30 15:48:31 brouard |
| |
Summary: Graphs improvements |
| |
|
| |
Revision 1.275 2017/06/30 13:39:33 brouard |
| |
Summary: Saito's color |
| |
|
| |
Revision 1.274 2017/06/29 09:47:08 brouard |
| |
Summary: Version 0.99r14 |
| |
|
| |
Revision 1.273 2017/06/27 11:06:02 brouard |
| |
Summary: More documentation on projections |
| |
|
| |
Revision 1.272 2017/06/27 10:22:40 brouard |
| |
Summary: Color of backprojection changed from 6 to 5(yellow) |
| |
|
| |
Revision 1.271 2017/06/27 10:17:50 brouard |
| |
Summary: Some bug with rint |
| |
|
| |
Revision 1.270 2017/05/24 05:45:29 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.269 2017/05/23 08:39:25 brouard |
| |
Summary: Code into subroutine, cleanings |
| |
|
| Revision 1.268 2017/05/18 20:09:32 brouard |
Revision 1.268 2017/05/18 20:09:32 brouard |
| Summary: backprojection and confidence intervals of backprevalence |
Summary: backprojection and confidence intervals of backprevalence |
| |
|
|
Line 1011 typedef struct {
|
Line 1058 typedef struct {
|
| /* $State$ */ |
/* $State$ */ |
| #include "version.h" |
#include "version.h" |
| char version[]=__IMACH_VERSION__; |
char version[]=__IMACH_VERSION__; |
| char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; |
char copyright[]="April 2018,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; |
| char fullversion[]="$Revision$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
| char strstart[80]; |
char strstart[80]; |
| char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
|
Line 1084 FILE *ficrescveij;
|
Line 1131 FILE *ficrescveij;
|
| char filerescve[FILENAMELENGTH]; |
char filerescve[FILENAMELENGTH]; |
| FILE *ficresvij; |
FILE *ficresvij; |
| char fileresv[FILENAMELENGTH]; |
char fileresv[FILENAMELENGTH]; |
| FILE *ficresvpl; |
|
| char fileresvpl[FILENAMELENGTH]; |
|
| FILE *ficresvbl; |
|
| char fileresvbl[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 2489 void powell(double p[], double **xi, int
|
Line 2533 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 2932 double **pmij(double **ps, double *cov,
|
Line 2979 double **pmij(double **ps, double *cov,
|
| /* Diag(w_x) */ |
/* Diag(w_x) */ |
| /* Problem with prevacurrent which can be zero */ |
/* Problem with prevacurrent which can be zero */ |
| sumnew=0.; |
sumnew=0.; |
| /* for (ii=1;ii<=nlstate+ndeath;ii++){ */ |
/*for (ii=1;ii<=nlstate+ndeath;ii++){*/ |
| for (ii=1;ii<=nlstate;ii++){ /* Only on live states */ |
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]; |
sumnew+=prevacurrent[(int)agefin][ii][ij]; |
| } |
} |
| if(sumnew >0.01){ /* At least some value in the prevalence */ |
if(sumnew >0.01){ /* At least some value in the prevalence */ |
| for (ii=1;ii<=nlstate+ndeath;ii++){ |
for (ii=1;ii<=nlstate+ndeath;ii++){ |
| for (j=1;j<=nlstate+ndeath;j++) |
for (j=1;j<=nlstate+ndeath;j++) |
| doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0); |
doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0); |
| } |
} |
| }else{ |
}else{ |
| for (ii=1;ii<=nlstate+ndeath;ii++){ |
for (ii=1;ii<=nlstate+ndeath;ii++){ |
|
Line 3239 double ***hbxij(double ***po, int nhstep
|
Line 3287 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 3841 void likelione(FILE *ficres,double p[],
|
Line 3889 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 4323 void freqsummary(char fileres[], double
|
Line 4371 void freqsummary(char fileres[], double
|
| double ***freq; /* Frequencies */ |
double ***freq; /* Frequencies */ |
| double *x, *y, a=0.,b=0.,r=1., sa=0., sb=0.; /* for regression, y=b+m*x and r is the correlation coefficient */ |
double *x, *y, a=0.,b=0.,r=1., sa=0., sb=0.; /* for regression, y=b+m*x and r is the correlation coefficient */ |
| int no=0, linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb); |
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, *idq; |
| double **meanqt; |
double **meanqt; |
| double *pp, **prop, *posprop, *pospropt; |
double *pp, **prop, *posprop, *pospropt; |
| double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; |
double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; |
|
Line 4336 void freqsummary(char fileres[], double
|
Line 4384 void freqsummary(char fileres[], double
|
| pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ |
pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ |
| /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ |
/* prop=matrix(1,nlstate,iagemin,iagemax+3); */ |
| meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ |
meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ |
| |
idq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ |
| meanqt=matrix(1,lastpass,1,nqtveff); |
meanqt=matrix(1,lastpass,1,nqtveff); |
| strcpy(fileresp,"P_"); |
strcpy(fileresp,"P_"); |
| strcat(fileresp,fileresu); |
strcat(fileresp,fileresu); |
|
Line 4444 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4493 Title=%s <br>Datafile=%s Firstpass=%d La
|
| posprop[i]=0; |
posprop[i]=0; |
| pospropt[i]=0; |
pospropt[i]=0; |
| } |
} |
| /* for (z1=1; z1<= nqfveff; z1++) { */ |
for (z1=1; z1<= nqfveff; z1++) { /* zeroing for each combination j1 as well as for the total */ |
| /* meanq[z1]+=0.; */ |
idq[z1]+=0.; |
| |
meanq[z1]+=0.; |
| |
} |
| |
/* for (z1=1; z1<= nqtveff; z1++) { */ |
| /* for(m=1;m<=lastpass;m++){ */ |
/* for(m=1;m<=lastpass;m++){ */ |
| /* meanqt[m][z1]=0.; */ |
/* meanqt[m][z1]=0.; */ |
| /* } */ |
/* } */ |
| /* } */ |
/* } */ |
| |
|
| /* dateintsum=0; */ |
/* dateintsum=0; */ |
| /* k2cpt=0; */ |
/* k2cpt=0; */ |
| |
|
|
Line 4460 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4511 Title=%s <br>Datafile=%s Firstpass=%d La
|
| if(j !=0){ |
if(j !=0){ |
| if(anyvaryingduminmodel==0){ /* If All fixed covariates */ |
if(anyvaryingduminmodel==0){ /* If All fixed covariates */ |
| if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
| /* for (z1=1; z1<= nqfveff; z1++) { */ |
/* for (z1=1; z1<= nqfveff; z1++) { */ |
| /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ |
/* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ |
| /* } */ |
/* } */ |
| for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ |
for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ |
|
Line 4526 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4577 Title=%s <br>Datafile=%s Firstpass=%d La
|
| /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ |
/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ |
| freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */ |
freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */ |
| } |
} |
| |
for (z1=1; z1<= nqfveff; z1++) { |
| |
idq[z1]++; |
| |
meanq[z1]+=covar[ncovcol+z1][iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */ |
| |
/* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ |
| |
} |
| } /* end if between passes */ |
} /* end if between passes */ |
| if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99) && (j==0)) { |
if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99) && (j==0)) { |
| dateintsum=dateintsum+k2; /* on all covariates ?*/ |
dateintsum=dateintsum+k2; /* on all covariates ?*/ |
|
Line 4573 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4629 Title=%s <br>Datafile=%s Firstpass=%d La
|
| fprintf(ficresphtmfr, "**********</h3>\n"); |
fprintf(ficresphtmfr, "**********</h3>\n"); |
| fprintf(ficlog, "**********\n"); |
fprintf(ficlog, "**********\n"); |
| } |
} |
| |
/* |
| |
Printing means of quantitative variables if any |
| |
*/ |
| |
for (z1=1; z1<= nqfveff; z1++) { |
| |
fprintf(ficresphtmfr,"V quantitative id %d, number of idividuals= %f, sum=%f", z1, idq[z1], meanq[z1]); |
| |
fprintf(ficresphtmfr,", mean=%f<p>\n",meanq[z1]/idq[z1]); |
| |
} |
| |
/* for (z1=1; z1<= nqtveff; z1++) { */ |
| |
/* for(m=1;m<=lastpass;m++){ */ |
| |
/* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f<p>\n", z1, m, meanqt[m][z1]); */ |
| |
/* } */ |
| |
/* } */ |
| |
|
| fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
| if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ |
if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ |
| fprintf(ficresp, " Age"); |
fprintf(ficresp, " Age"); |
|
Line 4853 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4922 Title=%s <br>Datafile=%s Firstpass=%d La
|
| fclose(ficresp); |
fclose(ficresp); |
| fclose(ficresphtm); |
fclose(ficresphtm); |
| fclose(ficresphtmfr); |
fclose(ficresphtmfr); |
| |
free_vector(idq,1,nqfveff); |
| free_vector(meanq,1,nqfveff); |
free_vector(meanq,1,nqfveff); |
| free_matrix(meanqt,1,lastpass,1,nqtveff); |
free_matrix(meanqt,1,lastpass,1,nqtveff); |
| free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); |
free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); |
|
Line 5463 void concatwav(int wav[], int **dh, int
|
Line 5533 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 5694 void concatwav(int wav[], int **dh, int
|
Line 5764 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 5709 void concatwav(int wav[], int **dh, int
|
Line 5780 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 5742 void concatwav(int wav[], int **dh, int
|
Line 5814 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 5753 void concatwav(int wav[], int **dh, int
|
Line 5826 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 5818 void concatwav(int wav[], int **dh, int
|
Line 5891 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 5873 void concatwav(int wav[], int **dh, int
|
Line 5946 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 5885 void concatwav(int wav[], int **dh, int
|
Line 5961 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 5934 void concatwav(int wav[], int **dh, int
|
Line 6015 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 5957 void concatwav(int wav[], int **dh, int
|
Line 6042 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 5974 void concatwav(int wav[], int **dh, int
|
Line 6065 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 6062 void concatwav(int wav[], int **dh, int
|
Line 6157 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);*/ |
|
Line 6187 void concatwav(int wav[], int **dh, int
|
Line 6282 void concatwav(int wav[], int **dh, int
|
| |
|
| |
|
| /************ Variance of backprevalence limit ******************/ |
/************ Variance of backprevalence limit ******************/ |
| void varbrevlim(char fileres[], double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres) |
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) |
| { |
{ |
| /* Variance of backward prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ |
/* Variance of backward prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
|
Line 6583 To be simple, these graphs help to under
|
Line 6678 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 6614 To be simple, these graphs help to under
|
Line 6714 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 6652 void printinghtml(char fileresu[], char
|
Line 6752 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 6790 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 6890 divided by h: <sub>h</sub>P<sub>ij</sub>
|
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\ |
fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\ |
| Or probability to survive in various states (1 to %d) being in state %d at different ages. \ |
Or probability to survive in various states (1 to %d) being in state %d at different ages. \ |
| <a href=\"%s_%d-%d-%d.svg\">%s_%d%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); |
<a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); |
| } |
} |
| /* Period (stable) prevalence in each health state */ |
/* Period (stable) prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
|
Line 6807 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 6907 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 an observed weighted state (from 1 to %d). <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){ |
if(backcast==1){ |
| /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */ |
/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
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) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \ |
| <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); |
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 6921 true period expectancies (those weighted
|
Line 7024 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, int offbyear){ |
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 6930 void printinggnuplot(char fileresu[], ch
|
Line 7033 void printinggnuplot(char fileresu[], ch
|
| 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 6943 void printinggnuplot(char fileresu[], ch
|
Line 7048 void printinggnuplot(char fileresu[], ch
|
| /*#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 7012 void printinggnuplot(char fileresu[], ch
|
Line 7131 void printinggnuplot(char fileresu[], ch
|
| |
|
| 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 7034 void printinggnuplot(char fileresu[], ch
|
Line 7154 void printinggnuplot(char fileresu[], ch
|
| |
|
| 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 7092 void printinggnuplot(char fileresu[], ch
|
Line 7212 void printinggnuplot(char fileresu[], ch
|
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); |
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 ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); |
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 ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"\" w l lt 1"); |
fprintf(ficgp,"\" t\"\" w l lt 4"); |
| } /* end if backprojcast */ |
} /* 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 7503 set ter svg size 640, 480\nunset log y\n
|
Line 7624 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 ); |
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 7542 set ter svg size 640, 480\nunset log y\n
|
Line 7661 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 ); |
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 7618 set ter svg size 640, 480\nunset log y\n
|
Line 7737 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 'bw%d' with line lc variable ", \ |
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 ); |
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/0):5 with labels center not ", \ |
fprintf(ficgp," (($1-$2) == %d ) ? $%d : 1/0):1 with labels center not ", \ |
| offbyear, \ |
offbyear, \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1) ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1) ); |
| }else |
}else |
| fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ", \ |
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 ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i ); |
| }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 7657 set ter svg size 640, 480\nunset log y\n
|
Line 7774 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/0):5 t 'bw%d' with line lc variable", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d : 1/0):%d t 'bw%d' with line lc variable", gplotcondition, \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1),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 && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */ |
| fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", gplotcondition, \ |
fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d : 1/0):%d with labels center not ", gplotcondition, \ |
| offbyear, \ |
iyearc,iagec,offbyear, \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1) ); |
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*/ |
/* '' 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 7723 set ter svg size 640, 480\nunset log y\n
|
Line 7840 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 7740 set ter svg size 640, 480\nunset log y\n
|
Line 7857 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 7860 set ter svg size 640, 480\nunset log y\n
|
Line 7979 set ter svg size 640, 480\nunset log y\n
|
| } |
} |
| 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 7874 set ter svg size 640, 480\nunset log y\n
|
Line 7993 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 7898 set ter svg size 640, 480\nunset log y\n
|
Line 8018 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 8112 set ter svg size 640, 480\nunset log y\n
|
Line 8232 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 8151 void prevforecast(char fileres[], double
|
Line 8271 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 8196 void prevforecast(char fileres[], double
|
Line 8321 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); |
|
Line 8217 void prevforecast(char fileres[], double
|
Line 8343 void prevforecast(char fileres[], double
|
| for(j=1; j<=nlstate+ndeath;j++) { |
for(j=1; j<=nlstate+ndeath;j++) { |
| ppij=0.; |
ppij=0.; |
| for(i=1; i<=nlstate;i++) { |
for(i=1; i<=nlstate;i++) { |
| /* if (mobilav>=1) */ |
if (mobilav>=1) |
| ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k]; |
ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k]; |
| /* else { */ /* even if mobilav==-1 we use mobaverage */ |
else { /* even if mobilav==-1 we use mobaverage, probs may not sums to 1 */ |
| /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ |
ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; |
| /* } */ |
} |
| fprintf(ficresf," %.3f", p3mat[i][j][h]); |
fprintf(ficresf," %.3f", p3mat[i][j][h]); |
| } /* end i */ |
} /* end i */ |
| fprintf(ficresf," %.3f", ppij); |
fprintf(ficresf," %.3f", ppij); |
|
Line 8239 void prevforecast(char fileres[], double
|
Line 8365 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 */ |
|
Line 8283 void prevbackforecast(char fileres[], do
|
Line 8410 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 8304 void prevbackforecast(char fileres[], do
|
Line 8436 void prevbackforecast(char fileres[], do
|
| |
|
| 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 8333 void prevbackforecast(char fileres[], do
|
Line 8461 void prevbackforecast(char fileres[], do
|
| /* 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); |
| printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); |
/* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */ |
| /* for (agec=fage; agec>=(ageminpar-1); agec--){ */ |
/* for (agec=bage; agec<=agemax-1; agec++){ /\* testing *\/ */ |
| /* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */ |
for (agec=bage; agec<=fage; agec++){ /* testing */ |
| for (agec=bage; agec<=agemax-1; agec++){ /* testing */ |
|
| /* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/ |
/* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/ |
| nhstepm=(int) rint((agec-agelim)*YEARM/stepm); |
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 */ |
/* 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); */ |
/* 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 */ |
/* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */ |
|
Line 8360 void prevbackforecast(char fileres[], do
|
Line 8488 void prevbackforecast(char fileres[], do
|
| ppij=0.;ppi=0.; |
ppij=0.;ppi=0.; |
| for(j=1; j<=nlstate;j++) { |
for(j=1; j<=nlstate;j++) { |
| /* if (mobilav==1) */ |
/* if (mobilav==1) */ |
| ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k]; |
ppij=ppij+p3mat[i][j][h]*prevacurrent[(int)agec][j][k]; |
| ppi=ppi+mobaverage[(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]; */ |
| /* } */ |
/* } */ |
|
Line 8386 void prevbackforecast(char fileres[], do
|
Line 8516 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 10442 int main(int argc, char *argv[])
|
Line 10689 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 10494 int main(int argc, char *argv[])
|
Line 10743 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 **varbpl; /* Variances of back 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 10576 int main(int argc, char *argv[])
|
Line 10826 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 10655 int main(int argc, char *argv[])
|
Line 10910 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 10665 int main(int argc, char *argv[])
|
Line 10918 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 10685 int main(int argc, char *argv[])
|
Line 10973 int main(int argc, char *argv[])
|
| title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){ |
title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){ |
| if (num_filled != 5) { |
if (num_filled != 5) { |
| printf("Should be 5 parameters\n"); |
printf("Should be 5 parameters\n"); |
| |
fprintf(ficlog,"Should be 5 parameters\n"); |
| } |
} |
| numlinepar++; |
numlinepar++; |
| printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); |
printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); |
| |
fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); |
| |
fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); |
| |
fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); |
| } |
} |
| /* Second parameter line */ |
/* Second parameter line */ |
| while(fgets(line, MAXLINE, ficpar)) { |
while(fgets(line, MAXLINE, ficpar)) { |
| /* If line starts with a # it is a comment */ |
/* while(fscanf(ficpar,"%[^\n]", line)) { */ |
| |
/* If line starts with a # it is a comment. Strangely fgets reads the EOL and fputs doesn't */ |
| if (line[0] == '#') { |
if (line[0] == '#') { |
| numlinepar++; |
numlinepar++; |
| fputs(line,stdout); |
printf("%s",line); |
| fputs(line,ficparo); |
fprintf(ficres,"%s",line); |
| fputs(line,ficlog); |
fprintf(ficparo,"%s",line); |
| |
fprintf(ficlog,"%s",line); |
| continue; |
continue; |
| }else |
}else |
| break; |
break; |
|
Line 10706 int main(int argc, char *argv[])
|
Line 11000 int main(int argc, char *argv[])
|
| if (num_filled != 11) { |
if (num_filled != 11) { |
| printf("Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); |
printf("Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); |
| printf("but line=%s\n",line); |
printf("but line=%s\n",line); |
| |
fprintf(ficlog,"Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); |
| |
fprintf(ficlog,"but line=%s\n",line); |
| } |
} |
| printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); |
printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); |
| |
fprintf(ficparo,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); |
| |
fprintf(ficres,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); |
| |
fprintf(ficlog,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); |
| } |
} |
| /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ |
/* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ |
| /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ |
/*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ |
|
Line 10716 int main(int argc, char *argv[])
|
Line 11015 int main(int argc, char *argv[])
|
| /* If line starts with a # it is a comment */ |
/* If line starts with a # it is a comment */ |
| if (line[0] == '#') { |
if (line[0] == '#') { |
| numlinepar++; |
numlinepar++; |
| fputs(line,stdout); |
printf("%s",line); |
| fputs(line,ficparo); |
fprintf(ficres,"%s",line); |
| fputs(line,ficlog); |
fprintf(ficparo,"%s",line); |
| |
fprintf(ficlog,"%s",line); |
| 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 10744 int main(int argc, char *argv[])
|
Line 11039 int main(int argc, char *argv[])
|
| } |
} |
| /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */ |
/* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */ |
| printf("model=1+age+%s\n",model);fflush(stdout); |
printf("model=1+age+%s\n",model);fflush(stdout); |
| |
fprintf(ficparo,"model=1+age+%s\n",model);fflush(stdout); |
| |
fprintf(ficres,"model=1+age+%s\n",model);fflush(stdout); |
| |
fprintf(ficlog,"model=1+age+%s\n",model);fflush(stdout); |
| } |
} |
| /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */ |
/* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */ |
| /* numlinepar=numlinepar+3; /\* In general *\/ */ |
/* numlinepar=numlinepar+3; /\* In general *\/ */ |
| /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ |
/* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ |
| fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); |
/* fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); */ |
| fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); |
/* fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); */ |
| 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 10979 Please run with mle=-1 to get a correct
|
Line 11277 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 11326 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 11614 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 */ |
|
Line 11601 Please run with mle=-1 to get a correct
|
Line 11910 Please run with mle=-1 to get a correct
|
| printf("\n"); |
printf("\n"); |
| |
|
| /*--------- results files --------------*/ |
/*--------- results files --------------*/ |
| fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, weightopt,model); |
/* fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, weightopt,model); */ |
| |
|
| |
|
| fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
|
Line 11889 Please run with mle=-1 to get a correct
|
Line 12198 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 11904 Please run with mle=-1 to get a correct
|
Line 12216 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 11954 Please run with mle=-1 to get a correct
|
Line 12268 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)agemin, (int)anback1-(int)agemax+1); |
/* 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 11994 Please run with mle=-1 to get a correct
|
Line 12309 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+ndeath;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.; |
|
Line 12015 Please run with mle=-1 to get a correct
|
Line 12331 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 12048 Please run with mle=-1 to get a correct
|
Line 12360 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); |
| |
|
| /*------- Variance of back (stable) prevalence------*/ |
|
| |
|
| strcpy(fileresvbl,"VBL_"); |
free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
| strcat(fileresvbl,fileresu); |
|
| if((ficresvbl=fopen(fileresvbl,"w"))==NULL) { |
|
| printf("Problem with variance of back (stable) prevalence resultfile: %s\n", fileresvbl); |
|
| exit(0); |
|
| } |
|
| printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout); |
|
| fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog); |
|
| |
|
| /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
|
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
|
| |
|
| i1=pow(2,cptcoveff); |
|
| if (cptcovn < 1){i1=1;} |
|
| |
|
| for(nres=1; nres <= nresult; nres++) /* For each resultline */ |
|
| for(k=1; k<=i1;k++){ |
|
| if(i1 != 1 && TKresult[nres]!= k) |
|
| continue; |
|
| fprintf(ficresvbl,"\n#****** "); |
|
| printf("\n#****** "); |
|
| fprintf(ficlog,"\n#****** "); |
|
| for(j=1;j<=cptcoveff;j++) { |
|
| fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
|
| fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
|
| printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
|
| } |
|
| for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ |
|
| printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); |
|
| fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); |
|
| fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); |
|
| } |
|
| fprintf(ficresvbl,"******\n"); |
|
| printf("******\n"); |
|
| fprintf(ficlog,"******\n"); |
|
| |
|
| varbpl=matrix(1,nlstate,(int) bage, (int) fage); |
|
| oldm=oldms;savm=savms; |
|
| |
|
| varbrevlim(fileres, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, &ncvyear, k, strstart, nres); |
|
| free_matrix(varbpl,1,nlstate,(int) bage, (int)fage); |
|
| /*}*/ |
|
| } |
|
| |
|
| fclose(ficresvbl); |
|
| printf("done variance-covariance of back prevalence\n");fflush(stdout); |
|
| fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog); |
|
| |
|
| 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 */ |
| free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
|
| |
|
| |
|
| /* ------ Other prevalence ratios------------ */ |
/* ------ Other prevalence ratios------------ */ |
|
Line 12165 Please run with mle=-1 to get a correct
|
Line 12429 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 12177 Please run with mle=-1 to get a correct
|
Line 12441 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 12206 Please run with mle=-1 to get a correct
|
Line 12468 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 12270 Please run with mle=-1 to get a correct
|
Line 12529 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 12286 Please run with mle=-1 to get a correct
|
Line 12545 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 12324 Please run with mle=-1 to get a correct
|
Line 12582 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); |
|
Line 12402 Please run with mle=-1 to get a correct
|
Line 12613 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 12420 Please run with mle=-1 to get a correct
|
Line 12631 Please run with mle=-1 to get a correct
|
| /*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 12496 Please run with mle=-1 to get a correct
|
Line 12708 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 12564 end:
|
Line 12778 end:
|
| printf("\nType q for exiting: "); fflush(stdout); |
printf("\nType q for exiting: "); fflush(stdout); |
| scanf("%s",z); |
scanf("%s",z); |
| } |
} |
| |
printf("End\n"); |
| |
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