--- imach/src/imach.c 2006/01/24 19:37:15 1.109
+++ imach/src/imach.c 2006/03/06 10:29:27 1.116
@@ -1,6 +1,33 @@
-/* $Id: imach.c,v 1.109 2006/01/24 19:37:15 brouard Exp $
+/* $Id: imach.c,v 1.116 2006/03/06 10:29:27 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.116 2006/03/06 10:29:27 brouard
+ (Module): Variance-covariance wrong links and
+ varian-covariance of ej. is needed (Saito).
+
+ Revision 1.115 2006/02/27 12:17:45 brouard
+ (Module): One freematrix added in mlikeli! 0.98c
+
+ Revision 1.114 2006/02/26 12:57:58 brouard
+ (Module): Some improvements in processing parameter
+ filename with strsep.
+
+ Revision 1.113 2006/02/24 14:20:24 brouard
+ (Module): Memory leaks checks with valgrind and:
+ datafile was not closed, some imatrix were not freed and on matrix
+ allocation too.
+
+ Revision 1.112 2006/01/30 09:55:26 brouard
+ (Module): Back to gnuplot.exe instead of wgnuplot.exe
+
+ Revision 1.111 2006/01/25 20:38:18 brouard
+ (Module): Lots of cleaning and bugs added (Gompertz)
+ (Module): Comments can be added in data file. Missing date values
+ can be a simple dot '.'.
+
+ Revision 1.110 2006/01/25 00:51:50 brouard
+ (Module): Lots of cleaning and bugs added (Gompertz)
+
Revision 1.109 2006/01/24 19:37:15 brouard
(Module): Comments (lines starting with a #) are allowed in data.
@@ -281,11 +308,11 @@ extern int errno;
#define ODIRSEPARATOR '/'
#endif
-/* $Id: imach.c,v 1.109 2006/01/24 19:37:15 brouard Exp $ */
+/* $Id: imach.c,v 1.116 2006/03/06 10:29:27 brouard Exp $ */
/* $State: Exp $ */
-char version[]="Imach version 0.98a, January 2006, INED-EUROREVES ";
-char fullversion[]="$Revision: 1.109 $ $Date: 2006/01/24 19:37:15 $";
+char version[]="Imach version 0.98c, February 2006, INED-EUROREVES ";
+char fullversion[]="$Revision: 1.116 $ $Date: 2006/03/06 10:29:27 $";
int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
int nvar;
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -298,6 +325,7 @@ int popbased=0;
int *wav; /* Number of waves for this individuual 0 is possible */
int maxwav; /* Maxim number of waves */
int jmin, jmax; /* min, max spacing between 2 waves */
+int ijmin, ijmax; /* Individuals having jmin and jmax */
int gipmx, gsw; /* Global variables on the number of contributions
to the likelihood and the sum of weights (done by funcone)*/
int mle, weightopt;
@@ -328,7 +356,7 @@ FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];
char title[MAXLINE];
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
-char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
+char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH];
char command[FILENAMELENGTH];
int outcmd=0;
@@ -1647,6 +1675,7 @@ void mlikeli(FILE *ficres,double p[], in
powell(p,xi,npar,ftol,&iter,&fret,func);
+ free_matrix(xi,1,npar,1,npar);
fclose(ficrespow);
printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
@@ -2201,11 +2230,11 @@ void concatwav(int wav[], int **dh, int
if(mi==0){
nbwarn++;
if(first==0){
- printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i);
+ printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
first=1;
}
if(first==1){
- fprintf(ficlog,"Warning! None valid information for:%ld line=%d (skipped)\n",num[i],i);
+ fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
}
} /* end mi==0 */
} /* End individuals */
@@ -2228,8 +2257,14 @@ void concatwav(int wav[], int **dh, int
fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
}
k=k+1;
- if (j >= jmax) jmax=j;
- if (j <= jmin) jmin=j;
+ if (j >= jmax){
+ jmax=j;
+ ijmax=i;
+ }
+ if (j <= jmin){
+ jmin=j;
+ ijmin=i;
+ }
sum=sum+j;
/*if (j<0) printf("j=%d num=%d \n",j,i);*/
/* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
@@ -2240,8 +2275,14 @@ void concatwav(int wav[], int **dh, int
/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
k=k+1;
- if (j >= jmax) jmax=j;
- else if (j <= jmin)jmin=j;
+ if (j >= jmax) {
+ jmax=j;
+ ijmax=i;
+ }
+ else if (j <= jmin){
+ jmin=j;
+ ijmin=i;
+ }
/* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
/*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
if(j<0){
@@ -2284,8 +2325,8 @@ void concatwav(int wav[], int **dh, int
} /* end wave */
}
jmean=sum/k;
- printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
- fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
+ printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);
+ fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
}
/*********** Tricode ****************************/
@@ -3210,6 +3251,8 @@ To be simple, these graphs help to under
}
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_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
+ free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar);
free_vector(xp,1,npar);
fclose(ficresprob);
fclose(ficresprobcov);
@@ -3241,9 +3284,18 @@ void printinghtml(char fileres[], char t
- Stable prevalence in each health state: %s
\n",
subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
fprintf(fichtm,"\
- - Life expectancies by age and initial health status (estepm=%2d months): \
+ - Life expectancies by age and initial health status (estepm=%2d months) WRONG LINK (to be made): \
+ %s
\n",
+ estepm,subdirf2(fileres,"le"),subdirf2(fileres,"le"));
+ fprintf(fichtm,"\
+ - Health expectancies by age and initial health status with standard errors (estepm=%2d months): \
%s
\n",
estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
+ fprintf(fichtm,"\
+ - Variances and covariances of health expectancies by age and initial health status (estepm=%2d months) TO BE MADE: \
+ %s
\n",
+ estepm,subdirf2(fileres,"vch"),subdirf2(fileres,"vch"));
+
fprintf(fichtm," \n- Graphs
");
@@ -3295,10 +3347,10 @@ fprintf(fichtm," \n
- Graphs
- Correlation matrix of one-step probabilities: %s
\n",
subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
fprintf(fichtm,"\
- - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): %s
\n",
+ - Variances and covariances of health expectancies by age (estepm=%d months): %s
\n",
estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
fprintf(fichtm,"\
- - Health expectancies with their variances (no covariance): %s
\n",
+ - Life and health expectancies with their standard errors: %s
\n",
subdirf2(fileres,"t"),subdirf2(fileres,"t"));
fprintf(fichtm,"\
- Standard deviation of stable prevalences: %s
\n",\
@@ -3996,6 +4048,7 @@ double gompertz(double x[])
{
double A,B,L=0.0,sump=0.,num=0.;
int i,n=0; /* n is the size of the sample */
+
for (i=0;i<=imx-1 ; i++) {
sump=sump+weight[i];
/* sump=sump+1;*/
@@ -4008,14 +4061,15 @@ double gompertz(double x[])
for (i=1;i<=imx ; i++)
{
- if (cens[i]==1 & wav[i]>1)
+ if (cens[i] == 1 && wav[i]>1)
A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
- if (cens[i]==0 & wav[i]>1)
+ if (cens[i] == 0 && wav[i]>1)
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);
- if (wav[i]>1 & agecens[i]>15) {
+ /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
+ if (wav[i] > 1 ) { /* ??? */
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]);*/
}
@@ -4084,13 +4138,14 @@ int main(int argc, char *argv[])
{
int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
- int linei;
+ int linei, month, year,iout;
int jj, ll, li, lj, lk, imk;
int numlinepar=0; /* Current linenumber of parameter file */
int itimes;
int NDIM=2;
char ca[32], cb[32], cc[32];
+ char dummy[]=" ";
/* FILE *fichtm; *//* Html File */
/* FILE *ficgp;*/ /*Gnuplot File */
struct stat info;
@@ -4107,6 +4162,7 @@ int main(int argc, char *argv[])
char line[MAXLINE], linepar[MAXLINE];
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
char pathr[MAXLINE], pathimach[MAXLINE];
+ char **bp, *tok, *val; /* pathtot */
int firstobs=1, lastobs=10;
int sdeb, sfin; /* Status at beginning and end */
int c, h , cpt,l;
@@ -4184,7 +4240,17 @@ int main(int argc, char *argv[])
printf("\n%s\n%s",version,fullversion);
if(argc <=1){
printf("\nEnter the parameter file name: ");
- scanf("%s",pathtot);
+ fgets(pathr,FILENAMELENGTH,stdin);
+ i=strlen(pathr);
+ if(pathr[i-1]=='\n')
+ pathr[i-1]='\0';
+ for (tok = pathr; tok != NULL; ){
+ printf("Pathr |%s|\n",pathr);
+ while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
+ printf("val= |%s| pathr=%s\n",val,pathr);
+ strcpy (pathtot, val);
+ if(pathr[0] == '\0') break; /* Un peu sale */
+ }
}
else{
strcpy(pathtot,argv[1]);
@@ -4492,131 +4558,116 @@ int main(int argc, char *argv[])
i=1;
linei=0;
- while ((fgets(line, MAXLINE, fic) != NULL) ||((i >= firstobs) && (i <=lastobs))) {
+ while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
linei=linei+1;
- printf("IIIII= %d linei=%d\n",i,linei);
for(j=strlen(line); j>=0;j--){ /* Untabifies line */
- if(line[j] == '\t')
- line[j] = ' ';
- }
- for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10);j--){;};line[j+1]=0; /* Trims blanks at end of line */
- if(line[0]=='#'){
- fprintf(ficlog,"Comment line\n%s\n",line);
- printf("Comment line\n%s\n",line);
- continue;
- }
- for (j=maxwav;j>=1;j--){
- cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
- /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %d %s for individual %d\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n",lval, i,line,linei,j,maxwav);
- exit(1);
- }
- s[j][i]=lval;
-
- strcpy(line,stra);
- cutv(stra, strb,line,'/');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d. Exiting.\n",lval, i,line,linei,j);
- exit(1);
- }
- anint[j][i]=(double)(lval);
-
- strcpy(line,stra);
- cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a month of exam at wave %d. Exiting.\n",lval, i,line, linei,j);
- exit(1);
- }
- mint[j][i]=(double)(lval);
- strcpy(line,stra);
- }
-
- cutv(stra, strb,line,'/');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a year of death. Exiting.\n",lval, i,line,linei);
- exit(1);
- }
- andc[i]=(double)(lval);
- strcpy(line,stra);
+ if(line[j] == '\t')
+ line[j] = ' ';
+ }
+ for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){
+ ;
+ };
+ line[j+1]=0; /* Trims blanks at end of line */
+ if(line[0]=='#'){
+ fprintf(ficlog,"Comment line\n%s\n",line);
+ printf("Comment line\n%s\n",line);
+ continue;
+ }
+ for (j=maxwav;j>=1;j--){
cutv(stra, strb,line,' ');
errno=0;
lval=strtol(strb,&endptr,10);
+ /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a month of death. Exiting.\n",lval,i,line, linei);
+ printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
exit(1);
}
- moisdc[i]=(double)(lval);
-
- strcpy(line,stra);
- cutv(stra, strb,line,'/');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a year of birth. Exiting.\n",lval, i,line, linei);
- exit(1);
- }
- annais[i]=(double)(lval);
-
+ s[j][i]=lval;
+
strcpy(line,stra);
cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a month of birth. Exiting.\n",lval,i,line,linei);
+ if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+ }
+ else if(iout=sscanf(strb,"%s.") != 0){
+ month=99;
+ year=9999;
+ }else{
+ printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
exit(1);
}
- moisnais[i]=(double)(lval);
+ anint[j][i]= (double) year;
+ mint[j][i]= (double)month;
strcpy(line,stra);
-
+ } /* ENd Waves */
+
+ cutv(stra, strb,line,' ');
+ if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+ }
+ else if(iout=sscanf(strb,"%s.",dummy) != 0){
+ month=99;
+ year=9999;
+ }else{
+ printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
+ exit(1);
+ }
+ andc[i]=(double) year;
+ moisdc[i]=(double) month;
+ strcpy(line,stra);
+
+ cutv(stra, strb,line,' ');
+ if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+ }
+ else if(iout=sscanf(strb,"%s.") != 0){
+ month=99;
+ year=9999;
+ }else{
+ printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);
+ exit(1);
+ }
+ annais[i]=(double)(year);
+ moisnais[i]=(double)(month);
+ strcpy(line,stra);
+
+ cutv(stra, strb,line,' ');
+ errno=0;
+ lval=strtol(strb,&endptr,10);
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight. Exiting.\n",lval, i,line,linei);
+ exit(1);
+ }
+ weight[i]=(double)(lval);
+ strcpy(line,stra);
+
+ for (j=ncovcol;j>=1;j--){
cutv(stra, strb,line,' ');
errno=0;
lval=strtol(strb,&endptr,10);
if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight. Exiting.\n",lval, i,line,linei);
+ printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line);
exit(1);
}
- weight[i]=(double)(lval);
- strcpy(line,stra);
-
- for (j=ncovcol;j>=1;j--){
- cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, i,line,linei);
- exit(1);
- }
- if(lval <0 || lval >1){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,i,line,linei,j);
- exit(1);
- }
- covar[j][i]=(double)(lval);
- strcpy(line,stra);
- }
- lstra=strlen(stra);
-
- if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
- stratrunc = &(stra[lstra-9]);
- num[i]=atol(stratrunc);
+ if(lval <-1 || lval >1){
+ printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,linei, i,line,j);
+ exit(1);
}
- else
- num[i]=atol(stra);
- printf ("num [i] %ld %d\n",i, num[i]);fflush(stdout);
- /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
- printf("%ld %.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;
+ covar[j][i]=(double)(lval);
+ strcpy(line,stra);
+ }
+ lstra=strlen(stra);
+
+ if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
+ stratrunc = &(stra[lstra-9]);
+ num[i]=atol(stratrunc);
+ }
+ else
+ num[i]=atol(stra);
+ /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
+ printf("%ld %.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;
} /* End loop reading data */
+ fclose(fic);
/* printf("ii=%d", ij);
scanf("%d",i);*/
imx=i-1; /* Number of individuals */
@@ -4714,8 +4765,7 @@ int main(int argc, char *argv[])
printf("cptcovprod=%d ", cptcovprod);
fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
- scanf("%d ",i);
- fclose(fic);*/
+ scanf("%d ",i);*/
/* if(mle==1){*/
if (weightopt != 1) { /* Maximisation without weights*/
@@ -4940,6 +4990,7 @@ Interval (in months) between two waves:
p=param[1][1]; /* *(*(*(param +1)+1)+0) */
globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
+
if (mle==-3){
ximort=matrix(1,NDIM,1,NDIM);
cens=ivector(1,n);
@@ -4949,9 +5000,9 @@ Interval (in months) between two waves:
for (i=1; i<=imx; i++){
dcwave[i]=-1;
- for (j=1; j<=lastpass; j++)
- if (s[j][i]>nlstate) {
- dcwave[i]=j;
+ for (m=firstpass; m<=lastpass; m++)
+ if (s[m][i]>nlstate) {
+ dcwave[i]=m;
/* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
break;
}
@@ -4960,12 +5011,16 @@ Interval (in months) between two waves:
for (i=1; i<=imx; i++) {
if (wav[i]>0){
ageexmed[i]=agev[mw[1][i]][i];
- j=wav[i];agecens[i]=1.;
- if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i];
- cens[i]=1;
-
- if (ageexmed[i]<1) cens[i]=-1;
- if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ;
+ j=wav[i];
+ agecens[i]=1.;
+
+ if (ageexmed[i]> 1 && wav[i] > 0){
+ agecens[i]=agev[mw[j][i]][i];
+ cens[i]= 1;
+ }else if (ageexmed[i]< 1)
+ cens[i]= -1;
+ if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
+ cens[i]=0 ;
}
else cens[i]=-1;
}
@@ -4974,29 +5029,29 @@ Interval (in months) between two waves:
for (j=1;j<=NDIM;j++)
ximort[i][j]=(i == j ? 1.0 : 0.0);
}
-
- p[1]=0.1; p[2]=0.1;
+
+ p[1]=0.0268; p[NDIM]=0.083;
/*printf("%lf %lf", p[1], p[2]);*/
- printf("Powell\n"); fprintf(ficlog,"Powell\n");
- strcpy(filerespow,"pow-mort");
- strcat(filerespow,fileres);
- if((ficrespow=fopen(filerespow,"w"))==NULL) {
- printf("Problem with resultfile: %s\n", filerespow);
- fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
- }
- fprintf(ficrespow,"# Powell\n# iter -2*LL");
- /* for (i=1;i<=nlstate;i++)
- for(j=1;j<=nlstate+ndeath;j++)
- if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
- */
- fprintf(ficrespow,"\n");
-
+ printf("Powell\n"); fprintf(ficlog,"Powell\n");
+ strcpy(filerespow,"pow-mort");
+ strcat(filerespow,fileres);
+ if((ficrespow=fopen(filerespow,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", filerespow);
+ fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
+ }
+ fprintf(ficrespow,"# Powell\n# iter -2*LL");
+ /* for (i=1;i<=nlstate;i++)
+ for(j=1;j<=nlstate+ndeath;j++)
+ if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
+ */
+ fprintf(ficrespow,"\n");
+
powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
fclose(ficrespow);
- hesscov(matcov, p, NDIM,delti, 1e-4, gompertz);
+ hesscov(matcov, p, NDIM, delti, 1e-4, gompertz);
for(i=1; i <=NDIM; i++)
for(j=i+1;j<=NDIM;j++)
@@ -5014,48 +5069,48 @@ Interval (in months) between two waves:
for (i=1;i<=NDIM;i++)
printf("%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(1,AGESUP);
lpop=vector(1,AGESUP);
tpop=vector(1,AGESUP);
lsurv[agegomp]=100000;
-
- for (k=agegomp;k<=AGESUP;k++) {
+
+ for (k=agegomp;k<=AGESUP;k++) {
agemortsup=k;
if (p[1]*exp(p[2]*(k-agegomp))>1) break;
}
-
- for (k=agegomp;k=1 */
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
@@ -5594,11 +5649,11 @@ lsurv=vector(1,AGESUP);
free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
} /* mle==-3 arrives here for freeing */
+ free_matrix(prlim,1,nlstate,1,nlstate);
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(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
-
free_matrix(covar,0,NCOVMAX,1,n);
free_matrix(matcov,1,npar,1,npar);
/*free_vector(delti,1,npar);*/
@@ -5613,7 +5668,8 @@ lsurv=vector(1,AGESUP);
free_ivector(Tage,1,15);
free_ivector(Tcode,1,100);
-
+ free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
+ free_imatrix(codtab,1,100,1,10);
fflush(fichtm);
fflush(ficgp);
@@ -5648,32 +5704,29 @@ lsurv=vector(1,AGESUP);
/*------ End -----------*/
chdir(path);
-#ifndef UNIX
- /* strcpy(plotcmd,"\""); */
-#endif
- strcpy(plotcmd,pathimach);
/*strcat(plotcmd,CHARSEPARATOR);*/
- strcat(plotcmd,GNUPLOTPROGRAM);
+ sprintf(plotcmd,"gnuplot");
#ifndef UNIX
- strcat(plotcmd,".exe");
- /* strcat(plotcmd,"\"");*/
+ sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);
#endif
- if(stat(plotcmd,&info)){
+ if(!stat(plotcmd,&info)){
printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
- }
-
-#ifndef UNIX
- strcpy(plotcmd,"\"");
-#endif
- strcat(plotcmd,pathimach);
- strcat(plotcmd,GNUPLOTPROGRAM);
-#ifndef UNIX
- strcat(plotcmd,".exe");
- strcat(plotcmd,"\"");
+ if(!stat(getenv("GNUPLOTBIN"),&info)){
+ printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);
+ }else
+ strcpy(pplotcmd,plotcmd);
+#ifdef UNIX
+ strcpy(plotcmd,GNUPLOTPROGRAM);
+ if(!stat(plotcmd,&info)){
+ printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
+ }else
+ strcpy(pplotcmd,plotcmd);
#endif
- strcat(plotcmd," ");
- strcat(plotcmd,optionfilegnuplot);
- printf("Starting graphs with: %s",plotcmd);fflush(stdout);
+ }else
+ strcpy(pplotcmd,plotcmd);
+
+ sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot);
+ printf("Starting graphs with: %s\n",plotcmd);fflush(stdout);
if((outcmd=system(plotcmd)) != 0){
printf("\n Problem with gnuplot\n");