|
version 1.269, 2017/05/23 08:39:25
|
version 1.284, 2018/04/20 05:22:13
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.284 2018/04/20 05:22:13 brouard |
| |
Summary: Computing mean and stdeviation of fixed quantitative variables |
| |
|
| |
Revision 1.283 2018/04/19 14:49:16 brouard |
| |
Summary: Some minor bugs fixed |
| |
|
| |
Revision 1.282 2018/02/27 22:50:02 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.281 2018/02/27 19:25:23 brouard |
| |
Summary: Adding second argument for quitting |
| |
|
| |
Revision 1.280 2018/02/21 07:58:13 brouard |
| |
Summary: 0.99r15 |
| |
|
| |
New Makefile with recent VirtualBox 5.26. Bug in sqrt negatve in imach.c |
| |
|
| |
Revision 1.279 2017/07/20 13:35:01 brouard |
| |
Summary: temporary working |
| |
|
| |
Revision 1.278 2017/07/19 14:09:02 brouard |
| |
Summary: Bug for mobil_average=0 and prevforecast fixed(?) |
| |
|
| |
Revision 1.277 2017/07/17 08:53:49 brouard |
| |
Summary: BOM files can be read now |
| |
|
| |
Revision 1.276 2017/06/30 15:48:31 brouard |
| |
Summary: Graphs improvements |
| |
|
| |
Revision 1.275 2017/06/30 13:39:33 brouard |
| |
Summary: Saito's color |
| |
|
| |
Revision 1.274 2017/06/29 09:47:08 brouard |
| |
Summary: Version 0.99r14 |
| |
|
| |
Revision 1.273 2017/06/27 11:06:02 brouard |
| |
Summary: More documentation on projections |
| |
|
| |
Revision 1.272 2017/06/27 10:22:40 brouard |
| |
Summary: Color of backprojection changed from 6 to 5(yellow) |
| |
|
| |
Revision 1.271 2017/06/27 10:17:50 brouard |
| |
Summary: Some bug with rint |
| |
|
| |
Revision 1.270 2017/05/24 05:45:29 brouard |
| |
*** empty log message *** |
| |
|
| Revision 1.269 2017/05/23 08:39:25 brouard |
Revision 1.269 2017/05/23 08:39:25 brouard |
| Summary: Code into subroutine, cleanings |
Summary: Code into subroutine, cleanings |
| |
|
|
Line 1014 typedef struct {
|
Line 1061 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 2489 void powell(double p[], double **xi, int
|
Line 2536 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 3240 double ***hbxij(double ***po, int nhstep
|
Line 3290 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 3842 void likelione(FILE *ficres,double p[],
|
Line 3892 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 4324 void freqsummary(char fileres[], double
|
Line 4374 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, *stdq, *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 4337 void freqsummary(char fileres[], double
|
Line 4387 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 */ |
| |
stdq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ |
| |
idq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ |
| meanqt=matrix(1,lastpass,1,nqtveff); |
meanqt=matrix(1,lastpass,1,nqtveff); |
| strcpy(fileresp,"P_"); |
strcpy(fileresp,"P_"); |
| strcat(fileresp,fileresu); |
strcat(fileresp,fileresu); |
|
Line 4445 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4497 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.; |
| |
stdq[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 4461 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4516 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++) { */ |
|
| /* 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 */ |
| /* if(Tvaraff[z1] ==-20){ */ |
/* if(Tvaraff[z1] ==-20){ */ |
| /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
/* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
|
Line 4484 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4536 Title=%s <br>Datafile=%s Firstpass=%d La
|
| }/* end j==0 */ |
}/* end j==0 */ |
| if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ |
if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ |
| /* for(m=firstpass; m<=lastpass; m++){ */ |
/* for(m=firstpass; m<=lastpass; m++){ */ |
| for(mi=1; mi<wav[iind];mi++){ /* For that wave */ |
for(mi=1; mi<wav[iind];mi++){ /* For each wave */ |
| m=mw[mi][iind]; |
m=mw[mi][iind]; |
| if(j!=0){ |
if(j!=0){ |
| if(anyvaryingduminmodel==1){ /* Some are varying covariates */ |
if(anyvaryingduminmodel==1){ /* Some are varying covariates */ |
|
Line 4504 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4556 Title=%s <br>Datafile=%s Firstpass=%d La
|
| }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ |
}/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ |
| } /* end j==0 */ |
} /* end j==0 */ |
| /* bool =0 we keep that guy which corresponds to the combination of dummy values */ |
/* bool =0 we keep that guy which corresponds to the combination of dummy values */ |
| if(bool==1){ |
if(bool==1){ /*Selected */ |
| /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] |
/* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] |
| and mw[mi+1][iind]. dh depends on stepm. */ |
and mw[mi+1][iind]. dh depends on stepm. */ |
| agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ |
agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ |
|
Line 4522 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4574 Title=%s <br>Datafile=%s Firstpass=%d La
|
| if(s[m][iind]==-1) |
if(s[m][iind]==-1) |
| printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.)); |
printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.)); |
| freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ |
freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ |
| |
for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */ |
| |
idq[z1]=idq[z1]+weight[iind]; |
| |
meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /* Computes mean of quantitative with selected filter */ |
| |
stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */ |
| |
} |
| /* if((int)agev[m][iind] == 55) */ |
/* if((int)agev[m][iind] == 55) */ |
| /* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ |
/* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ |
| /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ |
/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ |
|
Line 4537 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4594 Title=%s <br>Datafile=%s Firstpass=%d La
|
| bool=1; |
bool=1; |
| }/* end bool 2 */ |
}/* end bool 2 */ |
| } /* end m */ |
} /* end m */ |
| |
/* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */ |
| |
/* idq[z1]=idq[z1]+weight[iind]; */ |
| |
/* meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /\* Computes mean of quantitative with selected filter *\/ */ |
| |
/* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/ /\* Computes mean of quantitative with selected filter *\/ */ |
| |
/* } */ |
| } /* end bool */ |
} /* end bool */ |
| } /* end iind = 1 to imx */ |
} /* end iind = 1 to imx */ |
| /* prop[s][age] is feeded for any initial and valid live state as well as |
/* prop[s][age] is feeded for any initial and valid live state as well as |
|
Line 4574 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4636 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(ficlog,"Mean of quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]); |
| |
fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]); |
| |
if(weightopt==1){ |
| |
printf(" Weighted mean and standard deviation of"); |
| |
fprintf(ficlog," Weighted mean and standard deviation of"); |
| |
fprintf(ficresphtmfr," Weighted mean and standard deviation of"); |
| |
} |
| |
printf(" quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
| |
fprintf(ficlog," quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
| |
fprintf(ficresphtmfr," quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)<p>\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
| |
} |
| |
/* for (z1=1; z1<= nqtveff; z1++) { */ |
| |
/* for(m=1;m<=lastpass;m++){ */ |
| |
/* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f<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 4808 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4891 Title=%s <br>Datafile=%s Firstpass=%d La
|
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| } |
} |
| } |
} |
| } |
} /* end of state i */ |
| printf("#Freqsummary\n"); |
printf("#Freqsummary\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| for(s1=-1; s1 <=nlstate+ndeath; s1++){ |
for(s1=-1; s1 <=nlstate+ndeath; s1++){ |
|
Line 4854 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4937 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_vector(stdq,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); |
| free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); |
free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); |
|
Line 5464 void concatwav(int wav[], int **dh, int
|
Line 5549 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 5745 void concatwav(int wav[], int **dh, int
|
Line 5830 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 5756 void concatwav(int wav[], int **dh, int
|
Line 5842 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 5821 void concatwav(int wav[], int **dh, int
|
Line 5907 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 5876 void concatwav(int wav[], int **dh, int
|
Line 5962 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 5888 void concatwav(int wav[], int **dh, int
|
Line 5977 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 5937 void concatwav(int wav[], int **dh, int
|
Line 6031 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 5960 void concatwav(int wav[], int **dh, int
|
Line 6058 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 5977 void concatwav(int wav[], int **dh, int
|
Line 6081 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 6586 To be simple, these graphs help to under
|
Line 6694 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 6617 To be simple, these graphs help to under
|
Line 6730 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 6655 void printinghtml(char fileresu[], char
|
Line 6768 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 6793 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 6906 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 6810 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 6923 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 6924 true period expectancies (those weighted
|
Line 7040 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 6933 void printinggnuplot(char fileresu[], ch
|
Line 7049 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 6946 void printinggnuplot(char fileresu[], ch
|
Line 7064 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 7015 void printinggnuplot(char fileresu[], ch
|
Line 7147 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 7037 void printinggnuplot(char fileresu[], ch
|
Line 7170 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 7095 void printinggnuplot(char fileresu[], ch
|
Line 7228 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 7506 set ter svg size 640, 480\nunset log y\n
|
Line 7640 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 7545 set ter svg size 640, 480\nunset log y\n
|
Line 7677 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 7621 set ter svg size 640, 480\nunset log y\n
|
Line 7753 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 7660 set ter svg size 640, 480\nunset log y\n
|
Line 7790 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 7726 set ter svg size 640, 480\nunset log y\n
|
Line 7856 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 7743 set ter svg size 640, 480\nunset log y\n
|
Line 7873 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 7863 set ter svg size 640, 480\nunset log y\n
|
Line 7995 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 7877 set ter svg size 640, 480\nunset log y\n
|
Line 8009 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 7901 set ter svg size 640, 480\nunset log y\n
|
Line 8034 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 8154 set ter svg size 640, 480\nunset log y\n
|
Line 8287 set ter svg size 640, 480\nunset log y\n
|
| 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 8199 set ter svg size 640, 480\nunset log y\n
|
Line 8337 set ter svg size 640, 480\nunset log y\n
|
| 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 8220 set ter svg size 640, 480\nunset log y\n
|
Line 8359 set ter svg size 640, 480\nunset log y\n
|
| 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]*prev[(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 8287 set ter svg size 640, 480\nunset log y\n
|
Line 8426 set ter svg size 640, 480\nunset log y\n
|
| 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 8334 set ter svg size 640, 480\nunset log y\n
|
Line 8478 set ter svg size 640, 480\nunset log y\n
|
| 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=bage; agec<=agemax-1; agec++){ /* testing */ |
/* for (agec=bage; agec<=agemax-1; agec++){ /\* testing *\/ */ |
| |
for (agec=bage; agec<=fage; 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 10559 int main(int argc, char *argv[])
|
Line 10705 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 10614 int main(int argc, char *argv[])
|
Line 10762 int main(int argc, char *argv[])
|
| |
|
| 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 10693 int main(int argc, char *argv[])
|
Line 10842 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 10772 int main(int argc, char *argv[])
|
Line 10926 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 10782 int main(int argc, char *argv[])
|
Line 10934 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 10802 int main(int argc, char *argv[])
|
Line 10989 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 10823 int main(int argc, char *argv[])
|
Line 11016 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 10833 int main(int argc, char *argv[])
|
Line 11031 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 10861 int main(int argc, char *argv[])
|
Line 11055 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 11096 Please run with mle=-1 to get a correct
|
Line 11293 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 11443 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 11630 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 11718 Please run with mle=-1 to get a correct
|
Line 11926 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 12006 Please run with mle=-1 to get a correct
|
Line 12214 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 12021 Please run with mle=-1 to get a correct
|
Line 12232 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 12071 Please run with mle=-1 to get a correct
|
Line 12284 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 12510 Please run with mle=-1 to get a correct
|
Line 12724 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 12578 end:
|
Line 12794 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