#define ODIRSEPARATOR '\\'
#endif
-char version[80]="Imach version 0.8k, July 2002, INED-EUROREVES ";
+char version[80]="Imach version 0.9, November 2002, INED-EUROREVES ";
int erreur; /* Error number */
int nvar;
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
int mle, weightopt;
int **mw; /* mw[mi][i] is number of the mi wave for this individual */
int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
+int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between
+ * wave mi and wave mi+1 is not an exact multiple of stepm. */
double jmean; /* Mean space between 2 waves */
double **oldm, **newm, **savm; /* Working pointers to matrices */
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
/**************** split *************************/
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 );
- s= strrchr( path, DIRSEPARATOR ); /* find last / */
- if ( s == NULL ) { /* no directory, so use current */
- /*if(strrchr(path, ODIRSEPARATOR )==NULL)
- printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
+ char *ss; /* pointer */
+ int l1, l2; /* length counters */
+
+ l1 = strlen(path ); /* length of path */
+ if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
+ ss= strrchr( path, DIRSEPARATOR ); /* find last / */
+ if ( ss == NULL ) { /* no directory, so use current */
+ /*if(strrchr(path, ODIRSEPARATOR )==NULL)
+ printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
#if defined(__bsd__) /* get current working directory */
- extern char *getwd( );
+ extern char *getwd( );
- if ( getwd( dirc ) == NULL ) {
+ if ( getwd( dirc ) == NULL ) {
#else
- extern char *getcwd( );
+ extern char *getcwd( );
- if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
+ if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
#endif
- return( GLOCK_ERROR_GETCWD );
- }
- strcpy( name, path ); /* we've got it */
- } else { /* strip direcotry from path */
- s++; /* after this, the filename */
- l2 = strlen( s ); /* length of filename */
- if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
- strcpy( name, s ); /* save file name */
- strncpy( dirc, path, l1 - l2 ); /* now the directory */
- dirc[l1-l2] = 0; /* add zero */
- }
- l1 = strlen( dirc ); /* length of directory */
+ return( GLOCK_ERROR_GETCWD );
+ }
+ strcpy( name, path ); /* we've got it */
+ } else { /* strip direcotry from path */
+ ss++; /* after this, the filename */
+ l2 = strlen( ss ); /* length of filename */
+ if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
+ strcpy( name, ss ); /* save file name */
+ strncpy( dirc, path, l1 - l2 ); /* now the directory */
+ dirc[l1-l2] = 0; /* add zero */
+ }
+ l1 = strlen( dirc ); /* length of directory */
#ifdef windows
- if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
+ if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
#else
- if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
+ 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 */
+ ss = strrchr( name, '.' ); /* find last / */
+ ss++;
+ strcpy(ext,ss); /* save extension */
+ l1= strlen( name);
+ l2= strlen(ss)+1;
+ strncpy( finame, name, l1-l2);
+ finame[l1-l2]= 0;
+ return( 0 ); /* we're done */
}
{
fprintf(stderr,"ERREUR ...\n");
fprintf(stderr,"%s\n",error_text);
- exit(1);
+ exit(EXIT_FAILURE);
}
/*********************** vector *******************/
double *vector(int nl, int nh)
double **out;
double sw; /* Sum of weights */
double lli; /* Individual log likelihood */
+ int s1, s2;
+ double bbh;
long ipmx;
/*extern weight */
/* We are differentiating ll according to initial status */
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
for(mi=1; mi<= wav[i]-1; mi++){
for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ savm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
for(d=0; d<dh[mi][i]; d++){
newm=savm;
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
} /* end mult */
- lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
- /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
+ /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+ /* But now since version 0.9 we anticipate for bias and large stepm.
+ * If stepm is larger than one month (smallest stepm) and if the exact delay
+ * (in months) between two waves is not a multiple of stepm, we rounded to
+ * the nearest (and in case of equal distance, to the lowest) interval but now
+ * we keep into memory the bias bh[mi][i] and also the previous matrix product
+ * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
+ * probability in order to take into account the bias as a fraction of the way
+ * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
+ * -stepm/2 to stepm/2 .
+ * For stepm=1 the results are the same as for previous versions of Imach.
+ * For stepm > 1 the results are less biased than in previous versions.
+ */
+ s1=s[mw[mi][i]][i];
+ s2=s[mw[mi+1][i]][i];
+ bbh=(double)bh[mi][i]/(double)stepm;
+ lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));
+ /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));*/
+ /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
+ /*if(lli ==000.0)*/
+ /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
ipmx +=1;
sw += weight[i];
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
}
/************ Frequencies ********************/
-void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
+void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
{ /* Some frequencies */
int i, m, jk, k1,i1, j1, bool, z1,z2,j;
/************* Waves Concatenation ***************/
-void concatwav(int wav[], int **dh, 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.
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
- dh[m][i] of 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.
*/
jk= j/stepm;
jl= j -jk*stepm;
ju= j -(jk+1)*stepm;
- if(jl <= -ju)
+ if(jl <= -ju){
dh[mi][i]=jk;
- else
+ bh[mi][i]=jl;
+ }
+ else{
dh[mi][i]=jk+1;
- if(dh[mi][i]==0)
+ bh[mi][i]=ju;
+ }
+ if(dh[mi][i]==0){
dh[mi][i]=1; /* At least one step */
+ bh[mi][i]=ju; /* At least one step */
+ printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);
+ }
+ if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);
}
}
}
void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij)
{
/* Variance of prevalence limit */
- /* 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 **dnewm,**doldm;
int i, j, nhstepm, hstepm;
}
-
cov[1]=1;
tj=cptcoveff;
if (cptcovn<1) {tj=1;ncodemax[1]=1;}
for(t=1; t<=tj;t++){
for(i1=1; i1<=ncodemax[t];i1++){
j1++;
-
if (cptcovn>0) {
fprintf(ficresprob, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
-
+ free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
+ free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
+ free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+ free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+
pmij(pmmij,cov,ncovmodel,x,nlstate);
k=0;
varpij[i][j][(int)age] = doldm[i][j];
/*printf("\n%d ",(int)age);
- for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
- printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
- fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
- }*/
+ for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+ printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+ fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+ }*/
fprintf(ficresprob,"\n%d ",(int)age);
fprintf(ficresprobcov,"\n%d ",(int)age);
/* Confidence intervalle of pij */
/*
- fprintf(ficgp,"\nset noparametric;unset label");
- fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
- fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
- fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
- fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
- fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
- fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
+ fprintf(ficgp,"\nset noparametric;unset label");
+ fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
+ fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+
fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
+
fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
+ fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
+ fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
*/
/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
} /*l1 */
}/* k1 */
} /* loop covariates */
- free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
- free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
- free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
- free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
- free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
- free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
}
+ free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
+ free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
free_vector(xp,1,npar);
fclose(ficresprob);
fclose(ficresprobcov);
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
int *popage;
- double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double calagedate, agelim, kk1, kk2;
double *popeffectif,*popcount;
double ***p3mat,***tabpop,***tabpopprev;
double ***mobaverage;
int c, h , cpt,l;
int ju,jl, mi;
int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
- int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
+ int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,*tab;
int mobilav=0,popforecast=0;
int hstepm, nhstepm;
double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;
double *epj, vepp;
double kk1, kk2;
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
-
+ /*int *movingaverage; */
char *alph[]={"a","a","b","c","d","e"}, str[4];
char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
/* long total_usecs;
- struct timeval start_time, end_time;
+ struct timeval start_time, end_time;
- gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
+ gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
getcwd(pathcd, size);
printf("\n%s",version);
/* cutv(path,optionfile,pathtot,'\\');*/
split(pathtot,path,optionfile,optionfilext,optionfilefiname);
- printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",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);
-/*-------- arguments in the command line --------*/
+ /*-------- arguments in the command line --------*/
/* Log file */
strcat(filelog, optionfilefiname);
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=%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=%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 nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
-while((c=getc(ficpar))=='#' && c!= EOF){
+ while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);
puts(line);
ungetc(c,ficpar);
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
- for(i=1; i <=nlstate; i++)
+ for(i=1; i <=nlstate; i++)
for(j=1; j <=nlstate+ndeath-1; j++){
fscanf(ficpar,"%1d%1d",&i1,&j1);
fprintf(ficparo,"%1d%1d",i1,j1);
fprintf(ficparo,"\n");
}
- npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
+ npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
p=param[1][1];
fprintf(ficlog,"\n");
- /*-------- Rewriting paramater 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", fileres);goto end;
- fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
- }
- fprintf(ficres,"#%s\n",version);
+ /*-------- Rewriting paramater 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", fileres);goto end;
+ fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
+ }
+ fprintf(ficres,"#%s\n",version);
- /*-------- data file ----------*/
- if((fic=fopen(datafile,"r"))==NULL) {
- printf("Problem with datafile: %s\n", datafile);goto end;
- fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end;
- }
+ /*-------- data file ----------*/
+ if((fic=fopen(datafile,"r"))==NULL) {
+ printf("Problem with datafile: %s\n", datafile);goto end;
+ fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end;
+ }
- n= lastobs;
- severity = vector(1,maxwav);
- outcome=imatrix(1,maxwav+1,1,n);
- num=ivector(1,n);
- moisnais=vector(1,n);
- annais=vector(1,n);
- moisdc=vector(1,n);
- andc=vector(1,n);
- agedc=vector(1,n);
- cod=ivector(1,n);
- weight=vector(1,n);
- for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
- mint=matrix(1,maxwav,1,n);
- anint=matrix(1,maxwav,1,n);
- s=imatrix(1,maxwav+1,1,n);
- adl=imatrix(1,maxwav+1,1,n);
- tab=ivector(1,NCOVMAX);
- ncodemax=ivector(1,8);
-
- i=1;
- while (fgets(line, MAXLINE, fic) != NULL) {
- if ((i >= firstobs) && (i <=lastobs)) {
+ n= lastobs;
+ severity = vector(1,maxwav);
+ outcome=imatrix(1,maxwav+1,1,n);
+ num=ivector(1,n);
+ moisnais=vector(1,n);
+ annais=vector(1,n);
+ moisdc=vector(1,n);
+ andc=vector(1,n);
+ agedc=vector(1,n);
+ cod=ivector(1,n);
+ weight=vector(1,n);
+ for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
+ mint=matrix(1,maxwav,1,n);
+ anint=matrix(1,maxwav,1,n);
+ s=imatrix(1,maxwav+1,1,n);
+ tab=ivector(1,NCOVMAX);
+ ncodemax=ivector(1,8);
+
+ i=1;
+ while (fgets(line, MAXLINE, fic) != NULL) {
+ if ((i >= firstobs) && (i <=lastobs)) {
- for (j=maxwav;j>=1;j--){
- cutv(stra, strb,line,' '); s[j][i]=atoi(strb);
- strcpy(line,stra);
- cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
- cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
- }
+ for (j=maxwav;j>=1;j--){
+ cutv(stra, strb,line,' '); s[j][i]=atoi(strb);
+ strcpy(line,stra);
+ cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
+ cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
+ }
- cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra);
- cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra);
+ cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra);
+ cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra);
- cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra);
- cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);
+ cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra);
+ cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);
- cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);
- for (j=ncovcol;j>=1;j--){
- cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
- }
- num[i]=atol(stra);
+ cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);
+ for (j=ncovcol;j>=1;j--){
+ cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
+ }
+ num[i]=atol(stra);
- /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
- printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
+ /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
+ printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
- i=i+1;
- }
- }
- /* printf("ii=%d", ij);
- scanf("%d",i);*/
+ i=i+1;
+ }
+ }
+ /* printf("ii=%d", ij);
+ scanf("%d",i);*/
imx=i-1; /* Number of individuals */
/* for (i=1; i<=imx; i++){
goto end;
}
- /* This loop fill the array Tvar from the string 'model'.*/
+ /* This loop fills the array Tvar from the string 'model'.*/
for(i=(j+1); i>=1;i--){
cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */
- if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */
+ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
/*scanf("%d",i);*/
if (strchr(strb,'*')) { /* Model includes a product */
fclose(fic);*/
/* if(mle==1){*/
- if (weightopt != 1) { /* Maximisation without weights*/
- for(i=1;i<=n;i++) weight[i]=1.0;
- }
+ if (weightopt != 1) { /* Maximisation without weights*/
+ for(i=1;i<=n;i++) weight[i]=1.0;
+ }
/*-calculation of age at interview from date of interview and age at death -*/
- agev=matrix(1,maxwav,1,imx);
+ agev=matrix(1,maxwav,1,imx);
- for (i=1; i<=imx; i++) {
- for(m=2; (m<= maxwav); m++) {
- if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
- anint[m][i]=9999;
- s[m][i]=-1;
- }
- if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;
+ for (i=1; i<=imx; i++) {
+ for(m=2; (m<= maxwav); m++) {
+ if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
+ anint[m][i]=9999;
+ s[m][i]=-1;
}
+ if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;
}
+ }
- for (i=1; i<=imx; i++) {
- agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
- for(m=1; (m<= maxwav); m++){
- if(s[m][i] >0){
- if (s[m][i] >= nlstate+1) {
- if(agedc[i]>0)
- if(moisdc[i]!=99 && andc[i]!=9999)
- agev[m][i]=agedc[i];
- /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
- else {
+ for (i=1; i<=imx; i++) {
+ agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
+ for(m=1; (m<= maxwav); m++){
+ if(s[m][i] >0){
+ if (s[m][i] >= nlstate+1) {
+ if(agedc[i]>0)
+ if(moisdc[i]!=99 && andc[i]!=9999)
+ agev[m][i]=agedc[i];
+ /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
+ else {
if (andc[i]!=9999){
- printf("Warning negative age at death: %d line:%d\n",num[i],i);
- fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
- agev[m][i]=-1;
+ printf("Warning negative age at death: %d line:%d\n",num[i],i);
+ fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
+ agev[m][i]=-1;
}
}
- }
- else if(s[m][i] !=9){ /* Should no more exist */
- agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
- if(mint[m][i]==99 || anint[m][i]==9999)
- agev[m][i]=1;
- else if(agev[m][i] <agemin){
- agemin=agev[m][i];
- /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/
- }
- else if(agev[m][i] >agemax){
- agemax=agev[m][i];
- /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/
- }
- /*agev[m][i]=anint[m][i]-annais[i];*/
- /* agev[m][i] = age[i]+2*m;*/
- }
- else { /* =9 */
+ }
+ else if(s[m][i] !=9){ /* Should no more exist */
+ agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
+ if(mint[m][i]==99 || anint[m][i]==9999)
agev[m][i]=1;
- s[m][i]=-1;
+ else if(agev[m][i] <agemin){
+ agemin=agev[m][i];
+ /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/
}
+ else if(agev[m][i] >agemax){
+ agemax=agev[m][i];
+ /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/
+ }
+ /*agev[m][i]=anint[m][i]-annais[i];*/
+ /* agev[m][i] = age[i]+2*m;*/
}
- else /*= 0 Unknown */
+ else { /* =9 */
agev[m][i]=1;
+ s[m][i]=-1;
+ }
}
-
+ else /*= 0 Unknown */
+ agev[m][i]=1;
}
- for (i=1; i<=imx; i++) {
- for(m=1; (m<= maxwav); m++){
- if (s[m][i] > (nlstate+ndeath)) {
- printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
- fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
- goto end;
- }
+
+ }
+ for (i=1; i<=imx; i++) {
+ for(m=1; (m<= maxwav); m++){
+ if (s[m][i] > (nlstate+ndeath)) {
+ printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
+ fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
+ goto end;
}
}
+ }
- printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
- fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
+ printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
+ fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
- free_vector(severity,1,maxwav);
- free_imatrix(outcome,1,maxwav+1,1,n);
- free_vector(moisnais,1,n);
- free_vector(annais,1,n);
- /* free_matrix(mint,1,maxwav,1,n);
- free_matrix(anint,1,maxwav,1,n);*/
- free_vector(moisdc,1,n);
- free_vector(andc,1,n);
+ free_vector(severity,1,maxwav);
+ free_imatrix(outcome,1,maxwav+1,1,n);
+ free_vector(moisnais,1,n);
+ free_vector(annais,1,n);
+ /* free_matrix(mint,1,maxwav,1,n);
+ free_matrix(anint,1,maxwav,1,n);*/
+ free_vector(moisdc,1,n);
+ free_vector(andc,1,n);
- wav=ivector(1,imx);
- dh=imatrix(1,lastpass-firstpass+1,1,imx);
- mw=imatrix(1,lastpass-firstpass+1,1,imx);
+ wav=ivector(1,imx);
+ dh=imatrix(1,lastpass-firstpass+1,1,imx);
+ bh=imatrix(1,lastpass-firstpass+1,1,imx);
+ mw=imatrix(1,lastpass-firstpass+1,1,imx);
- /* Concatenates waves */
- concatwav(wav, dh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm);
+ /* Concatenates waves */
+ concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm);
- /* 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 */
- Tcode=ivector(1,100);
- nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
- ncodemax[1]=1;
- if (cptcovn > 0) tricode(Tvar,nbcode,imx);
+ Tcode=ivector(1,100);
+ nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
+ ncodemax[1]=1;
+ if (cptcovn > 0) tricode(Tvar,nbcode,imx);
- codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of
- the estimations*/
- h=0;
- m=pow(2,cptcoveff);
+ codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of
+ the estimations*/
+ h=0;
+ m=pow(2,cptcoveff);
- for(k=1;k<=cptcoveff; k++){
- for(i=1; i <=(m/pow(2,k));i++){
- for(j=1; j <= ncodemax[k]; j++){
- for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
- h++;
- if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
- /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/
- }
- }
- }
- }
- /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]);
- codtab[1][2]=1;codtab[2][2]=2; */
- /* 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");
+ for(k=1;k<=cptcoveff; k++){
+ for(i=1; i <=(m/pow(2,k));i++){
+ for(j=1; j <= ncodemax[k]; j++){
+ for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
+ h++;
+ if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
+ /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/
+ }
}
- scanf("%d",i);*/
+ }
+ }
+ /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]);
+ codtab[1][2]=1;codtab[2][2]=2; */
+ /* 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'. */
+ /* Calculates basic frequencies. Computes observed prevalence at single age
+ and prints on file fileres'p'. */
- pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
-
- /* For Powell, parameters are in a vector p[] starting at p[1]
- so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
- p=param[1][1]; /* *(*(*(param +1)+1)+0) */
+ /* For Powell, parameters are in a vector p[] starting at p[1]
+ so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
+ p=param[1][1]; /* *(*(*(param +1)+1)+0) */
- if(mle==1){
+ if(mle==1){
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
- }
+ }
- /*--------- results files --------------*/
- fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
+ /*--------- results files --------------*/
+ fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
- jk=1;
- fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
- printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
- fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
- for(i=1,jk=1; i <=nlstate; i++){
- for(k=1; k <=(nlstate+ndeath); k++){
- if (k != i)
- {
- printf("%d%d ",i,k);
- fprintf(ficlog,"%d%d ",i,k);
- fprintf(ficres,"%1d%1d ",i,k);
- for(j=1; j <=ncovmodel; j++){
- printf("%f ",p[jk]);
- fprintf(ficlog,"%f ",p[jk]);
- fprintf(ficres,"%f ",p[jk]);
- jk++;
- }
- printf("\n");
- fprintf(ficlog,"\n");
- fprintf(ficres,"\n");
- }
- }
- }
- if(mle==1){
- /* Computing hessian and covariance matrix */
- ftolhess=ftol; /* Usually correct */
- hesscov(matcov, p, npar, delti, ftolhess, func);
- }
- fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
- printf("# Scales (for hessian or gradient estimation)\n");
- fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
- for(i=1,jk=1; i <=nlstate; i++){
- for(j=1; j <=nlstate+ndeath; j++){
- if (j!=i) {
- fprintf(ficres,"%1d%1d",i,j);
- printf("%1d%1d",i,j);
- fprintf(ficlog,"%1d%1d",i,j);
- for(k=1; k<=ncovmodel;k++){
- printf(" %.5e",delti[jk]);
- fprintf(ficlog," %.5e",delti[jk]);
- fprintf(ficres," %.5e",delti[jk]);
- jk++;
- }
- printf("\n");
- fprintf(ficlog,"\n");
- fprintf(ficres,"\n");
- }
- }
- }
+ jk=1;
+ fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+ printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+ fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+ for(i=1,jk=1; i <=nlstate; i++){
+ for(k=1; k <=(nlstate+ndeath); k++){
+ if (k != i)
+ {
+ printf("%d%d ",i,k);
+ fprintf(ficlog,"%d%d ",i,k);
+ fprintf(ficres,"%1d%1d ",i,k);
+ for(j=1; j <=ncovmodel; j++){
+ printf("%f ",p[jk]);
+ fprintf(ficlog,"%f ",p[jk]);
+ fprintf(ficres,"%f ",p[jk]);
+ jk++;
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
+ fprintf(ficres,"\n");
+ }
+ }
+ }
+ if(mle==1){
+ /* Computing hessian and covariance matrix */
+ ftolhess=ftol; /* Usually correct */
+ hesscov(matcov, p, npar, delti, ftolhess, func);
+ }
+ fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
+ printf("# Scales (for hessian or gradient estimation)\n");
+ fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
+ for(i=1,jk=1; i <=nlstate; i++){
+ for(j=1; j <=nlstate+ndeath; j++){
+ if (j!=i) {
+ fprintf(ficres,"%1d%1d",i,j);
+ printf("%1d%1d",i,j);
+ fprintf(ficlog,"%1d%1d",i,j);
+ for(k=1; k<=ncovmodel;k++){
+ printf(" %.5e",delti[jk]);
+ fprintf(ficlog," %.5e",delti[jk]);
+ fprintf(ficres," %.5e",delti[jk]);
+ jk++;
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
+ fprintf(ficres,"\n");
+ }
+ }
+ }
- k=1;
- fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
- if(mle==1)
- printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
- fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
- for(i=1;i<=npar;i++){
- /* if (k>nlstate) k=1;
- i1=(i-1)/(ncovmodel*nlstate)+1;
- fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
- printf("%s%d%d",alph[k],i1,tab[i]);*/
- fprintf(ficres,"%3d",i);
- if(mle==1)
- printf("%3d",i);
- fprintf(ficlog,"%3d",i);
- for(j=1; j<=i;j++){
- fprintf(ficres," %.5e",matcov[i][j]);
- if(mle==1)
- printf(" %.5e",matcov[i][j]);
- fprintf(ficlog," %.5e",matcov[i][j]);
- }
- fprintf(ficres,"\n");
- if(mle==1)
- printf("\n");
- fprintf(ficlog,"\n");
- k++;
- }
+ fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
+ if(mle==1)
+ printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
+ fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
+ for(i=1,k=1;i<=npar;i++){
+ /* if (k>nlstate) k=1;
+ i1=(i-1)/(ncovmodel*nlstate)+1;
+ fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
+ printf("%s%d%d",alph[k],i1,tab[i]);
+ */
+ fprintf(ficres,"%3d",i);
+ if(mle==1)
+ printf("%3d",i);
+ fprintf(ficlog,"%3d",i);
+ for(j=1; j<=i;j++){
+ fprintf(ficres," %.5e",matcov[i][j]);
+ if(mle==1)
+ printf(" %.5e",matcov[i][j]);
+ fprintf(ficlog," %.5e",matcov[i][j]);
+ }
+ fprintf(ficres,"\n");
+ if(mle==1)
+ printf("\n");
+ fprintf(ficlog,"\n");
+ k++;
+ }
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
- puts(line);
- fputs(line,ficparo);
- }
- ungetc(c,ficpar);
- estepm=0;
- fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
- if (estepm==0 || estepm < stepm) estepm=stepm;
- if (fage <= 2) {
- bage = ageminpar;
- fage = agemaxpar;
- }
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ puts(line);
+ fputs(line,ficparo);
+ }
+ ungetc(c,ficpar);
+
+ estepm=0;
+ fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
+ if (estepm==0 || estepm < stepm) estepm=stepm;
+ if (fage <= 2) {
+ bage = ageminpar;
+ fage = agemaxpar;
+ }
- 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 estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
- fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+ 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 estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+ fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
- puts(line);
- fputs(line,ficparo);
- }
- ungetc(c,ficpar);
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ puts(line);
+ fputs(line,ficparo);
+ }
+ ungetc(c,ficpar);
- fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
- fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
- fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+ fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
+ fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+ fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
- puts(line);
- fputs(line,ficparo);
- }
- ungetc(c,ficpar);
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ puts(line);
+ fputs(line,ficparo);
+ }
+ ungetc(c,ficpar);
- dateprev1=anprev1+mprev1/12.+jprev1/365.;
- dateprev2=anprev2+mprev2/12.+jprev2/365.;
+ dateprev1=anprev1+mprev1/12.+jprev1/365.;
+ dateprev2=anprev2+mprev2/12.+jprev2/365.;
fscanf(ficpar,"pop_based=%d\n",&popbased);
fprintf(ficparo,"pop_based=%d\n",popbased);
ungetc(c,ficpar);
fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);
-fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
-fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
+ fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
+ fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
-while((c=getc(ficpar))=='#' && c!= EOF){
+ while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);
puts(line);
fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
- freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
+ freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
-/*------------ gnuplot -------------*/
- strcpy(optionfilegnuplot,optionfilefiname);
- strcat(optionfilegnuplot,".gp");
- if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
- printf("Problem with file %s",optionfilegnuplot);
- }
- else{
- fprintf(ficgp,"\n# %s\n", version);
- fprintf(ficgp,"# %s\n", optionfilegnuplot);
- fprintf(ficgp,"set missing 'NaNq'\n");
-}
- fclose(ficgp);
- printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
-/*--------- index.htm --------*/
+ /*------------ gnuplot -------------*/
+ strcpy(optionfilegnuplot,optionfilefiname);
+ strcat(optionfilegnuplot,".gp");
+ if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
+ printf("Problem with file %s",optionfilegnuplot);
+ }
+ else{
+ fprintf(ficgp,"\n# %s\n", version);
+ fprintf(ficgp,"# %s\n", optionfilegnuplot);
+ fprintf(ficgp,"set missing 'NaNq'\n");
+ }
+ fclose(ficgp);
+ printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
+ /*--------- index.htm --------*/
strcpy(optionfilehtm,optionfile);
strcat(optionfilehtm,".htm");
- Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);
fclose(fichtm);
- printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
+ printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
-/*------------ free_vector -------------*/
- chdir(path);
+ /*------------ free_vector -------------*/
+ chdir(path);
- free_ivector(wav,1,imx);
- free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
- free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
- free_ivector(num,1,n);
- free_vector(agedc,1,n);
- /*free_matrix(covar,1,NCOVMAX,1,n);*/
- fclose(ficparo);
- fclose(ficres);
+ free_ivector(wav,1,imx);
+ free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
+ free_imatrix(bh,1,lastpass-firstpass+1,1,imx);
+ free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
+ free_ivector(num,1,n);
+ free_vector(agedc,1,n);
+ free_matrix(covar,0,NCOVMAX,1,n);
+ /*free_matrix(covar,1,NCOVMAX,1,n);*/
+ fclose(ficparo);
+ fclose(ficres);
/*--------------- Prevalence limit (stable prevalence) --------------*/
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
- k=0;
+
agebase=ageminpar;
agelim=agemaxpar;
ftolpl=1.e-10;
i1=cptcoveff;
if (cptcovn < 1){i1=1;}
- for(cptcov=1;cptcov<=i1;cptcov++){
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
- fprintf(ficrespl,"\n#******");
- printf("\n#******");
- fprintf(ficlog,"\n#******");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- }
- fprintf(ficrespl,"******\n");
- printf("******\n");
- fprintf(ficlog,"******\n");
+ k=k+1;
+ /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
+ fprintf(ficrespl,"\n#******");
+ printf("\n#******");
+ fprintf(ficlog,"\n#******");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ }
+ fprintf(ficrespl,"******\n");
+ printf("******\n");
+ fprintf(ficlog,"******\n");
- for (age=agebase; age<=agelim; age++){
- prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
- fprintf(ficrespl,"%.0f",age );
- for(i=1; i<=nlstate;i++)
+ for (age=agebase; age<=agelim; age++){
+ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+ fprintf(ficrespl,"%.0f",age );
+ for(i=1; i<=nlstate;i++)
fprintf(ficrespl," %.5f", prlim[i][i]);
- fprintf(ficrespl,"\n");
- }
+ fprintf(ficrespl,"\n");
}
}
+ }
fclose(ficrespl);
/*------------- h Pij x at various ages ------------*/
/* hstepm=1; aff par mois*/
- k=0;
- for(cptcov=1;cptcov<=i1;cptcov++){
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
- fprintf(ficrespij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- fprintf(ficrespij,"******\n");
+ fprintf(ficrespij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ fprintf(ficrespij,"******\n");
- for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
- nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+ for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
+ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
- /* nhstepm=nhstepm*YEARM; aff par mois*/
+ /* nhstepm=nhstepm*YEARM; aff par mois*/
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- oldm=oldms;savm=savms;
- hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
- fprintf(ficrespij,"# Age");
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+ fprintf(ficrespij,"# Age");
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate+ndeath;j++)
+ fprintf(ficrespij," %1d-%1d",i,j);
+ fprintf(ficrespij,"\n");
+ for (h=0; h<=nhstepm; h++){
+ fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespij," %1d-%1d",i,j);
- fprintf(ficrespij,"\n");
- for (h=0; h<=nhstepm; h++){
- fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespij," %.5f", p3mat[i][j][h]);
- fprintf(ficrespij,"\n");
- }
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespij," %.5f", p3mat[i][j][h]);
fprintf(ficrespij,"\n");
}
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespij,"\n");
+ }
}
}
}
}
- k=0;
- for(cptcov=1;cptcov<=i1;cptcov++){
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
fprintf(ficrest,"\n#****** ");
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav);
if(popbased==1){
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav);
- }
+ }
fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
}
fprintf(ficrest,"\n");
}
+ free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ free_vector(epj,1,nlstate+1);
}
}
-free_matrix(mint,1,maxwav,1,n);
- free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
- free_vector(weight,1,n);
+ free_vector(weight,1,n);
+ free_imatrix(Tvard,1,15,1,2);
+ free_imatrix(s,1,maxwav+1,1,n);
+ free_matrix(anint,1,maxwav,1,n);
+ free_matrix(mint,1,maxwav,1,n);
+ free_ivector(cod,1,n);
+ free_ivector(tab,1,NCOVMAX);
fclose(ficreseij);
fclose(ficresvij);
fclose(ficrest);
fclose(ficpar);
- free_vector(epj,1,nlstate+1);
/*------- Variance of stable prevalence------*/
}
printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
- k=0;
- for(cptcov=1;cptcov<=i1;cptcov++){
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
fprintf(ficresvpl,"\n#****** ");
varpl=matrix(1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
- varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
+ varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
+ free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
}
- }
+ }
fclose(ficresvpl);
/*---------- End : free ----------------*/
- free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
-
- free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
-
-
free_matrix(pmmij,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(agev,1,maxwav,1,imx);
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ free_ivector(ncodemax,1,8);
+ free_ivector(Tvar,1,15);
+ free_ivector(Tprod,1,15);
+ free_ivector(Tvaraff,1,15);
+ free_ivector(Tage,1,15);
+ free_ivector(Tcode,1,100);
fprintf(fichtm,"\n</body>");
fclose(fichtm);
/*printf("Total time was %d uSec.\n", total_usecs);*/
/*------ End -----------*/
-
- end:
+ end:
#ifdef windows
/* chdir(pathcd);*/
#endif
/*system("../gp37mgw/wgnuplot graph.plt");*/
/*system("cd ../gp37mgw");*/
/* system("..\\gp37mgw\\wgnuplot graph.plt");*/
- strcpy(plotcmd,GNUPLOTPROGRAM);
- strcat(plotcmd," ");
- strcat(plotcmd,optionfilegnuplot);
- printf("Starting: %s\n",plotcmd);fflush(stdout);
- system(plotcmd);
+ strcpy(plotcmd,GNUPLOTPROGRAM);
+ strcat(plotcmd," ");
+ strcat(plotcmd,optionfilegnuplot);
+ printf("Starting: %s\n",plotcmd);fflush(stdout);
+ system(plotcmd);
/*#ifdef windows*/
while (z[0] != 'q') {
git://henry.ined.fr / .git / commitdiff