|
version 1.277, 2017/07/17 08:53:49
|
version 1.303, 2021/02/11 19:50:15
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.303 2021/02/11 19:50:15 brouard |
| |
* (Module): imach.c Someone entered 'results:' instead of 'result:'. Now it is an error which is printed. |
| |
|
| |
Revision 1.302 2020/02/22 21:00:05 brouard |
| |
* (Module): imach.c Update mle=-3 (for computing Life expectancy |
| |
and life table from the data without any state) |
| |
|
| |
Revision 1.301 2019/06/04 13:51:20 brouard |
| |
Summary: Error in 'r'parameter file backcast yearsbproj instead of yearsfproj |
| |
|
| |
Revision 1.300 2019/05/22 19:09:45 brouard |
| |
Summary: version 0.99r19 of May 2019 |
| |
|
| |
Revision 1.299 2019/05/22 18:37:08 brouard |
| |
Summary: Cleaned 0.99r19 |
| |
|
| |
Revision 1.298 2019/05/22 18:19:56 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.297 2019/05/22 17:56:10 brouard |
| |
Summary: Fix bug by moving date2dmy and nhstepm which gaefin=-1 |
| |
|
| |
Revision 1.296 2019/05/20 13:03:18 brouard |
| |
Summary: Projection syntax simplified |
| |
|
| |
|
| |
We can now start projections, forward or backward, from the mean date |
| |
of inteviews up to or down to a number of years of projection: |
| |
prevforecast=1 yearsfproj=15.3 mobil_average=0 |
| |
or |
| |
prevforecast=1 starting-proj-date=1/1/2007 final-proj-date=12/31/2017 mobil_average=0 |
| |
or |
| |
prevbackcast=1 yearsbproj=12.3 mobil_average=1 |
| |
or |
| |
prevbackcast=1 starting-back-date=1/10/1999 final-back-date=1/1/1985 mobil_average=1 |
| |
|
| |
Revision 1.295 2019/05/18 09:52:50 brouard |
| |
Summary: doxygen tex bug |
| |
|
| |
Revision 1.294 2019/05/16 14:54:33 brouard |
| |
Summary: There was some wrong lines added |
| |
|
| |
Revision 1.293 2019/05/09 15:17:34 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.292 2019/05/09 14:17:20 brouard |
| |
Summary: Some updates |
| |
|
| |
Revision 1.291 2019/05/09 13:44:18 brouard |
| |
Summary: Before ncovmax |
| |
|
| |
Revision 1.290 2019/05/09 13:39:37 brouard |
| |
Summary: 0.99r18 unlimited number of individuals |
| |
|
| |
The number n which was limited to 20,000 cases is now unlimited, from firstobs to lastobs. If the number is too for the virtual memory, probably an error will occur. |
| |
|
| |
Revision 1.289 2018/12/13 09:16:26 brouard |
| |
Summary: Bug for young ages (<-30) will be in r17 |
| |
|
| |
Revision 1.288 2018/05/02 20:58:27 brouard |
| |
Summary: Some bugs fixed |
| |
|
| |
Revision 1.287 2018/05/01 17:57:25 brouard |
| |
Summary: Bug fixed by providing frequencies only for non missing covariates |
| |
|
| |
Revision 1.286 2018/04/27 14:27:04 brouard |
| |
Summary: some minor bugs |
| |
|
| |
Revision 1.285 2018/04/21 21:02:16 brouard |
| |
Summary: Some bugs fixed, valgrind tested |
| |
|
| |
Revision 1.284 2018/04/20 05:22:13 brouard |
| |
Summary: Computing mean and stdeviation of fixed quantitative variables |
| |
|
| |
Revision 1.283 2018/04/19 14:49:16 brouard |
| |
Summary: Some minor bugs fixed |
| |
|
| |
Revision 1.282 2018/02/27 22:50:02 brouard |
| |
*** empty log message *** |
| |
|
| |
Revision 1.281 2018/02/27 19:25:23 brouard |
| |
Summary: Adding second argument for quitting |
| |
|
| |
Revision 1.280 2018/02/21 07:58:13 brouard |
| |
Summary: 0.99r15 |
| |
|
| |
New Makefile with recent VirtualBox 5.26. Bug in sqrt negatve in imach.c |
| |
|
| |
Revision 1.279 2017/07/20 13:35:01 brouard |
| |
Summary: temporary working |
| |
|
| |
Revision 1.278 2017/07/19 14:09:02 brouard |
| |
Summary: Bug for mobil_average=0 and prevforecast fixed(?) |
| |
|
| Revision 1.277 2017/07/17 08:53:49 brouard |
Revision 1.277 2017/07/17 08:53:49 brouard |
| Summary: BOM files can be read now |
Summary: BOM files can be read now |
| |
|
|
Line 1011 typedef struct {
|
Line 1105 typedef struct {
|
| #define NINTERVMAX 8 |
#define NINTERVMAX 8 |
| #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ |
#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ |
| #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
| #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
| #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| /*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/ |
/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/ |
| #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 |
#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 |
| #define MAXN 20000 |
/*#define MAXN 20000 */ /* Should by replaced by nobs, real number of observations and unlimited */ |
| #define YEARM 12. /**< Number of months per year */ |
#define YEARM 12. /**< Number of months per year */ |
| /* #define AGESUP 130 */ |
/* #define AGESUP 130 */ |
| #define AGESUP 150 |
/* #define AGESUP 150 */ |
| |
#define AGESUP 200 |
| #define AGEINF 0 |
#define AGEINF 0 |
| #define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */ |
#define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */ |
| #define AGEBASE 40 |
#define AGEBASE 40 |
|
Line 1038 typedef struct {
|
Line 1133 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[]="May 2019,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020"; |
| 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 1062 int nqfveff=0; /**< nqfveff Number of Qu
|
Line 1157 int nqfveff=0; /**< nqfveff Number of Qu
|
| int ntveff=0; /**< ntveff number of effective time varying variables */ |
int ntveff=0; /**< ntveff number of effective time varying variables */ |
| int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */ |
int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */ |
| int cptcov=0; /* Working variable */ |
int cptcov=0; /* Working variable */ |
| |
int nobs=10; /* Number of observations in the data lastobs-firstobs */ |
| int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */ |
int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */ |
| int npar=NPARMAX; |
int npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */ |
| int nlstate=2; /* Number of live states */ |
int nlstate=2; /* Number of live states */ |
| int ndeath=1; /* Number of dead states */ |
int ndeath=1; /* Number of dead states */ |
| int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ |
int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ |
|
Line 1202 double **pmmij, ***probs; /* Global poin
|
Line 1298 double **pmmij, ***probs; /* Global poin
|
| double ***mobaverage, ***mobaverages; /* New global variable */ |
double ***mobaverage, ***mobaverages; /* New global variable */ |
| double *ageexmed,*agecens; |
double *ageexmed,*agecens; |
| double dateintmean=0; |
double dateintmean=0; |
| |
double anprojd, mprojd, jprojd; /* For eventual projections */ |
| |
double anprojf, mprojf, jprojf; |
| |
|
| |
double anbackd, mbackd, jbackd; /* For eventual backprojections */ |
| |
double anbackf, mbackf, jbackf; |
| |
double jintmean,mintmean,aintmean; |
| double *weight; |
double *weight; |
| int **s; /* Status */ |
int **s; /* Status */ |
| double *agedc; |
double *agedc; |
|
Line 2321 void powell(double p[], double **xi, int
|
Line 2422 void powell(double p[], double **xi, int
|
| /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ |
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ |
| /* New value of last point Pn is not computed, P(n-1) */ |
/* New value of last point Pn is not computed, P(n-1) */ |
| for(j=1;j<=n;j++) { |
for(j=1;j<=n;j++) { |
| if(flatdir[j] >0){ |
if(flatdir[j] >0){ |
| printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
| fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
| } |
} |
| /* printf("\n"); */ |
/* printf("\n"); */ |
| /* fprintf(ficlog,"\n"); */ |
/* fprintf(ficlog,"\n"); */ |
| } |
} |
| /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */ |
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */ |
| if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */ |
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */ |
| /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ |
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ |
|
Line 2513 void powell(double p[], double **xi, int
|
Line 2614 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 2542 void powell(double p[], double **xi, int
|
Line 2646 void powell(double p[], double **xi, int
|
| double **newm; |
double **newm; |
| double agefin, delaymax=200. ; /* 100 Max number of years to converge */ |
double agefin, delaymax=200. ; /* 100 Max number of years to converge */ |
| int ncvloop=0; |
int ncvloop=0; |
| |
int first=0; |
| |
|
| min=vector(1,nlstate); |
min=vector(1,nlstate); |
| max=vector(1,nlstate); |
max=vector(1,nlstate); |
|
Line 2638 void powell(double p[], double **xi, int
|
Line 2743 void powell(double p[], double **xi, int
|
| free_vector(meandiff,1,nlstate); |
free_vector(meandiff,1,nlstate); |
| return prlim; |
return prlim; |
| } |
} |
| } /* age loop */ |
} /* agefin loop */ |
| /* After some age loop it doesn't converge */ |
/* After some age loop it doesn't converge */ |
| printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ |
if(!first){ |
| Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
first=1; |
| |
printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d). Others in log file only...\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
| |
} |
| |
fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
| |
|
| /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
| free_vector(min,1,nlstate); |
free_vector(min,1,nlstate); |
| free_vector(max,1,nlstate); |
free_vector(max,1,nlstate); |
|
Line 2707 Earliest age to start was %d-%d=%d, ncvl
|
Line 2816 Earliest age to start was %d-%d=%d, ncvl
|
| /* Even if hstepm = 1, at least one multiplication by the unit matrix */ |
/* Even if hstepm = 1, at least one multiplication by the unit matrix */ |
| /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */ |
/* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */ |
| /* for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */ |
/* for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */ |
| for(agefin=age; agefin<AGESUP; agefin=agefin+stepm/YEARM){ /* A changer en age */ |
/* for(agefin=age; agefin<AGESUP; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */ |
| |
for(agefin=age; agefin<FMIN(AGESUP,age+delaymax); agefin=agefin+stepm/YEARM){ /* A changer en age */ |
| ncvloop++; |
ncvloop++; |
| newm=savm; /* oldm should be kept from previous iteration or unity at start */ |
newm=savm; /* oldm should be kept from previous iteration or unity at start */ |
| /* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */ |
/* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */ |
|
Line 2821 Earliest age to start was %d-%d=%d, ncvl
|
Line 2931 Earliest age to start was %d-%d=%d, ncvl
|
| free_vector(meandiff,1,nlstate); |
free_vector(meandiff,1,nlstate); |
| return bprlim; |
return bprlim; |
| } |
} |
| } /* age loop */ |
} /* agefin loop */ |
| /* After some age loop it doesn't converge */ |
/* After some age loop it doesn't converge */ |
| if(first){ |
if(!first){ |
| first=1; |
first=1; |
| printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\ |
printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\ |
| Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
|
Line 2906 double **pmij(double **ps, double *cov,
|
Line 3016 double **pmij(double **ps, double *cov,
|
| ps[ii][ii]=1; |
ps[ii][ii]=1; |
| } |
} |
| } |
} |
| |
|
| |
|
| /* for(ii=1; ii<= nlstate+ndeath; ii++){ */ |
/* for(ii=1; ii<= nlstate+ndeath; ii++){ */ |
| /* for(jj=1; jj<= nlstate+ndeath; jj++){ */ |
/* for(jj=1; jj<= nlstate+ndeath; jj++){ */ |
| /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */ |
/* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */ |
|
Line 2927 double **pmij(double **ps, double *cov,
|
Line 3037 double **pmij(double **ps, double *cov,
|
| /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ |
/* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ |
| double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) |
double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) |
| { |
{ |
| /* Computes the backward probability at age agefin and covariate combination ij. In fact cov is already filled and x too. |
/* Computes the backward probability at age agefin, cov[2], and covariate combination 'ij'. In fact cov is already filled and x too. |
| * Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij. |
* Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij. |
| */ |
*/ |
| int i, ii, j,k; |
int i, ii, j,k; |
|
Line 2935 double **pmij(double **ps, double *cov,
|
Line 3045 double **pmij(double **ps, double *cov,
|
| double **out, **pmij(); |
double **out, **pmij(); |
| double sumnew=0.; |
double sumnew=0.; |
| double agefin; |
double agefin; |
| double k3=0.; /* constant of the w_x diagonal matrixe (in order for B to sum to 1 even for death state) */ |
double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */ |
| double **dnewm, **dsavm, **doldm; |
double **dnewm, **dsavm, **doldm; |
| double **bbmij; |
double **bbmij; |
| |
|
|
Line 2954 double **pmij(double **ps, double *cov,
|
Line 3064 double **pmij(double **ps, double *cov,
|
| /* outputs pmmij which is a stochastic matrix in row */ |
/* outputs pmmij which is a stochastic matrix in row */ |
| |
|
| /* Diag(w_x) */ |
/* Diag(w_x) */ |
| /* Problem with prevacurrent which can be zero */ |
/* Rescaling the cross-sectional prevalence: 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]); */ |
/* 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 */ |
|
Line 2981 double **pmij(double **ps, double *cov,
|
Line 3091 double **pmij(double **ps, double *cov,
|
| } |
} |
| /* End doldm, At the end doldm is diag[(w_i)] */ |
/* End doldm, At the end doldm is diag[(w_i)] */ |
| |
|
| /* left Product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm) */ |
/* Left product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm): diag[(w_i)*Px */ |
| bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* Bug Valgrind */ |
bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* was a Bug Valgrind */ |
| |
|
| /* Diag(Sum_i w^i_x p^ij_x */ |
/* Diag(Sum_i w^i_x p^ij_x, should be the prevalence at age x+stepm */ |
| /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */ |
/* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */ |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| sumnew=0.; |
sumnew=0.; |
|
Line 3002 double **pmij(double **ps, double *cov,
|
Line 3112 double **pmij(double **ps, double *cov,
|
| } /*End ii */ |
} /*End ii */ |
| } /* End j, At the end dsavm is diag[1/(w_1p1i+w_2 p2i)] for ALL states even if the sum is only for live states */ |
} /* End j, At the end dsavm is diag[1/(w_1p1i+w_2 p2i)] for ALL states even if the sum is only for live states */ |
| |
|
| ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* Bug Valgrind */ |
ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* was a Bug Valgrind */ |
| /* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ |
/* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ |
| /* end bmij */ |
/* end bmij */ |
| return ps; /*pointer is unchanged */ |
return ps; /*pointer is unchanged */ |
|
Line 3074 double **bpmij(double **ps, double *cov,
|
Line 3184 double **bpmij(double **ps, double *cov,
|
| ps[ii][ii]=1; |
ps[ii][ii]=1; |
| } |
} |
| } |
} |
| /* Added for backcast */ /* Transposed matrix too */ |
/* Added for prevbcast */ /* Transposed matrix too */ |
| for(jj=1; jj<= nlstate+ndeath; jj++){ |
for(jj=1; jj<= nlstate+ndeath; jj++){ |
| s1=0.; |
s1=0.; |
| for(ii=1; ii<= nlstate+ndeath; ii++){ |
for(ii=1; ii<= nlstate+ndeath; ii++){ |
|
Line 3834 return -l;
|
Line 3944 return -l;
|
| |
|
| |
|
| /*************** function likelione ***********/ |
/*************** function likelione ***********/ |
| void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) |
void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*func)(double [])) |
| { |
{ |
| /* This routine should help understanding what is done with |
/* This routine should help understanding what is done with |
| the selection of individuals/waves and |
the selection of individuals/waves and |
|
Line 3858 void likelione(FILE *ficres,double p[],
|
Line 3968 void likelione(FILE *ficres,double p[],
|
| fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n"); |
fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n"); |
| } |
} |
| |
|
| *fretone=(*funcone)(p); |
*fretone=(*func)(p); |
| if(*globpri !=0){ |
if(*globpri !=0){ |
| fclose(ficresilk); |
fclose(ficresilk); |
| if (mle ==0) |
if (mle ==0) |
|
Line 4333 void pstamp(FILE *fichier)
|
Line 4443 void pstamp(FILE *fichier)
|
| fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); |
fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); |
| } |
} |
| |
|
| |
void date2dmy(double date,double *day, double *month, double *year){ |
| |
double yp=0., yp1=0., yp2=0.; |
| |
|
| |
yp1=modf(date,&yp);/* extracts integral of date in yp and |
| |
fractional in yp1 */ |
| |
*year=yp; |
| |
yp2=modf((yp1*12),&yp); |
| |
*month=yp; |
| |
yp1=modf((yp2*30.5),&yp); |
| |
*day=yp; |
| |
if(*day==0) *day=1; |
| |
if(*month==0) *month=1; |
| |
} |
| |
|
| |
|
| |
|
| /************ Frequencies ********************/ |
/************ Frequencies ********************/ |
|
Line 4348 void freqsummary(char fileres[], double
|
Line 4472 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 4361 void freqsummary(char fileres[], double
|
Line 4485 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 4469 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4595 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 4485 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4614 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 4508 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4634 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 4528 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4654 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 4546 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4672 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 4561 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4692 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 4598 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4734 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 fixed quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]); |
| |
fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]); |
| |
if(weightopt==1){ |
| |
printf(" Weighted mean and standard deviation of"); |
| |
fprintf(ficlog," Weighted mean and standard deviation of"); |
| |
fprintf(ficresphtmfr," Weighted mean and standard deviation of"); |
| |
} |
| |
printf(" fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
| |
fprintf(ficlog," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
| |
fprintf(ficresphtmfr," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)<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 4832 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4989 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 4874 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5031 Title=%s <br>Datafile=%s Firstpass=%d La
|
| } |
} |
| } /* end mle=-2 */ |
} /* end mle=-2 */ |
| dateintmean=dateintsum/k2cpt; |
dateintmean=dateintsum/k2cpt; |
| |
date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); |
| |
|
| 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 4985 void prevalence(double ***probs, double
|
Line 5145 void prevalence(double ***probs, double
|
| /*j=cptcoveff;*/ |
/*j=cptcoveff;*/ |
| if (cptcovn<1) {j=1;ncodemax[1]=1;} |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
| |
|
| first=1; |
first=0; |
| for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ |
for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ |
| for (i=1; i<=nlstate; i++) |
for (i=1; i<=nlstate; i++) |
| for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) |
for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) |
|
Line 5043 void prevalence(double ***probs, double
|
Line 5203 void prevalence(double ***probs, double
|
| if(posprop>=1.e-5){ |
if(posprop>=1.e-5){ |
| probs[i][jk][j1]= prop[jk][i]/posprop; |
probs[i][jk][j1]= prop[jk][i]/posprop; |
| } else{ |
} else{ |
| if(first==1){ |
if(!first){ |
| first=0; |
first=1; |
| printf("Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); |
printf("Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); |
| fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); |
|
| }else{ |
}else{ |
| fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); |
fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases.\n",jk,i,jk, j1,probs[i][jk][j1]); |
| } |
} |
| } |
} |
| } |
} |
|
Line 5066 void prevalence(double ***probs, double
|
Line 5225 void prevalence(double ***probs, double
|
| |
|
| void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm) |
void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm) |
| { |
{ |
| /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i. |
/* Concatenates waves: wav[i] is the number of effective (useful waves in the sense that a non interview is useless) of individual i. |
| Death is a valid wave (if date is known). |
Death is a valid wave (if date is known). |
| mw[mi][i] is the mi (mi=1 to wav[i]) effective wave of individual i |
mw[mi][i] is the mi (mi=1 to wav[i]) effective wave of individual i |
| dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] |
dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] |
| and mw[mi+1][i]. dh depends on stepm. |
and mw[mi+1][i]. dh depends on stepm. s[m][i] exists for any wave from firstpass to lastpass |
| */ |
*/ |
| |
|
| int i=0, mi=0, m=0, mli=0; |
int i=0, mi=0, m=0, mli=0; |
|
Line 5108 void concatwav(int wav[], int **dh, int
|
Line 5267 void concatwav(int wav[], int **dh, int
|
| #else |
#else |
| if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ |
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ |
| if(firsthree == 0){ |
if(firsthree == 0){ |
| printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
| firsthree=1; |
firsthree=1; |
| } |
} |
| fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
| mw[++mi][i]=m; |
mw[++mi][i]=m; |
| mli=m; |
mli=m; |
| } |
} |
|
Line 5294 void concatwav(int wav[], int **dh, int
|
Line 5453 void concatwav(int wav[], int **dh, int
|
| /* *cptcov=0; */ |
/* *cptcov=0; */ |
| |
|
| for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ |
for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ |
| |
for (k=1; k <= maxncov; k++) |
| |
for(j=1; j<=2; j++) |
| |
nbcode[k][j]=0; /* Valgrind */ |
| |
|
| /* Loop on covariates without age and products and no quantitative variable */ |
/* Loop on covariates without age and products and no quantitative variable */ |
| /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */ |
|
| for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ |
for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ |
| for (j=-1; (j < maxncov); j++) Ndum[j]=0; |
for (j=-1; (j < maxncov); j++) Ndum[j]=0; |
| if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
|
Line 5314 void concatwav(int wav[], int **dh, int
|
Line 5475 void concatwav(int wav[], int **dh, int
|
| modmaxcovj=ij; |
modmaxcovj=ij; |
| else if (ij < modmincovj) |
else if (ij < modmincovj) |
| modmincovj=ij; |
modmincovj=ij; |
| if ((ij < -1) && (ij > NCOVMAX)){ |
if (ij <0 || ij >1 ){ |
| |
printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); |
| |
fprintf(ficlog,"Information, currently IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); |
| |
} |
| |
if ((ij < -1) || (ij > NCOVMAX)){ |
| printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
| exit(1); |
exit(1); |
| }else |
}else |
|
Line 5360 void concatwav(int wav[], int **dh, int
|
Line 5525 void concatwav(int wav[], int **dh, int
|
| /* nbcode[Tvar[j]][3]=2; */ |
/* nbcode[Tvar[j]][3]=2; */ |
| /* To be continued (not working yet). */ |
/* To be continued (not working yet). */ |
| ij=0; /* ij is similar to i but can jump over null modalities */ |
ij=0; /* ij is similar to i but can jump over null modalities */ |
| for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ |
|
| |
/* for (i=modmincovj; i<=modmaxcovj; i++) { */ /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ |
| |
/* Skipping the case of missing values by reducing nbcode to 0 and 1 and not -1, 0, 1 */ |
| |
/* model=V1+V2+V3, if V2=-1, 0 or 1, then nbcode[2][1]=0 and nbcode[2][2]=1 instead of |
| |
* nbcode[2][1]=-1, nbcode[2][2]=0 and nbcode[2][3]=1 */ |
| |
/*, could be restored in the future */ |
| |
for (i=0; i<=1; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ |
| if (Ndum[i] == 0) { /* If nobody responded to this modality k */ |
if (Ndum[i] == 0) { /* If nobody responded to this modality k */ |
| break; |
break; |
| } |
} |
| ij++; |
ij++; |
| nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/ |
nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1 . Could be -1*/ |
| cptcode = ij; /* New max modality for covar j */ |
cptcode = ij; /* New max modality for covar j */ |
| } /* end of loop on modality i=-1 to 1 or more */ |
} /* end of loop on modality i=-1 to 1 or more */ |
| break; |
break; |
|
Line 5381 void concatwav(int wav[], int **dh, int
|
Line 5552 void concatwav(int wav[], int **dh, int
|
| break; |
break; |
| } /* end switch */ |
} /* end switch */ |
| } /* end dummy test */ |
} /* end dummy test */ |
| |
} /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ |
| /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ |
|
| /* /\*recode from 0 *\/ */ |
|
| /* k is a modality. If we have model=V1+V1*sex */ |
|
| /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ |
|
| /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ |
|
| /* } */ |
|
| /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ |
|
| /* if (ij > ncodemax[j]) { */ |
|
| /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
|
| /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
|
| /* break; */ |
|
| /* } */ |
|
| /* } /\* end of loop on modality k *\/ */ |
|
| } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ |
|
| |
|
| for (k=-1; k< maxncov; k++) Ndum[k]=0; |
for (k=-1; k< maxncov; k++) Ndum[k]=0; |
| /* Look at fixed dummy (single or product) covariates to check empty modalities */ |
/* Look at fixed dummy (single or product) covariates to check empty modalities */ |
|
Line 5769 void concatwav(int wav[], int **dh, int
|
Line 5926 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; |
| double **dnewmp,**doldmp; |
double **dnewmp,**doldmp; |
| int i, j, nhstepm, hstepm, h, nstepm ; |
int i, j, nhstepm, hstepm, h, nstepm ; |
| |
int first=0; |
| 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 5845 void concatwav(int wav[], int **dh, int
|
Line 6004 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 5900 void concatwav(int wav[], int **dh, int
|
Line 6059 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 5912 void concatwav(int wav[], int **dh, int
|
Line 6074 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); |
| |
|
| 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) { |
if (popbased==1) { |
|
Line 5961 void concatwav(int wav[], int **dh, int
|
Line 6128 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 5984 void concatwav(int wav[], int **dh, int
|
Line 6155 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 6001 void concatwav(int wav[], int **dh, int
|
Line 6178 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 6104 void concatwav(int wav[], int **dh, int
|
Line 6285 void concatwav(int wav[], int **dh, int
|
| int theta; |
int theta; |
| |
|
| pstamp(ficresvpl); |
pstamp(ficresvpl); |
| fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); |
fprintf(ficresvpl,"# Standard deviation of period (forward stable) prevalences \n"); |
| fprintf(ficresvpl,"# Age "); |
fprintf(ficresvpl,"# Age "); |
| if(nresult >=1) |
if(nresult >=1) |
| fprintf(ficresvpl," Result# "); |
fprintf(ficresvpl," Result# "); |
|
Line 6133 void concatwav(int wav[], int **dh, int
|
Line 6314 void concatwav(int wav[], int **dh, int
|
| for(i=1; i<=npar; i++){ /* Computes gradient */ |
for(i=1; i<=npar; i++){ /* Computes gradient */ |
| xp[i] = x[i] + (i==theta ?delti[theta]:0); |
xp[i] = x[i] + (i==theta ?delti[theta]:0); |
| } |
} |
| if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) |
/* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ |
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); |
/* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ |
| else |
/* else */ |
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); |
| for(i=1;i<=nlstate;i++){ |
for(i=1;i<=nlstate;i++){ |
| gp[i] = prlim[i][i]; |
gp[i] = prlim[i][i]; |
| mgp[theta][i] = prlim[i][i]; |
mgp[theta][i] = prlim[i][i]; |
| } |
} |
| for(i=1; i<=npar; i++) /* Computes gradient */ |
for(i=1; i<=npar; i++) /* Computes gradient */ |
| xp[i] = x[i] - (i==theta ?delti[theta]:0); |
xp[i] = x[i] - (i==theta ?delti[theta]:0); |
| if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) |
/* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ |
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); |
/* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ |
| else |
/* else */ |
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); |
| for(i=1;i<=nlstate;i++){ |
for(i=1;i<=nlstate;i++){ |
| gm[i] = prlim[i][i]; |
gm[i] = prlim[i][i]; |
| mgm[theta][i] = prlim[i][i]; |
mgm[theta][i] = prlim[i][i]; |
|
Line 6195 void concatwav(int wav[], int **dh, int
|
Line 6376 void concatwav(int wav[], int **dh, int
|
| fprintf(ficresvpl,"%.0f ",age ); |
fprintf(ficresvpl,"%.0f ",age ); |
| if(nresult >=1) |
if(nresult >=1) |
| fprintf(ficresvpl,"%d ",nres ); |
fprintf(ficresvpl,"%d ",nres ); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++){ |
| fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); |
fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); |
| |
/* for(j=1;j<=nlstate;j++) */ |
| |
/* fprintf(ficresvpl," %d %.5f ",j,prlim[j][i]); */ |
| |
} |
| fprintf(ficresvpl,"\n"); |
fprintf(ficresvpl,"\n"); |
| free_vector(gp,1,nlstate); |
free_vector(gp,1,nlstate); |
| free_vector(gm,1,nlstate); |
free_vector(gm,1,nlstate); |
|
Line 6413 void varprob(char optionfilefiname[], do
|
Line 6597 void varprob(char optionfilefiname[], do
|
| fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
| fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
| |
|
| fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. %s</li>\n",optionfilehtmcov,optionfilehtmcov); |
fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. File %s</li>\n",optionfilehtmcov,optionfilehtmcov); |
| fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov); |
fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov); |
| fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \ |
fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \ |
| and drawn. It helps understanding how is the covariance between two incidences.\ |
and drawn. It helps understanding how is the covariance between two incidences.\ |
|
Line 6610 To be simple, these graphs help to under
|
Line 6794 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 6641 To be simple, these graphs help to under
|
Line 6830 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 6678 To be simple, these graphs help to under
|
Line 6867 To be simple, these graphs help to under
|
| void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ |
void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ |
| 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 prevbcast, int estepm , \ |
| double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \ |
double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \ |
| double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){ |
double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, 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 6701 void printinghtml(char fileresu[], char
|
Line 6890 void printinghtml(char fileresu[], char
|
| - Estimated back transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", |
- Estimated back transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", |
| stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_")); |
stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
- Period (forward) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); |
subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Period (stable) back prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
- Backward prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_")); |
subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ |
- (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ |
|
Line 6810 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 6999 divided by h: <sub>h</sub>P<sub>ij</sub>
|
| <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); |
<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); |
| /* Survival functions (period) in state j */ |
/* Survival functions (period) in state j */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); |
<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); |
| } |
} |
| /* State specific survival functions (period) */ |
/* State specific survival functions (period) */ |
| 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 in state %d and in any other live state (total).\ |
| Or probability to survive in various states (1 to %d) being in state %d at different ages. \ |
And probability to be observed in various states (up 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 (forward stable) prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <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,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); |
<img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); |
| } |
} |
| if(backcast==1){ |
if(prevbcast==1){ |
| /* Period (stable) back prevalence in each health state */ |
/* Backward prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
| <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); |
<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); |
| } |
} |
| } |
} |
| if(prevfcast==1){ |
if(prevfcast==1){ |
| /* Projection of prevalence up to period (stable) prevalence in each health state */ |
/* Projection of prevalence up to period (forward 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), 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> \ |
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) forward 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, dateproj1, dateproj2, 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, dateprojd, dateprojf, 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(prevbcast==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), \ |
fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \ |
|
Line 6887 See page 'Matrix of variance-covariance
|
Line 7076 See page 'Matrix of variance-covariance
|
| <a href=\"%s\">%s</a> <br>\n</li>", |
<a href=\"%s\">%s</a> <br>\n</li>", |
| estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); |
estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n", |
- Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the forward (period) prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n", |
| estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_")); |
estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n", |
- Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n", |
| estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_")); |
estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\ |
- Standard deviation of forward (period) prevalences: <a href=\"%s\">%s</a> <br>\n",\ |
| subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); |
subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); |
| |
|
| /* if(popforecast==1) fprintf(fichtm,"\n */ |
/* if(popforecast==1) fprintf(fichtm,"\n */ |
|
Line 6951 true period expectancies (those weighted
|
Line 7140 true period expectancies (those weighted
|
| } |
} |
| |
|
| /******************* Gnuplot file **************/ |
/******************* Gnuplot file **************/ |
| 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){ |
void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, 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 7029 void printinggnuplot(char fileresu[], ch
|
Line 7218 void printinggnuplot(char fileresu[], ch
|
| continue; |
continue; |
| /* We are interested in selected combination by the resultline */ |
/* We are interested in selected combination by the resultline */ |
| /* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */ |
/* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */ |
| fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); |
fprintf(ficgp,"\n# 1st: Forward (stable period) prevalence with CI: 'VPL_' files and live state =%d ", cpt); |
| strcpy(gplotlabel,"("); |
strcpy(gplotlabel,"("); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ |
|
Line 7067 void printinggnuplot(char fileresu[], ch
|
Line 7256 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\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); |
fprintf(ficgp,"\" t\"Forward prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),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)"); |
|
Line 7105 void printinggnuplot(char fileresu[], ch
|
Line 7294 void printinggnuplot(char fileresu[], ch
|
| } /* end covariate */ |
} /* end covariate */ |
| } /* end if no covariate */ |
} /* end if no covariate */ |
| |
|
| if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
if(prevbcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
| /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
| fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */ |
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */ |
| if(cptcoveff ==0){ |
if(cptcoveff ==0){ |
|
Line 7132 void printinggnuplot(char fileresu[], ch
|
Line 7321 void printinggnuplot(char fileresu[], ch
|
| } |
} |
| } /* end covariate */ |
} /* end covariate */ |
| } /* end if no covariate */ |
} /* end if no covariate */ |
| if(backcast == 1){ |
if(prevbcast == 1){ |
| fprintf(ficgp,", \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); |
fprintf(ficgp,", \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); |
| /* k1-1 error should be nres-1*/ |
/* k1-1 error should be nres-1*/ |
| for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
|
Line 7151 void printinggnuplot(char fileresu[], ch
|
Line 7340 void printinggnuplot(char fileresu[], ch
|
| } |
} |
| fprintf(ficgp,"\" t\"\" w l lt 4"); |
fprintf(ficgp,"\" t\"\" w l lt 4"); |
| } /* end if backprojcast */ |
} /* end if backprojcast */ |
| } /* end if backcast */ |
} /* end if prevbcast */ |
| /* fprintf(ficgp,"\nset out ;unset label;\n"); */ |
/* fprintf(ficgp,"\nset out ;unset label;\n"); */ |
| fprintf(ficgp,"\nset out ;unset title;\n"); |
fprintf(ficgp,"\nset out ;unset title;\n"); |
| } /* nres */ |
} /* nres */ |
|
Line 7396 set ter svg size 640, 480\nunset log y\n
|
Line 7585 set ter svg size 640, 480\nunset log y\n
|
| continue; |
continue; |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */ |
| strcpy(gplotlabel,"("); |
strcpy(gplotlabel,"("); |
| fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n#CV preval stable (forward): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| /* 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 7439 set ter svg size 640, 480\nunset log y\n
|
Line 7628 set ter svg size 640, 480\nunset log y\n
|
| |
|
| |
|
| /* 7eme */ |
/* 7eme */ |
| if(backcast == 1){ |
if(prevbcast == 1){ |
| /* CV back preval stable (period) for each covariate */ |
/* CV backward prevalence for each covariate */ |
| for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ |
for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ |
| for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
| if(m != 1 && TKresult[nres]!= k1) |
if(m != 1 && TKresult[nres]!= k1) |
| continue; |
continue; |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */ |
| strcpy(gplotlabel,"("); |
strcpy(gplotlabel,"("); |
| fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n#CV Backward stable prevalence: 'pijb' files, covariatecombination#=%d state=%d",k1, cpt); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| /* 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 7491 set ter svg size 640, 480\nunset log y\n
|
Line 7680 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficgp,"\nset out; unset label;\n"); |
fprintf(ficgp,"\nset out; unset label;\n"); |
| } /* end cpt state*/ |
} /* end cpt state*/ |
| } /* end covariate */ |
} /* end covariate */ |
| } /* End if backcast */ |
} /* End if prevbcast */ |
| |
|
| /* 8eme */ |
/* 8eme */ |
| if(prevfcast==1){ |
if(prevfcast==1){ |
| /* Projection from cross-sectional to stable (period) for each covariate */ |
/* Projection from cross-sectional to forward stable (period) prevalence for each covariate */ |
| |
|
| for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ |
for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ |
| for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
|
Line 7503 set ter svg size 640, 480\nunset log y\n
|
Line 7692 set ter svg size 640, 480\nunset log y\n
|
| continue; |
continue; |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| strcpy(gplotlabel,"("); |
strcpy(gplotlabel,"("); |
| fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n#Projection of prevalence to forward stable prevalence (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); |
| 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 7607 set ter svg size 640, 480\nunset log y\n
|
Line 7796 set ter svg size 640, 480\nunset log y\n
|
| } /* end covariate */ |
} /* end covariate */ |
| } /* End if prevfcast */ |
} /* End if prevfcast */ |
| |
|
| if(backcast==1){ |
if(prevbcast==1){ |
| /* Back projection from cross-sectional to stable (mixed) for each covariate */ |
/* Back projection from cross-sectional to stable (mixed) for each covariate */ |
| |
|
| for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ |
for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ |
|
Line 7720 set ter svg size 640, 480\nunset log y\n
|
Line 7909 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficgp,"\nset out; unset label;\n"); |
fprintf(ficgp,"\nset out; unset label;\n"); |
| } /* end cpt state*/ |
} /* end cpt state*/ |
| } /* end covariate */ |
} /* end covariate */ |
| } /* End if backcast */ |
} /* End if prevbcast */ |
| |
|
| |
|
| /* 9eme writing MLE parameters */ |
/* 9eme writing MLE parameters */ |
|
Line 7937 set ter svg size 640, 480\nunset log y\n
|
Line 8126 set ter svg size 640, 480\nunset log y\n
|
| int modcovmax =1; |
int modcovmax =1; |
| int mobilavrange, mob; |
int mobilavrange, mob; |
| int iage=0; |
int iage=0; |
| |
int firstA1=0, firstA2=0; |
| |
|
| double sum=0., sumr=0.; |
double sum=0., sumr=0.; |
| double age; |
double age; |
|
Line 7945 set ter svg size 640, 480\nunset log y\n
|
Line 8135 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 8034 set ter svg size 640, 480\nunset log y\n
|
Line 8224 set ter svg size 640, 480\nunset log y\n
|
| } /* age */ |
} /* age */ |
| /* Thus we have agemingood and agemaxgood as well as goodr for raw (preobs) */ |
/* Thus we have agemingood and agemaxgood as well as goodr for raw (preobs) */ |
| /* but they will change */ |
/* but they will change */ |
| |
firstA1=0;firstA2=0; |
| for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, filling up to the youngest */ |
for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, filling up to the youngest */ |
| sumnewm[cptcod]=0.; |
sumnewm[cptcod]=0.; |
| sumnewmr[cptcod]=0.; |
sumnewmr[cptcod]=0.; |
|
Line 8066 set ter svg size 640, 480\nunset log y\n
|
Line 8257 set ter svg size 640, 480\nunset log y\n
|
| sumr+=probs[(int)age][i][cptcod]; |
sumr+=probs[(int)age][i][cptcod]; |
| } |
} |
| if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
| printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); |
if(!firstA1){ |
| |
firstA1=1; |
| |
printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage); |
| |
} |
| |
fprintf(ficlog,"Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); |
| } /* end bad */ |
} /* end bad */ |
| /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
| if(fabs(sumr - 1.) > 1.e-3) { /* bad */ |
if(fabs(sumr - 1.) > 1.e-3) { /* bad */ |
| printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); |
if(!firstA2){ |
| |
firstA2=1; |
| |
printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage); |
| |
} |
| |
fprintf(ficlog,"Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); |
| } /* end bad */ |
} /* end bad */ |
| }/* age */ |
}/* age */ |
| |
|
|
Line 8158 set ter svg size 640, 480\nunset log y\n
|
Line 8357 set ter svg size 640, 480\nunset log y\n
|
| }/* End movingaverage */ |
}/* End movingaverage */ |
| |
|
| |
|
| |
|
| /************** Forecasting ******************/ |
/************** Forecasting ******************/ |
| 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){ |
/* void prevforecast(char fileres[], double dateintmean, double anprojd, double mprojd, double jprojd, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double anprojf, double p[], int cptcoveff)*/ |
| /* proj1, year, month, day of starting projection |
void prevforecast(char fileres[], double dateintmean, double dateprojd, double dateprojf, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double p[], int cptcoveff){ |
| |
/* dateintemean, mean date of interviews |
| |
dateprojd, year, month, day of starting projection |
| |
dateprojf date of end of projection;year of end of projection (same day and month as proj1). |
| 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 |
| anproj2 year of en of projection (same day and month as proj1). |
|
| */ |
*/ |
| |
/* double anprojd, mprojd, jprojd; */ |
| |
/* double anprojf, mprojf, jprojf; */ |
| int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; |
int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; |
| double agec; /* generic age */ |
double agec; /* generic age */ |
| double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; |
double agelim, ppij, yp,yp1,yp2; |
| double *popeffectif,*popcount; |
double *popeffectif,*popcount; |
| double ***p3mat; |
double ***p3mat; |
| /* double ***mobaverage; */ |
/* double ***mobaverage; */ |
|
Line 8204 set ter svg size 640, 480\nunset log y\n
|
Line 8408 set ter svg size 640, 480\nunset log y\n
|
| if(estepm > stepm){ /* Yes every two year */ |
if(estepm > stepm){ /* Yes every two year */ |
| stepsize=2; |
stepsize=2; |
| } |
} |
| |
hstepm=hstepm/stepm; |
| |
|
| |
|
| |
/* yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp and */ |
| |
/* fractional in yp1 *\/ */ |
| |
/* aintmean=yp; */ |
| |
/* yp2=modf((yp1*12),&yp); */ |
| |
/* mintmean=yp; */ |
| |
/* yp1=modf((yp2*30.5),&yp); */ |
| |
/* jintmean=yp; */ |
| |
/* if(jintmean==0) jintmean=1; */ |
| |
/* if(mintmean==0) mintmean=1; */ |
| |
|
| hstepm=hstepm/stepm; |
|
| yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and |
|
| fractional in yp1 */ |
|
| anprojmean=yp; |
|
| yp2=modf((yp1*12),&yp); |
|
| mprojmean=yp; |
|
| yp1=modf((yp2*30.5),&yp); |
|
| jprojmean=yp; |
|
| if(jprojmean==0) jprojmean=1; |
|
| if(mprojmean==0) jprojmean=1; |
|
| |
|
| |
/* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */ |
| |
/* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */ |
| |
/* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */ |
| i1=pow(2,cptcoveff); |
i1=pow(2,cptcoveff); |
| if (cptcovn < 1){i1=1;} |
if (cptcovn < 1){i1=1;} |
| |
|
| fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); |
fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); |
| |
|
| fprintf(ficresf,"#****** Routine prevforecast **\n"); |
fprintf(ficresf,"#****** Routine prevforecast **\n"); |
| |
|
|
Line 8245 set ter svg size 640, 480\nunset log y\n
|
Line 8454 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficresf," p%d%d",i,j); |
fprintf(ficresf," p%d%d",i,j); |
| fprintf(ficresf," wp.%d",j); |
fprintf(ficresf," wp.%d",j); |
| } |
} |
| for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { |
for (yearp=0; yearp<=(anprojf-anprojd);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 ",jprojd,mprojd,anprojd+yearp); |
| /* for (agec=fage; agec>=(ageminpar-1); agec--){ */ |
/* for (agec=fage; agec>=(ageminpar-1); agec--){ */ |
| for (agec=fage; agec>=(bage); agec--){ |
for (agec=fage; agec>=(bage); agec--){ |
| nhstepm=(int) rint((agelim-agec)*YEARM/stepm); |
nhstepm=(int) rint((agelim-agec)*YEARM/stepm); |
|
Line 8265 set ter svg size 640, 480\nunset log y\n
|
Line 8474 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficresf,"\n"); |
fprintf(ficresf,"\n"); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); |
fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm); |
| |
|
| 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 8293 set ter svg size 640, 480\nunset log y\n
|
Line 8502 set ter svg size 640, 480\nunset log y\n
|
| } |
} |
| |
|
| /************** 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 |
void prevbackforecast(char fileres[], double ***prevacurrent, double dateintmean, double dateprojd, double dateprojf, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double p[], int cptcoveff){ |
| |
/* back1, year, month, day of starting backprojection |
| 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 end of backprojection (same day and month as back1). |
anback2 year of end of backprojection (same day and month as back1). |
|
Line 8302 set ter svg size 640, 480\nunset log y\n
|
Line 8512 set ter svg size 640, 480\nunset log y\n
|
| */ |
*/ |
| int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; |
int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; |
| double agec; /* generic age */ |
double agec; /* generic age */ |
| double agelim, ppij, ppi, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; |
double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/ |
| double *popeffectif,*popcount; |
double *popeffectif,*popcount; |
| double ***p3mat; |
double ***p3mat; |
| /* double ***mobaverage; */ |
/* double ***mobaverage; */ |
|
Line 8345 set ter svg size 640, 480\nunset log y\n
|
Line 8555 set ter svg size 640, 480\nunset log y\n
|
| } |
} |
| |
|
| 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 */ |
| fractional in yp1 */ |
/* fractional in yp1 *\/ */ |
| anprojmean=yp; |
/* aintmean=yp; */ |
| yp2=modf((yp1*12),&yp); |
/* yp2=modf((yp1*12),&yp); */ |
| mprojmean=yp; |
/* mintmean=yp; */ |
| yp1=modf((yp2*30.5),&yp); |
/* yp1=modf((yp2*30.5),&yp); */ |
| jprojmean=yp; |
/* jintmean=yp; */ |
| if(jprojmean==0) jprojmean=1; |
/* if(jintmean==0) jintmean=1; */ |
| if(mprojmean==0) jprojmean=1; |
/* if(mintmean==0) jintmean=1; */ |
| |
|
| i1=pow(2,cptcoveff); |
i1=pow(2,cptcoveff); |
| if (cptcovn < 1){i1=1;} |
if (cptcovn < 1){i1=1;} |
| |
|
| fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); |
fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); |
| printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); |
printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); |
| |
|
| fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); |
fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); |
| |
|
|
Line 8384 set ter svg size 640, 480\nunset log y\n
|
Line 8594 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficresfb," b%d%d",i,j); |
fprintf(ficresfb," b%d%d",i,j); |
| fprintf(ficresfb," b.%d",j); |
fprintf(ficresfb," b.%d",j); |
| } |
} |
| for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) { |
for (yearp=0; yearp>=(anbackf-anbackd);yearp -=stepsize) { |
| /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */ |
/* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */ |
| fprintf(ficresfb,"\n"); |
fprintf(ficresfb,"\n"); |
| fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); |
fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jbackd,mbackd,anbackd+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 */ |
for (agec=bage; agec<=fage; agec++){ /* testing */ |
|
Line 8410 set ter svg size 640, 480\nunset log y\n
|
Line 8620 set ter svg size 640, 480\nunset log y\n
|
| fprintf(ficresfb,"\n"); |
fprintf(ficresfb,"\n"); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec-h*hstepm/YEARM*stepm); |
fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm); |
| for(i=1; i<=nlstate+ndeath;i++) { |
for(i=1; i<=nlstate+ndeath;i++) { |
| ppij=0.;ppi=0.; |
ppij=0.;ppi=0.; |
| for(j=1; j<=nlstate;j++) { |
for(j=1; j<=nlstate;j++) { |
|
Line 8445 set ter svg size 640, 480\nunset log y\n
|
Line 8655 set ter svg size 640, 480\nunset log y\n
|
| |
|
| /* Variance of prevalence limit: varprlim */ |
/* 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){ |
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------*/ |
/*------- Variance of forward period (stable) prevalence------*/ |
| |
|
| char fileresvpl[FILENAMELENGTH]; |
char fileresvpl[FILENAMELENGTH]; |
| FILE *ficresvpl; |
FILE *ficresvpl; |
|
Line 8456 set ter svg size 640, 480\nunset log y\n
|
Line 8666 set ter svg size 640, 480\nunset log y\n
|
| strcpy(fileresvpl,"VPL_"); |
strcpy(fileresvpl,"VPL_"); |
| strcat(fileresvpl,fileresu); |
strcat(fileresvpl,fileresu); |
| if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { |
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { |
| printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); |
printf("Problem with variance of forward period (stable) prevalence resultfile: %s\n", fileresvpl); |
| exit(0); |
exit(0); |
| } |
} |
| printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); |
printf("Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); |
| fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); |
fprintf(ficlog, "Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); |
| |
|
| /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ |
|
Line 8497 set ter svg size 640, 480\nunset log y\n
|
Line 8707 set ter svg size 640, 480\nunset log y\n
|
| } |
} |
| |
|
| fclose(ficresvpl); |
fclose(ficresvpl); |
| printf("done variance-covariance of period prevalence\n");fflush(stdout); |
printf("done variance-covariance of forward period prevalence\n");fflush(stdout); |
| fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); |
fprintf(ficlog,"done variance-covariance of forward period prevalence\n");fflush(ficlog); |
| |
|
| } |
} |
| /* Variance of back prevalence: varbprlim */ |
/* Variance of back prevalence: varbprlim */ |
|
Line 8858 void prwizard(int ncovmodel, int nlstate
|
Line 9068 void prwizard(int ncovmodel, int nlstate
|
| /******************* Gompertz Likelihood ******************************/ |
/******************* Gompertz Likelihood ******************************/ |
| double gompertz(double x[]) |
double gompertz(double x[]) |
| { |
{ |
| double A,B,L=0.0,sump=0.,num=0.; |
double A=0.0,B=0.,L=0.0,sump=0.,num=0.; |
| int i,n=0; /* n is the size of the sample */ |
int i,n=0; /* n is the size of the sample */ |
| |
|
| for (i=1;i<=imx ; i++) { |
for (i=1;i<=imx ; i++) { |
|
Line 8866 double gompertz(double x[])
|
Line 9076 double gompertz(double x[])
|
| /* sump=sump+1;*/ |
/* sump=sump+1;*/ |
| num=num+1; |
num=num+1; |
| } |
} |
| |
L=0.0; |
| |
/* agegomp=AGEGOMP; */ |
| /* for (i=0; i<=imx; i++) |
/* for (i=0; i<=imx; i++) |
| if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ |
if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ |
| |
|
| for (i=1;i<=imx ; i++) |
for (i=1;i<=imx ; i++) { |
| { |
/* mu(a)=mu(agecomp)*exp(teta*(age-agegomp)) |
| if (cens[i] == 1 && wav[i]>1) |
mu(a)=x[1]*exp(x[2]*(age-agegomp)); x[1] and x[2] are per year. |
| A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); |
* L= Product mu(agedeces)exp(-\int_ageexam^agedc mu(u) du ) for a death between agedc (in month) |
| |
* and agedc +1 month, cens[i]=0: log(x[1]/YEARM) |
| if (cens[i] == 0 && wav[i]>1) |
* + |
| |
* exp(-\int_ageexam^agecens mu(u) du ) when censored, cens[i]=1 |
| |
*/ |
| |
if (wav[i] > 1 || agedc[i] < AGESUP) { |
| |
if (cens[i] == 1){ |
| |
A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); |
| |
} else if (cens[i] == 0){ |
| A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) |
A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) |
| +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); |
+log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); |
| |
} else |
| |
printf("Gompertz cens[%d] neither 1 nor 0\n",i); |
| /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ |
/*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ |
| if (wav[i] > 1 ) { /* ??? */ |
L=L+A*weight[i]; |
| L=L+A*weight[i]; |
|
| /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ |
/* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ |
| } |
} |
| } |
} |
| |
|
| /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ |
/*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ |
| |
|
| return -2*L*num/sump; |
return -2*L*num/sump; |
| } |
} |
|
Line 8896 double gompertz(double x[])
|
Line 9112 double gompertz(double x[])
|
| /******************* Gompertz_f Likelihood ******************************/ |
/******************* Gompertz_f Likelihood ******************************/ |
| double gompertz_f(const gsl_vector *v, void *params) |
double gompertz_f(const gsl_vector *v, void *params) |
| { |
{ |
| double A,B,LL=0.0,sump=0.,num=0.; |
double A=0.,B=0.,LL=0.0,sump=0.,num=0.; |
| double *x= (double *) v->data; |
double *x= (double *) v->data; |
| int i,n=0; /* n is the size of the sample */ |
int i,n=0; /* n is the size of the sample */ |
| |
|
|
Line 8989 int readdata(char datafile[], int firsto
|
Line 9205 int readdata(char datafile[], int firsto
|
| int i=0, j=0, n=0, iv=0, v; |
int i=0, j=0, n=0, iv=0, v; |
| int lstra; |
int lstra; |
| int linei, month, year,iout; |
int linei, month, year,iout; |
| |
int noffset=0; /* This is the offset if BOM data file */ |
| char line[MAXLINE], linetmp[MAXLINE]; |
char line[MAXLINE], linetmp[MAXLINE]; |
| char stra[MAXLINE], strb[MAXLINE]; |
char stra[MAXLINE], strb[MAXLINE]; |
| char *stratrunc; |
char *stratrunc; |
|
Line 9022 int readdata(char datafile[], int firsto
|
Line 9239 int readdata(char datafile[], int firsto
|
| fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1; |
fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1; |
| } |
} |
| |
|
| i=1; |
/* Is it a BOM UTF-8 Windows file? */ |
| |
/* First data line */ |
| linei=0; |
linei=0; |
| |
while(fgets(line, MAXLINE, fic)) { |
| |
noffset=0; |
| |
if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */ |
| |
{ |
| |
noffset=noffset+3; |
| |
printf("# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile);fflush(stdout); |
| |
fprintf(ficlog,"# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile); |
| |
fflush(ficlog); return 1; |
| |
} |
| |
/* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/ |
| |
else if( line[0] == (char)0xFF && line[1] == (char)0xFE) |
| |
{ |
| |
noffset=noffset+2; |
| |
printf("# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); |
| |
fprintf(ficlog,"# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); |
| |
fflush(ficlog); return 1; |
| |
} |
| |
else if( line[0] == 0 && line[1] == 0) |
| |
{ |
| |
if( line[2] == (char)0xFE && line[3] == (char)0xFF){ |
| |
noffset=noffset+4; |
| |
printf("# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); |
| |
fprintf(ficlog,"# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); |
| |
fflush(ficlog); return 1; |
| |
} |
| |
} else{ |
| |
;/*printf(" Not a BOM file\n");*/ |
| |
} |
| |
/* If line starts with a # it is a comment */ |
| |
if (line[noffset] == '#') { |
| |
linei=linei+1; |
| |
break; |
| |
}else{ |
| |
break; |
| |
} |
| |
} |
| |
fclose(fic); |
| |
if((fic=fopen(datafile,"r"))==NULL) { |
| |
printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout); |
| |
fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1; |
| |
} |
| |
/* Not a Bom file */ |
| |
|
| |
i=1; |
| while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { |
while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { |
| linei=linei+1; |
linei=linei+1; |
| for(j=strlen(line); j>=0;j--){ /* Untabifies line */ |
for(j=strlen(line); j>=0;j--){ /* Untabifies line */ |
|
Line 9144 int readdata(char datafile[], int firsto
|
Line 9406 int readdata(char datafile[], int firsto
|
| return 1; |
return 1; |
| } |
} |
| anint[j][i]= (double) year; |
anint[j][i]= (double) year; |
| mint[j][i]= (double)month; |
mint[j][i]= (double)month; |
| |
/* if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ */ |
| |
/* printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */ |
| |
/* fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */ |
| |
/* } */ |
| strcpy(line,stra); |
strcpy(line,stra); |
| } /* End loop on waves */ |
} /* End loop on waves */ |
| |
|
|
Line 9183 int readdata(char datafile[], int firsto
|
Line 9449 int readdata(char datafile[], int firsto
|
| |
|
| } |
} |
| annais[i]=(double)(year); |
annais[i]=(double)(year); |
| moisnais[i]=(double)(month); |
moisnais[i]=(double)(month); |
| |
for (j=1;j<=maxwav;j++){ |
| |
if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ |
| |
printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j,(int)moisnais[i],(int)annais[i]); |
| |
fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j, (int)moisnais[i],(int)annais[i]); |
| |
} |
| |
} |
| |
|
| strcpy(line,stra); |
strcpy(line,stra); |
| |
|
| /* Sample weight */ |
/* Sample weight */ |
|
Line 9641 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 9914 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ |
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ |
| Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ |
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ |
| Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); |
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); |
| for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} |
for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} |
| for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ |
for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ |
| if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ |
if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ |
| Fixed[k]= 0; |
Fixed[k]= 0; |
|
Line 9891 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 10164 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
| /* Searching for doublons in the model */ |
/* Searching for doublons in the model */ |
| for(k1=1; k1<= cptcovt;k1++){ |
for(k1=1; k1<= cptcovt;k1++){ |
| for(k2=1; k2 <k1;k2++){ |
for(k2=1; k2 <k1;k2++){ |
| if((Typevar[k1]==Typevar[k2]) && (Fixed[Tvar[k1]]==Fixed[Tvar[k2]]) && (Dummy[Tvar[k1]]==Dummy[Tvar[k2]] )){ |
/* if((Typevar[k1]==Typevar[k2]) && (Fixed[Tvar[k1]]==Fixed[Tvar[k2]]) && (Dummy[Tvar[k1]]==Dummy[Tvar[k2]] )){ */ |
| |
if((Typevar[k1]==Typevar[k2]) && (Fixed[k1]==Fixed[k2]) && (Dummy[k1]==Dummy[k2] )){ |
| if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */ |
if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */ |
| if(Tvar[k1]==Tvar[k2]){ |
if(Tvar[k1]==Tvar[k2]){ |
| printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); |
printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]); |
| fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog); |
fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]); fflush(ficlog); |
| return(1); |
return(1); |
| } |
} |
| }else if (Typevar[k1] ==2){ |
}else if (Typevar[k1] ==2){ |
|
Line 10067 BOOL IsWow64()
|
Line 10341 BOOL IsWow64()
|
| #endif |
#endif |
| |
|
| void syscompilerinfo(int logged) |
void syscompilerinfo(int logged) |
| { |
{ |
| /* #include "syscompilerinfo.h"*/ |
#include <stdint.h> |
| |
|
| |
/* #include "syscompilerinfo.h"*/ |
| /* command line Intel compiler 32bit windows, XP compatible:*/ |
/* command line Intel compiler 32bit windows, XP compatible:*/ |
| /* /GS /W3 /Gy |
/* /GS /W3 /Gy |
| /Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D |
/Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D |
|
Line 10103 void syscompilerinfo(int logged)
|
Line 10379 void syscompilerinfo(int logged)
|
| /ManifestFile:"Release\IMaCh.exe.intermediate.manifest" /OPT:ICF |
/ManifestFile:"Release\IMaCh.exe.intermediate.manifest" /OPT:ICF |
| /NOLOGO /TLBID:1 |
/NOLOGO /TLBID:1 |
| */ |
*/ |
| |
|
| |
|
| #if defined __INTEL_COMPILER |
#if defined __INTEL_COMPILER |
| #if defined(__GNUC__) |
#if defined(__GNUC__) |
| struct utsname sysInfo; /* For Intel on Linux and OS/X */ |
struct utsname sysInfo; /* For Intel on Linux and OS/X */ |
|
Line 10119 void syscompilerinfo(int logged)
|
Line 10397 void syscompilerinfo(int logged)
|
| } |
} |
| #endif |
#endif |
| |
|
| #include <stdint.h> |
|
| |
|
| printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:"); |
printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:"); |
| #if defined(__clang__) |
#if defined(__clang__) |
| printf(" Clang/LLVM");if(logged)fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ |
printf(" Clang/LLVM");if(logged)fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ |
|
Line 10206 void syscompilerinfo(int logged)
|
Line 10482 void syscompilerinfo(int logged)
|
| #endif |
#endif |
| #endif |
#endif |
| |
|
| // void main() |
// void main () |
| // { |
// { |
| #if defined(_MSC_VER) |
#if defined(_MSC_VER) |
| if (IsWow64()){ |
if (IsWow64()){ |
|
Line 10227 void syscompilerinfo(int logged)
|
Line 10503 void syscompilerinfo(int logged)
|
| } |
} |
| |
|
| int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){ |
int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){ |
| /*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Prevalence limit (forward period or forward stable prevalence) --------------*/ |
| int i, j, k, i1, k4=0, nres=0 ; |
int i, j, k, i1, k4=0, nres=0 ; |
| /* double ftolpl = 1.e-10; */ |
/* double ftolpl = 1.e-10; */ |
| double age, agebase, agelim; |
double age, agebase, agelim; |
|
Line 10236 int prevalence_limit(double *p, double *
|
Line 10512 int prevalence_limit(double *p, double *
|
| strcpy(filerespl,"PL_"); |
strcpy(filerespl,"PL_"); |
| strcat(filerespl,fileresu); |
strcat(filerespl,fileresu); |
| if((ficrespl=fopen(filerespl,"w"))==NULL) { |
if((ficrespl=fopen(filerespl,"w"))==NULL) { |
| printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; |
printf("Problem with forward period (stable) prevalence resultfile: %s\n", filerespl);return 1; |
| fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; |
fprintf(ficlog,"Problem with forward period (stable) prevalence resultfile: %s\n", filerespl);return 1; |
| } |
} |
| printf("\nComputing period (stable) prevalence: result on file '%s' \n", filerespl); |
printf("\nComputing forward period (stable) prevalence: result on file '%s' \n", filerespl); |
| fprintf(ficlog,"\nComputing period (stable) prevalence: result on file '%s' \n", filerespl); |
fprintf(ficlog,"\nComputing forward period (stable) prevalence: result on file '%s' \n", filerespl); |
| pstamp(ficrespl); |
pstamp(ficrespl); |
| fprintf(ficrespl,"# Period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl); |
fprintf(ficrespl,"# Forward period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl); |
| fprintf(ficrespl,"#Age "); |
fprintf(ficrespl,"#Age "); |
| for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
| fprintf(ficrespl,"\n"); |
fprintf(ficrespl,"\n"); |
|
Line 10317 int prevalence_limit(double *p, double *
|
Line 10593 int prevalence_limit(double *p, double *
|
| } |
} |
| |
|
| int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj){ |
int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj){ |
| /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Back Prevalence limit (backward stable prevalence) --------------*/ |
| |
|
| /* Computes the back prevalence limit for any combination of covariate values |
/* Computes the back prevalence limit for any combination of covariate values |
| * at any age between ageminpar and agemaxpar |
* at any age between ageminpar and agemaxpar |
|
Line 10332 int back_prevalence_limit(double *p, dou
|
Line 10608 int back_prevalence_limit(double *p, dou
|
| strcpy(fileresplb,"PLB_"); |
strcpy(fileresplb,"PLB_"); |
| strcat(fileresplb,fileresu); |
strcat(fileresplb,fileresu); |
| if((ficresplb=fopen(fileresplb,"w"))==NULL) { |
if((ficresplb=fopen(fileresplb,"w"))==NULL) { |
| printf("Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1; |
printf("Problem with backward prevalence resultfile: %s\n", fileresplb);return 1; |
| fprintf(ficlog,"Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1; |
fprintf(ficlog,"Problem with backward prevalence resultfile: %s\n", fileresplb);return 1; |
| } |
} |
| printf("Computing period (stable) back prevalence: result on file '%s' \n", fileresplb); |
printf("Computing backward prevalence: result on file '%s' \n", fileresplb); |
| fprintf(ficlog,"Computing period (stable) back prevalence: result on file '%s' \n", fileresplb); |
fprintf(ficlog,"Computing backward prevalence: result on file '%s' \n", fileresplb); |
| pstamp(ficresplb); |
pstamp(ficresplb); |
| fprintf(ficresplb,"# Period (stable) back prevalence. Precision given by ftolpl=%g \n", ftolpl); |
fprintf(ficresplb,"# Backward prevalence. Precision given by ftolpl=%g \n", ftolpl); |
| fprintf(ficresplb,"#Age "); |
fprintf(ficresplb,"#Age "); |
| for(i=1; i<=nlstate;i++) fprintf(ficresplb,"%d-%d ",i,i); |
for(i=1; i<=nlstate;i++) fprintf(ficresplb,"%d-%d ",i,i); |
| fprintf(ficresplb,"\n"); |
fprintf(ficresplb,"\n"); |
|
Line 10527 int hPijx(double *p, int bage, int fage)
|
Line 10803 int hPijx(double *p, int bage, int fage)
|
| /*if (stepm<=24) stepsize=2;*/ |
/*if (stepm<=24) stepsize=2;*/ |
| |
|
| /* agelim=AGESUP; */ |
/* agelim=AGESUP; */ |
| ageminl=30; |
ageminl=AGEINF; /* was 30 */ |
| hstepm=stepsize*YEARM; /* Every year of age */ |
hstepm=stepsize*YEARM; /* Every year of age */ |
| hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ |
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ |
| |
|
|
Line 10557 int hPijx(double *p, int bage, int fage)
|
Line 10833 int hPijx(double *p, int bage, int fage)
|
| /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */ |
/* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */ |
| for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */ |
for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */ |
| /* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */ |
/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */ |
| nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ |
nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm+0.1)-1; /* Typically 20 years = 20*12/6=40 or 55*12/24=27.5-1.1=>27 */ |
| nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */ |
nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 or 28*/ |
| |
|
| /* nhstepm=nhstepm*YEARM; aff par mois*/ |
/* nhstepm=nhstepm*YEARM; aff par mois*/ |
| |
|
|
Line 10606 int main(int argc, char *argv[])
|
Line 10882 int main(int argc, char *argv[])
|
| double ssval; |
double ssval; |
| #endif |
#endif |
| 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 i,j, k, n=MAXN,iter=0,m,size=100, cptcod; |
int i,j, k, iter=0,m,size=100, cptcod; /* Suppressing because nobs */ |
| |
/* int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; */ |
| int ncvyear=0; /* Number of years needed for the period prevalence to converge */ |
int ncvyear=0; /* Number of years needed for the period prevalence to converge */ |
| int jj, ll, li, lj, lk; |
int jj, ll, li, lj, lk; |
| int numlinepar=0; /* Current linenumber of parameter file */ |
int numlinepar=0; /* Current linenumber of parameter file */ |
|
Line 10641 int main(int argc, char *argv[])
|
Line 10918 int main(int argc, char *argv[])
|
| |
|
| char pathr[MAXLINE], pathimach[MAXLINE]; |
char pathr[MAXLINE], pathimach[MAXLINE]; |
| char *tok, *val; /* pathtot */ |
char *tok, *val; /* pathtot */ |
| int firstobs=1, lastobs=10; |
int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs declared globally ;*/ |
| int c, h , cpt, c2; |
int c, h , cpt, c2; |
| int jl=0; |
int jl=0; |
| int i1, j1, jk, stepsize=0; |
int i1, j1, jk, stepsize=0; |
|
Line 10649 int main(int argc, char *argv[])
|
Line 10926 int main(int argc, char *argv[])
|
| |
|
| int *tab; |
int *tab; |
| int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ |
int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ |
| int backcast=0; |
/* double anprojd, mprojd, jprojd; /\* For eventual projections *\/ */ |
| |
/* double anprojf, mprojf, jprojf; */ |
| |
/* double jintmean,mintmean,aintmean; */ |
| |
int prvforecast = 0; /* Might be 1 (date of beginning of projection is a choice or 2 is the dateintmean */ |
| |
int prvbackcast = 0; /* Might be 1 (date of beginning of projection is a choice or 2 is the dateintmean */ |
| |
double yrfproj= 10.0; /* Number of years of forward projections */ |
| |
double yrbproj= 10.0; /* Number of years of backward projections */ |
| |
int prevbcast=0; /* defined as global for mlikeli and mle, replacing backcast */ |
| int mobilav=0,popforecast=0; |
int mobilav=0,popforecast=0; |
| int hstepm=0, nhstepm=0; |
int hstepm=0, nhstepm=0; |
| int agemortsup; |
int agemortsup; |
|
Line 10674 int main(int argc, char *argv[])
|
Line 10958 int main(int argc, char *argv[])
|
| double *epj, vepp; |
double *epj, vepp; |
| |
|
| double dateprev1, dateprev2; |
double dateprev1, dateprev2; |
| double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0; |
double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0, dateprojd=0, dateprojf=0; |
| double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0; |
double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0, datebackd=0, datebackf=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 10753 int main(int argc, char *argv[])
|
Line 11038 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 10832 int main(int argc, char *argv[])
|
Line 11122 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 10842 int main(int argc, char *argv[])
|
Line 11130 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? */ |
/* Is it a BOM UTF-8 Windows file? */ |
|
Line 10854 int main(int argc, char *argv[])
|
Line 11155 int main(int argc, char *argv[])
|
| noffset=noffset+3; |
noffset=noffset+3; |
| printf("# File is an UTF8 Bom.\n"); // 0xBF |
printf("# File is an UTF8 Bom.\n"); // 0xBF |
| } |
} |
| else if( line[0] == (char)0xFE && line[1] == (char)0xFF) |
/* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/ |
| |
else if( line[0] == (char)0xFF && line[1] == (char)0xFE) |
| { |
{ |
| noffset=noffset+2; |
noffset=noffset+2; |
| printf("# File is an UTF16BE BOM file\n"); |
printf("# File is an UTF16BE BOM file\n"); |
|
Line 10874 int main(int argc, char *argv[])
|
Line 11176 int main(int argc, char *argv[])
|
| numlinepar++; |
numlinepar++; |
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| |
fputs(line,ficres); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| continue; |
continue; |
| }else |
}else |
|
Line 10883 int main(int argc, char *argv[])
|
Line 11186 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 10904 int main(int argc, char *argv[])
|
Line 11213 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); |
if( lastpass > maxwav){ |
| |
printf("Error (lastpass = %d) > (maxwav = %d)\n",lastpass, maxwav); |
| |
fprintf(ficlog,"Error (lastpass = %d) > (maxwav = %d)\n",lastpass, maxwav); |
| |
fflush(ficlog); |
| |
goto end; |
| |
} |
| |
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, 0, 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 10914 int main(int argc, char *argv[])
|
Line 11234 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+' instead of '%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+' instead of '%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 10942 int main(int argc, char *argv[])
|
Line 11258 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 10975 int main(int argc, char *argv[])
|
Line 11294 int main(int argc, char *argv[])
|
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| |
|
| covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ |
covar=matrix(0,NCOVMAX,firstobs,lastobs); /**< used in readdata */ |
| if(nqv>=1)coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */ |
if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs); /**< Fixed quantitative covariate */ |
| if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */ |
if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs); /**< Time varying quantitative covariate */ |
| if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/ |
if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs); /**< Time varying covariate (dummy and quantitative)*/ |
| cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/ |
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/ |
| /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5 |
/* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5 |
| v1+v2*age+v2*v3 makes cptcovn = 3 |
v1+v2*age+v2*v3 makes cptcovn = 3 |
|
Line 11041 int main(int argc, char *argv[])
|
Line 11360 int main(int argc, char *argv[])
|
| for(jj=1; jj <=nlstate+ndeath; jj++){ |
for(jj=1; jj <=nlstate+ndeath; jj++){ |
| if(jj==i) continue; |
if(jj==i) continue; |
| j++; |
j++; |
| |
while((c=getc(ficpar))=='#' && c!= EOF){ |
| |
ungetc(c,ficpar); |
| |
fgets(line, MAXLINE, ficpar); |
| |
numlinepar++; |
| |
fputs(line,stdout); |
| |
fputs(line,ficparo); |
| |
fputs(line,ficlog); |
| |
} |
| |
ungetc(c,ficpar); |
| fscanf(ficpar,"%1d%1d",&i1,&j1); |
fscanf(ficpar,"%1d%1d",&i1,&j1); |
| if ((i1 != i) || (j1 != jj)){ |
if ((i1 != i) || (j1 != jj)){ |
| printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \ |
printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \ |
|
Line 11177 Please run with mle=-1 to get a correct
|
Line 11505 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 |
| */ |
*/ |
| n= lastobs; |
nobs=lastobs-firstobs+1; /* was = lastobs;*/ |
| num=lvector(1,n); |
/* num=lvector(1,n); */ |
| moisnais=vector(1,n); |
/* moisnais=vector(1,n); */ |
| annais=vector(1,n); |
/* annais=vector(1,n); */ |
| moisdc=vector(1,n); |
/* moisdc=vector(1,n); */ |
| andc=vector(1,n); |
/* andc=vector(1,n); */ |
| weight=vector(1,n); |
/* weight=vector(1,n); */ |
| agedc=vector(1,n); |
/* agedc=vector(1,n); */ |
| cod=ivector(1,n); |
/* cod=ivector(1,n); */ |
| for(i=1;i<=n;i++){ |
/* for(i=1;i<=n;i++){ */ |
| |
num=lvector(firstobs,lastobs); |
| |
moisnais=vector(firstobs,lastobs); |
| |
annais=vector(firstobs,lastobs); |
| |
moisdc=vector(firstobs,lastobs); |
| |
andc=vector(firstobs,lastobs); |
| |
weight=vector(firstobs,lastobs); |
| |
agedc=vector(firstobs,lastobs); |
| |
cod=ivector(firstobs,lastobs); |
| |
for(i=firstobs;i<=lastobs;i++){ |
| num[i]=0; |
num[i]=0; |
| moisnais[i]=0; |
moisnais[i]=0; |
| annais[i]=0; |
annais[i]=0; |
|
Line 11210 Please run with mle=-1 to get a correct
|
Line 11537 Please run with mle=-1 to get a correct
|
| cod[i]=0; |
cod[i]=0; |
| weight[i]=1.0; /* Equal weights, 1 by default */ |
weight[i]=1.0; /* Equal weights, 1 by default */ |
| } |
} |
| mint=matrix(1,maxwav,1,n); |
mint=matrix(1,maxwav,firstobs,lastobs); |
| anint=matrix(1,maxwav,1,n); |
anint=matrix(1,maxwav,firstobs,lastobs); |
| s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ |
s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ |
| tab=ivector(1,NCOVMAX); |
tab=ivector(1,NCOVMAX); |
| ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
| ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
|
Line 11314 Please run with mle=-1 to get a correct
|
Line 11641 Please run with mle=-1 to get a correct
|
| |
|
| |
|
| agegomp=(int)agemin; |
agegomp=(int)agemin; |
| free_vector(moisnais,1,n); |
free_vector(moisnais,firstobs,lastobs); |
| free_vector(annais,1,n); |
free_vector(annais,firstobs,lastobs); |
| /* free_matrix(mint,1,maxwav,1,n); |
/* free_matrix(mint,1,maxwav,1,n); |
| free_matrix(anint,1,maxwav,1,n);*/ |
free_matrix(anint,1,maxwav,1,n);*/ |
| /* free_vector(moisdc,1,n); */ |
/* free_vector(moisdc,1,n); */ |
|
Line 11341 Please run with mle=-1 to get a correct
|
Line 11668 Please run with mle=-1 to get a correct
|
| concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); |
concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); |
| /* Concatenates waves */ |
/* Concatenates waves */ |
| |
|
| free_vector(moisdc,1,n); |
free_vector(moisdc,firstobs,lastobs); |
| free_vector(andc,1,n); |
free_vector(andc,firstobs,lastobs); |
| |
|
| /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ |
/* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ |
| nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
|
Line 11524 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 11851 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,"<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",\ |
fprintf(fichtm,"<h4>Parameter line 2</h4><ul><li>Tolerance for the convergence of the likelihood: ftol=%g \n<li>Interval for the elementary matrix (in month): stepm=%d",\ |
| ftol, stepm); |
ftol, stepm); |
| fprintf(fichtm,"\n<li>Number of fixed dummy covariates: ncovcol=%d ", ncovcol); |
fprintf(fichtm,"\n<li>Number of fixed dummy covariates: ncovcol=%d ", ncovcol); |
| ncurrv=1; |
ncurrv=1; |
|
Line 11532 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 11859 Title=%s <br>Datafile=%s Firstpass=%d La
|
| fprintf(fichtm,"\n<li> Number of fixed quantitative variables: nqv=%d ", nqv); |
fprintf(fichtm,"\n<li> Number of fixed quantitative variables: nqv=%d ", nqv); |
| ncurrv=i; |
ncurrv=i; |
| for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", 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); |
fprintf(fichtm,"\n<li> Number of time varying (wave varying) dummy covariates: ntv=%d ", ntv); |
| ncurrv=i; |
ncurrv=i; |
| for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", 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); |
fprintf(fichtm,"\n<li>Number of time varying quantitative covariates: nqtv=%d ", nqtv); |
| ncurrv=i; |
ncurrv=i; |
| for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", 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", \ |
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", \ |
|
Line 11567 Interval (in months) between two waves:
|
Line 11894 Interval (in months) between two waves:
|
| for(j=1;j<=NDIM;j++) |
for(j=1;j<=NDIM;j++) |
| ximort[i][j]=0.; |
ximort[i][j]=0.; |
| /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ |
/* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ |
| cens=ivector(1,n); |
cens=ivector(firstobs,lastobs); |
| ageexmed=vector(1,n); |
ageexmed=vector(firstobs,lastobs); |
| agecens=vector(1,n); |
agecens=vector(firstobs,lastobs); |
| dcwave=ivector(1,n); |
dcwave=ivector(firstobs,lastobs); |
| |
|
| for (i=1; i<=imx; i++){ |
for (i=1; i<=imx; i++){ |
| dcwave[i]=-1; |
dcwave[i]=-1; |
|
Line 11604 Interval (in months) between two waves:
|
Line 11931 Interval (in months) between two waves:
|
| ximort[i][j]=(i == j ? 1.0 : 0.0); |
ximort[i][j]=(i == j ? 1.0 : 0.0); |
| } |
} |
| |
|
| /*p[1]=0.0268; p[NDIM]=0.083;*/ |
p[1]=0.0268; p[NDIM]=0.083; |
| /*printf("%lf %lf", p[1], p[2]);*/ |
/* printf("%lf %lf", p[1], p[2]); */ |
| |
|
| |
|
| #ifdef GSL |
#ifdef GSL |
|
Line 11731 Interval (in months) between two waves:
|
Line 12058 Interval (in months) between two waves:
|
| printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
| fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
| } |
} |
| lsurv=vector(1,AGESUP); |
lsurv=vector(agegomp,AGESUP); |
| lpop=vector(1,AGESUP); |
lpop=vector(agegomp,AGESUP); |
| tpop=vector(1,AGESUP); |
tpop=vector(agegomp,AGESUP); |
| lsurv[agegomp]=100000; |
lsurv[agegomp]=100000; |
| |
|
| for (k=agegomp;k<=AGESUP;k++) { |
for (k=agegomp;k<=AGESUP;k++) { |
|
Line 11780 Please run with mle=-1 to get a correct
|
Line 12107 Please run with mle=-1 to get a correct
|
| stepm, weightopt,\ |
stepm, weightopt,\ |
| model,imx,p,matcov,agemortsup); |
model,imx,p,matcov,agemortsup); |
| |
|
| free_vector(lsurv,1,AGESUP); |
free_vector(lsurv,agegomp,AGESUP); |
| free_vector(lpop,1,AGESUP); |
free_vector(lpop,agegomp,AGESUP); |
| free_vector(tpop,1,AGESUP); |
free_vector(tpop,agegomp,AGESUP); |
| free_matrix(ximort,1,NDIM,1,NDIM); |
free_matrix(ximort,1,NDIM,1,NDIM); |
| free_ivector(cens,1,n); |
free_ivector(dcwave,firstobs,lastobs); |
| free_vector(agecens,1,n); |
free_vector(agecens,firstobs,lastobs); |
| free_ivector(dcwave,1,n); |
free_vector(ageexmed,firstobs,lastobs); |
| |
free_ivector(cens,firstobs,lastobs); |
| #ifdef GSL |
#ifdef GSL |
| #endif |
#endif |
| } /* Endof if mle==-3 mortality only */ |
} /* Endof if mle==-3 mortality only */ |
|
Line 11820 Please run with mle=-1 to get a correct
|
Line 12148 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 11982 Please run with mle=-1 to get a correct
|
Line 12310 Please run with mle=-1 to get a correct
|
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| |
fputs(line,ficres); |
| continue; |
continue; |
| }else |
}else |
| break; |
break; |
|
Line 12027 Please run with mle=-1 to get a correct
|
Line 12356 Please run with mle=-1 to get a correct
|
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| |
fputs(line,ficres); |
| continue; |
continue; |
| }else |
}else |
| break; |
break; |
|
Line 12052 Please run with mle=-1 to get a correct
|
Line 12382 Please run with mle=-1 to get a correct
|
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| |
fputs(line,ficres); |
| continue; |
continue; |
| }else |
}else |
| break; |
break; |
|
Line 12085 Please run with mle=-1 to get a correct
|
Line 12416 Please run with mle=-1 to get a correct
|
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| fputs(line,ficlog); |
fputs(line,ficlog); |
| |
fputs(line,ficres); |
| continue; |
continue; |
| }else if(sscanf(line,"prevforecast=%[^\n]\n",modeltemp)) |
}else if(sscanf(line,"prevforecast=%[^\n]\n",modeltemp)) |
| parameterline=11; |
parameterline=11; |
| else if(sscanf(line,"backcast=%[^\n]\n",modeltemp)) |
else if(sscanf(line,"prevbackcast=%[^\n]\n",modeltemp)) |
| parameterline=12; |
parameterline=12; |
| else if(sscanf(line,"result:%[^\n]\n",modeltemp)) |
else if(sscanf(line,"result:%[^\n]\n",modeltemp)) |
| parameterline=13; |
parameterline=13; |
|
Line 12097 Please run with mle=-1 to get a correct
|
Line 12429 Please run with mle=-1 to get a correct
|
| } |
} |
| switch (parameterline){ |
switch (parameterline){ |
| case 11: |
case 11: |
| if((num_filled=sscanf(line,"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)) !=EOF){ |
if((num_filled=sscanf(line,"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)) !=EOF && (num_filled == 8)){ |
| if (num_filled != 8) { |
fprintf(ficparo,"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); |
| printf("Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line); |
|
| fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line); |
|
| goto end; |
|
| } |
|
| fprintf(ficparo,"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); |
|
| printf("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); |
printf("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(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.; |
dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.; |
| dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.; |
dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.; |
| |
prvforecast = 1; |
| |
} |
| |
else if((num_filled=sscanf(line,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",&prevfcast,&yrfproj,&mobilavproj)) !=EOF){/* && (num_filled == 3))*/ |
| |
printf("prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
| |
fprintf(ficlog,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
| |
fprintf(ficres,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
| |
prvforecast = 2; |
| |
} |
| |
else { |
| |
printf("Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevforecast=1 yearsfproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line); |
| |
fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevforecast=1 yearproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line); |
| |
goto end; |
| } |
} |
| break; |
break; |
| case 12: |
case 12: |
| /*fscanf(ficpar,"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);*/ |
if((num_filled=sscanf(line,"prevbackcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&prevbcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj)) !=EOF && (num_filled == 8)){ |
| if((num_filled=sscanf(line,"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)) !=EOF){ |
fprintf(ficparo,"prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| if (num_filled != 8) { |
printf("prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| printf("Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 final-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line); |
fprintf(ficlog,"prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 final-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line); |
fprintf(ficres,"prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| goto end; |
/* day and month of back2 are not used but only year anback2.*/ |
| } |
|
| printf("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(ficparo,"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); |
|
| /* day and month of proj2 are not used but only year anproj2.*/ |
|
| dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.; |
dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.; |
| dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.; |
dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.; |
| |
prvbackcast = 1; |
| |
} |
| |
else if((num_filled=sscanf(line,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",&prevbcast,&yrbproj,&mobilavproj)) ==3){/* && (num_filled == 3))*/ |
| |
printf("prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
| |
fprintf(ficlog,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
| |
fprintf(ficres,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
| |
prvbackcast = 2; |
| |
} |
| |
else { |
| |
printf("Error: Not 8 (data)parameters in line but %d, for example:prevbackcast=1 starting-back-date=1/1/1990 final-back-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevbackcast=1 yearsbproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line); |
| |
fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevbackcast=1 starting-back-date=1/1/1990 final-back-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevbackcast=1 yearbproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line); |
| |
goto end; |
| } |
} |
| break; |
break; |
| case 13: |
case 13: |
|
Line 12153 Please run with mle=-1 to get a correct
|
Line 12497 Please run with mle=-1 to get a correct
|
| fprintf(ficres,"result: %s\n",resultline); |
fprintf(ficres,"result: %s\n",resultline); |
| fprintf(ficlog,"result: %s\n",resultline); |
fprintf(ficlog,"result: %s\n",resultline); |
| break; |
break; |
| case 14: |
case 14: |
| if(ncovmodel >2 && nresult==0 ){ |
printf("Error: Unknown command '%s'\n",line); |
| |
fprintf(ficlog,"Error: Unknown command '%s'\n",line); |
| |
if(ncovmodel >=2 && nresult==0 ){ |
| printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); |
printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); |
| goto end; |
fprintf(ficlog,"ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); |
| } |
} |
| |
goto end; |
| break; |
break; |
| default: |
default: |
| nresult=1; |
nresult=1; |
|
Line 12179 This is probably because your parameter
|
Line 12526 This is probably because your parameter
|
| 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); |
/* It seems that anprojd which is computed from the mean year at interview which is known yet because of freqsummary */ |
| |
/* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */ /* Done in freqsummary */ |
| |
if(prvforecast==1){ |
| |
dateprojd=(jproj1+12*mproj1+365*anproj1)/365; |
| |
jprojd=jproj1; |
| |
mprojd=mproj1; |
| |
anprojd=anproj1; |
| |
dateprojf=(jproj2+12*mproj2+365*anproj2)/365; |
| |
jprojf=jproj2; |
| |
mprojf=mproj2; |
| |
anprojf=anproj2; |
| |
} else if(prvforecast == 2){ |
| |
dateprojd=dateintmean; |
| |
date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); |
| |
dateprojf=dateintmean+yrfproj; |
| |
date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); |
| |
} |
| |
if(prvbackcast==1){ |
| |
datebackd=(jback1+12*mback1+365*anback1)/365; |
| |
jbackd=jback1; |
| |
mbackd=mback1; |
| |
anbackd=anback1; |
| |
datebackf=(jback2+12*mback2+365*anback2)/365; |
| |
jbackf=jback2; |
| |
mbackf=mback2; |
| |
anbackf=anback2; |
| |
} else if(prvbackcast == 2){ |
| |
datebackd=dateintmean; |
| |
date2dmy(datebackd,&jbackd, &mbackd, &anbackd); |
| |
datebackf=dateintmean-yrbproj; |
| |
date2dmy(datebackf,&jbackf, &mbackf, &anbackf); |
| |
} |
| |
|
| |
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, prevbcast, pathc,p, (int)anprojd-bage, (int)anbackd-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,prevbcast, estepm, \ |
| jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2); |
jprev1,mprev1,anprev1,dateprev1, dateprojd, datebackd,jprev2,mprev2,anprev2,dateprev2,dateprojf, datebackf); |
| |
|
| /*------------ free_vector -------------*/ |
/*------------ free_vector -------------*/ |
| /* chdir(path); */ |
/* chdir(path); */ |
|
Line 12192 Please run with mle=-1 to get a correct
|
Line 12572 Please run with mle=-1 to get a correct
|
| /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ |
/* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ |
| /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ |
/* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ |
| /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ |
/* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ |
| free_lvector(num,1,n); |
free_lvector(num,firstobs,lastobs); |
| free_vector(agedc,1,n); |
free_vector(agedc,firstobs,lastobs); |
| /*free_matrix(covar,0,NCOVMAX,1,n);*/ |
/*free_matrix(covar,0,NCOVMAX,1,n);*/ |
| /*free_matrix(covar,1,NCOVMAX,1,n);*/ |
/*free_matrix(covar,1,NCOVMAX,1,n);*/ |
| fclose(ficparo); |
fclose(ficparo); |
|
Line 12256 Please run with mle=-1 to get a correct
|
Line 12636 Please run with mle=-1 to get a correct
|
| }/* end if moving average */ |
}/* end if moving average */ |
| |
|
| /*---------- Forecasting ------------------*/ |
/*---------- Forecasting ------------------*/ |
| if(prevfcast==1){ |
if(prevfcast==1){ |
| /* if(stepm ==1){*/ |
/* /\* if(stepm ==1){*\/ */ |
| prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); |
/* /\* anproj1, mproj1, jproj1 either read explicitly or yrfproj *\/ */ |
| |
/*This done previously after freqsummary.*/ |
| |
/* dateprojd=(jproj1+12*mproj1+365*anproj1)/365; */ |
| |
/* dateprojf=(jproj2+12*mproj2+365*anproj2)/365; */ |
| |
|
| |
/* } else if (prvforecast==2){ */ |
| |
/* /\* if(stepm ==1){*\/ */ |
| |
/* /\* anproj1, mproj1, jproj1 either read explicitly or yrfproj *\/ */ |
| |
/* } */ |
| |
/*prevforecast(fileresu, dateintmean, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);*/ |
| |
prevforecast(fileresu,dateintmean, dateprojd, dateprojf, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, p, cptcoveff); |
| } |
} |
| |
|
| /* Backcasting */ |
/* Prevbcasting */ |
| if(backcast==1){ |
if(prevbcast==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); |
| ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
|
Line 12277 Please run with mle=-1 to get a correct
|
Line 12667 Please run with mle=-1 to get a correct
|
| hBijx(p, bage, fage, mobaverage); |
hBijx(p, bage, fage, mobaverage); |
| fclose(ficrespijb); |
fclose(ficrespijb); |
| |
|
| prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, |
/* /\* prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, *\/ */ |
| mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); |
/* /\* mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); *\/ */ |
| |
/* prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, */ |
| |
/* mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */ |
| |
prevbackforecast(fileresu, mobaverage, dateintmean, dateprojd, dateprojf, agemin, agemax, dateprev1, dateprev2, |
| |
mobilavproj, bage, fage, firstpass, lastpass, p, cptcoveff); |
| |
|
| |
|
| varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); |
varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); |
| |
|
| |
|
|
Line 12286 Please run with mle=-1 to get a correct
|
Line 12682 Please run with mle=-1 to get a correct
|
| 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 */ |
} /* end Prevbcasting */ |
| |
|
| |
|
| /* ------ Other prevalence ratios------------ */ |
/* ------ Other prevalence ratios------------ */ |
|
Line 12450 Please run with mle=-1 to get a correct
|
Line 12846 Please run with mle=-1 to get a correct
|
| if(vpopbased==1) |
if(vpopbased==1) |
| fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); |
fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); |
| else |
else |
| fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); |
fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n"); |
| fprintf(ficrest,"# Age popbased mobilav e.. (std) "); |
fprintf(ficrest,"# Age popbased mobilav e.. (std) "); |
| 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"); */ |
| printf("Computing age specific period (stable) prevalences in each health state \n"); |
printf("Computing age specific forward 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 forward period (stable) prevalences in each health state \n"); |
| for(age=bage; age <=fage ;age++){ |
for(age=bage; age <=fage ;age++){ |
| prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k, nres); /*ZZ Is it the correct prevalim */ |
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k, nres); /*ZZ Is it the correct prevalim */ |
| if (vpopbased==1) { |
if (vpopbased==1) { |
|
Line 12503 Please run with mle=-1 to get a correct
|
Line 12899 Please run with mle=-1 to get a correct
|
| 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-covariance of period prevalence*/ |
/* variance-covariance of forward period prevalence*/ |
| varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); |
varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); |
| |
|
| |
|
| free_vector(weight,1,n); |
free_vector(weight,firstobs,lastobs); |
| free_imatrix(Tvard,1,NCOVMAX,1,2); |
free_imatrix(Tvard,1,NCOVMAX,1,2); |
| free_imatrix(s,1,maxwav+1,1,n); |
free_imatrix(s,1,maxwav+1,firstobs,lastobs); |
| free_matrix(anint,1,maxwav,1,n); |
free_matrix(anint,1,maxwav,firstobs,lastobs); |
| free_matrix(mint,1,maxwav,1,n); |
free_matrix(mint,1,maxwav,firstobs,lastobs); |
| free_ivector(cod,1,n); |
free_ivector(cod,firstobs,lastobs); |
| free_ivector(tab,1,NCOVMAX); |
free_ivector(tab,1,NCOVMAX); |
| fclose(ficresstdeij); |
fclose(ficresstdeij); |
| fclose(ficrescveij); |
fclose(ficrescveij); |
|
Line 12532 Please run with mle=-1 to get a correct
|
Line 12928 Please run with mle=-1 to get a correct
|
| free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n); |
if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs); |
| if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); |
if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs); |
| if(nqv>=1)free_matrix(coqvar,1,nqv,1,n); |
if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs); |
| free_matrix(covar,0,NCOVMAX,1,n); |
free_matrix(covar,0,NCOVMAX,firstobs,lastobs); |
| free_matrix(matcov,1,npar,1,npar); |
free_matrix(matcov,1,npar,1,npar); |
| free_matrix(hess,1,npar,1,npar); |
free_matrix(hess,1,npar,1,npar); |
| /*free_vector(delti,1,npar);*/ |
/*free_vector(delti,1,npar);*/ |
|
Line 12618 Please run with mle=-1 to get a correct
|
Line 13014 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 12653 Please run with mle=-1 to get a correct
|
Line 13051 Please run with mle=-1 to get a correct
|
| |
|
| sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); |
sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); |
| printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); |
printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); |
| |
strcpy(pplotcmd,plotcmd); |
| |
|
| if((outcmd=system(plotcmd)) != 0){ |
if((outcmd=system(plotcmd)) != 0){ |
| printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); |
printf("Error in gnuplot, command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); |
| printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); |
printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); |
| sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); |
sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); |
| if((outcmd=system(plotcmd)) != 0) |
if((outcmd=system(plotcmd)) != 0){ |
| printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); |
printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); |
| |
strcpy(plotcmd,pplotcmd); |
| |
} |
| } |
} |
| printf(" Successful, please wait..."); |
printf(" Successful, please wait..."); |
| while (z[0] != 'q') { |
while (z[0] != 'q') { |
|
Line 12686 end:
|
Line 13087 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