From e412aaf5070046a41e178e334d87cfaa0c018da0 Mon Sep 17 00:00:00 2001
From: "N. Brouard" <brouard@ined.fr>
Date: Fri, 22 Feb 2002 17:54:20 +0000
Subject: [PATCH] =?utf8?q?First=20merge=20between=20Agn=EF=BF=BD=EF=BF=BDs?=
=?utf8?q?=20and=20Nicolas!!?=
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit
---
src/imach.c | 168 +++++++++++++++++++++++++++++++++-------------------
1 file changed, 106 insertions(+), 62 deletions(-)
diff --git a/src/imach.c b/src/imach.c
index 5845ea1..94e55d4 100644
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,34 +1,42 @@
-
-/*********************** Imach **************************************
- This program computes Healthy Life Expectancies from cross-longitudinal
- data. Cross-longitudinal consist in a first survey ("cross") where
- individuals from different ages are interviewed on their health status
- or degree of disability. At least a second wave of interviews
- ("longitudinal") should measure each new individual health status.
- Health expectancies are computed from the transistions observed between
- waves and are computed for each degree of severity of disability (number
- of life states). More degrees you consider, more time is necessary to
- reach the Maximum Likelihood of the parameters involved in the model.
- The simplest model is the multinomial logistic model where pij is
- the probabibility to be observed in state j at the second wave conditional
- to be observed in state i at the first wave. Therefore the model is:
- log(pij/pii)= aij + bij*age+ cij*sex + etc , where 'age' is age and 'sex'
- is a covariate. If you want to have a more complex model than "constant and
- age", you should modify the program where the markup
- *Covariates have to be included here again* invites you to do it.
- More covariates you add, less is the speed of the convergence.
-
- The advantage that this computer programme claims, comes from that if the
- delay between waves is not identical for each individual, or if some
- individual missed an interview, the information is not rounded or lost, but
- taken into account using an interpolation or extrapolation.
- hPijx is the probability to be
- observed in state i at age x+h conditional to the observed state i at age
- x. The delay 'h' can be split into an exact number (nh*stepm) of
- unobserved intermediate states. This elementary transition (by month or
- quarter trimester, semester or year) is model as a multinomial logistic.
- The hPx matrix is simply the matrix product of nh*stepm elementary matrices
- and the contribution of each individual to the likelihood is simply hPijx.
+/* $Id$
+ Interpolate Markov Chain
+
+ Short summary of the programme:
+
+ This program computes Healthy Life Expectancies from
+ cross-longitudinal data. Cross-longitudinal data consist in: -1- a
+ first survey ("cross") where individuals from different ages are
+ interviewed on their health status or degree of disability (in the
+ case of a health survey which is our main interest) -2- at least a
+ second wave of interviews ("longitudinal") which measure each change
+ (if any) in individual health status. Health expectancies are
+ computed from the time spent in each health state according to a
+ model. More health states you consider, more time is necessary to reach the
+ Maximum Likelihood of the parameters involved in the model. The
+ simplest model is the multinomial logistic model where pij is the
+ probabibility to be observed in state j at the second wave
+ conditional to be observed in state i at the first wave. Therefore
+ the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
+ 'age' is age and 'sex' is a covariate. If you want to have a more
+ complex model than "constant and age", you should modify the program
+ where the markup *Covariates have to be included here again* invites
+ you to do it. More covariates you add, slower the
+ convergence.
+
+ The advantage of this computer programme, compared to a simple
+ multinomial logistic model, is clear when the delay between waves is not
+ identical for each individual. Also, if a individual missed an
+ intermediate interview, the information is lost, but taken into
+ account using an interpolation or extrapolation.
+
+ hPijx is the probability to be observed in state i at age x+h
+ conditional to the observed state i at age x. The delay 'h' can be
+ split into an exact number (nh*stepm) of unobserved intermediate
+ states. This elementary transition (by month or quarter trimester,
+ semester or year) is model as a multinomial logistic. The hPx
+ matrix is simply the matrix product of nh*stepm elementary matrices
+ and the contribution of each individual to the likelihood is simply
+ hPijx.
Also this programme outputs the covariance matrix of the parameters but also
of the life expectancies. It also computes the prevalence limits.
@@ -48,6 +56,7 @@
#include <unistd.h>
#define MAXLINE 256
+#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"
#define FILENAMELENGTH 80
/*#define DEBUG*/
#define windows
@@ -141,14 +150,18 @@ double ftol=FTOL; /* Tolerance for computing Max Likelihood */
double ftolhess; /* Tolerance for computing hessian */
/**************** split *************************/
-static int split( char *path, char *dirc, char *name )
+static int split( char *path, char *dirc, char *name, char *ext, char *finame )
{
char *s; /* pointer */
int l1, l2; /* length counters */
l1 = strlen( path ); /* length of path */
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
+#ifdef windows
s = strrchr( path, '\\' ); /* find last / */
+#else
+ s = strrchr( path, '/' ); /* find last / */
+#endif
if ( s == NULL ) { /* no directory, so use current */
#if defined(__bsd__) /* get current working directory */
extern char *getwd( );
@@ -171,7 +184,18 @@ static int split( char *path, char *dirc, char *name )
dirc[l1-l2] = 0; /* add zero */
}
l1 = strlen( dirc ); /* length of directory */
+#ifdef windows
if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
+#else
+ if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
+#endif
+ s = strrchr( name, '.' ); /* find last / */
+ s++;
+ strcpy(ext,s); /* save extension */
+ l1= strlen( name);
+ l2= strlen( s)+1;
+ strncpy( finame, name, l1-l2);
+ finame[l1-l2]= 0;
return( 0 ); /* we're done */
}
@@ -719,7 +743,7 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
}
- ps[i][j]=(s2);
+ ps[i][j]=s2;
}
}
/*ps[3][2]=1;*/
@@ -1859,8 +1883,7 @@ fclose(ficresprob);
/**************** Main Program *****************/
/***********************************************/
-/*int main(int argc, char *argv[])*/
-int main()
+int main(int argc, char *argv[])
{
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
@@ -1876,7 +1899,10 @@ int main()
char line[MAXLINE], linepar[MAXLINE];
char title[MAXLINE];
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
- char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];
+ char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
+
+ char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];;
+
char filerest[FILENAMELENGTH];
char fileregp[FILENAMELENGTH];
char popfile[FILENAMELENGTH];
@@ -1909,7 +1935,7 @@ int main()
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate;
double yp,yp1,yp2;
- char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";
+ char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";
char *alph[]={"a","a","b","c","d","e"}, str[4];
@@ -1924,28 +1950,29 @@ int main()
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
- printf("\nIMACH, Version 0.7");
- printf("\nEnter the parameter file name: ");
-
-#ifdef windows
- scanf("%s",pathtot);
- getcwd(pathcd, size);
+ printf("\n%s",version);
+ if(argc <=1){
+ printf("\nEnter the parameter file name: ");
+ scanf("%s",pathtot);
+ }
+ else{
+ strcpy(pathtot,argv[1]);
+ }
+ /*if(getcwd(pathcd, 80)!= NULL)printf ("Error pathcd\n");*/
/*cygwin_split_path(pathtot,path,optionfile);
printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
/* cutv(path,optionfile,pathtot,'\\');*/
-split(pathtot, path,optionfile);
+ split(pathtot,path,optionfile,optionfilext,optionfilefiname);
+ printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
chdir(path);
replace(pathc,path);
-#endif
-#ifdef unix
- scanf("%s",optionfile);
-#endif
/*-------- arguments in the command line --------*/
strcpy(fileres,"r");
- strcat(fileres, optionfile);
+ strcat(fileres, optionfilefiname);
+ strcat(fileres,".txt"); /* Other files have txt extension */
/*---------arguments file --------*/
@@ -2324,6 +2351,15 @@ printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx,
}
}
}
+
+
+ /*for(i=1; i <=m ;i++){
+ for(k=1; k <=cptcovn; k++){
+ printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);
+ }
+ printf("\n");
+ }
+ scanf("%d",i);*/
/* Calculates basic frequencies. Computes observed prevalence at single age
and prints on file fileres'p'. */
@@ -2423,9 +2459,8 @@ printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx,
bage = agemin;
fage = agemax;
}
-
- fprintf(ficres,"# agemin agemax for life expectancy.\n");
-
+
+ fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
@@ -2470,13 +2505,16 @@ fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf fina
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2);
- /*------------ gnuplot -------------*/
-chdir(pathcd);
- if((ficgp=fopen("graph.plt","w"))==NULL) {
- printf("Problem with file graph.gp");goto end;
- }
+
+ /*------------ gnuplot -------------*/
+ /*chdir(pathcd);*/
+ strcpy(optionfilegnuplot,optionfilefiname);
+ strcat(optionfilegnuplot,".plt");
+ if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
+ printf("Problem with file %s",optionfilegnuplot);goto end;
+ }
#ifdef windows
- fprintf(ficgp,"cd \"%s\" \n",pathc);
+ fprintf(ficgp,"cd \"%s\" \n",pathc);
#endif
m=pow(2,cptcoveff);
@@ -3122,14 +3160,20 @@ strcpy(fileresvpl,"vpl");
end:
#ifdef windows
- chdir(pathcd);
+ /* chdir(pathcd);*/
#endif
-
- system("..\\gp37mgw\\wgnuplot graph.plt");
+ /*system("wgnuplot graph.plt");*/
+ /*system("../gp37mgw/wgnuplot graph.plt");*/
+ /*system("cd ../gp37mgw");*/
+ /* system("..\\gp37mgw\\wgnuplot graph.plt");*/
+ strcpy(plotcmd,GNUPLOTPROGRAM);
+ strcat(plotcmd," ");
+ strcat(plotcmd,optionfilegnuplot);
+ system(plotcmd);
#ifdef windows
while (z[0] != 'q') {
- chdir(pathcd);
+ chdir(path);
printf("\nType e to edit output files, c to start again, and q for exiting: ");
scanf("%s",z);
if (z[0] == 'c') system("./imach");
--
2.43.0