--- imach/src/imach.c 2016/08/26 14:23:35 1.238
+++ imach/src/imach.c 2016/08/29 07:53:18 1.240
@@ -1,6 +1,14 @@
-/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 brouard Exp $
+/* $Id: imach.c,v 1.240 2016/08/29 07:53:18 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.240 2016/08/29 07:53:18 brouard
+ Summary: Better
+
+ Revision 1.239 2016/08/26 15:51:03 brouard
+ Summary: Improvement in Powell output in order to copy and paste
+
+ Author:
+
Revision 1.238 2016/08/26 14:23:35 brouard
Summary: Starting tests of 0.99
@@ -914,12 +922,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 brouard Exp $ */
+/* $Id: imach.c,v 1.240 2016/08/29 07:53:18 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
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 fullversion[]="$Revision: 1.238 $ $Date: 2016/08/26 14:23:35 $";
+char fullversion[]="$Revision: 1.240 $ $Date: 2016/08/29 07:53:18 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1032,7 +1040,8 @@ double dval;
#define FTOL 1.0e-10
#define NRANSI
-#define ITMAX 200
+#define ITMAX 200
+#define ITPOWMAX 20 /* This is now multiplied by the number of parameters */
#define TOL 2.0e-4
@@ -2064,7 +2073,7 @@ void powell(double p[], double **xi, int
void linmin(double p[], double xi[], int n, double *fret,
double (*func)(double []),int *flat);
#endif
- int i,ibig,j;
+ int i,ibig,j,jk,k;
double del,t,*pt,*ptt,*xit;
double directest;
double fp,fptt;
@@ -2096,13 +2105,46 @@ void powell(double p[], double **xi, int
fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
for (i=1;i<=n;i++) {
- printf(" %d %.12f",i, p[i]);
- fprintf(ficlog," %d %.12lf",i, p[i]);
fprintf(ficrespow," %.12lf", p[i]);
}
+ fprintf(ficrespow,"\n");fflush(ficrespow);
+ printf("\n#model= 1 + age ");
+ fprintf(ficlog,"\n#model= 1 + age ");
+ if(nagesqr==1){
+ printf(" + age*age ",Tvar[j]);
+ fprintf(ficlog," + age*age ",Tvar[j]);
+ }
+ for(j=1;j <=ncovmodel-2;j++){
+ if(Typevar[j]==0) {
+ printf(" + V%d ",Tvar[j]);
+ fprintf(ficlog," + V%d ",Tvar[j]);
+ }else if(Typevar[j]==1) {
+ printf(" + V%d*age ",Tvar[j]);
+ fprintf(ficlog," + V%d*age ",Tvar[j]);
+ }else if(Typevar[j]==2) {
+ printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+ fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+ }
+ }
printf("\n");
+/* printf("12 47.0114589 0.0154322 33.2424412 0.3279905 2.3731903 */
+/* 13 -21.5392400 0.1118147 1.2680506 1.2973408 -1.0663662 */
fprintf(ficlog,"\n");
- fprintf(ficrespow,"\n");fflush(ficrespow);
+ 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);
+ for(j=1; j <=ncovmodel; j++){
+ printf("%12.7f ",p[jk]);
+ fprintf(ficlog,"%12.7f ",p[jk]);
+ jk++;
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
+ }
+ }
+ }
if(*iter <=3){
tml = *localtime(&rcurr_time);
strcpy(strcurr,asctime(&tml));
@@ -2220,7 +2262,7 @@ void powell(double p[], double **xi, int
free_vector(pt,1,n);
return;
} /* enough precision */
- if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");
+ if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations.");
for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
ptt[j]=2.0*p[j]-pt[j];
xit[j]=p[j]-pt[j];
@@ -4094,7 +4136,7 @@ void freqsummary(char fileres[], int ia
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
exit(0);
}
-
+
strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
@@ -4104,55 +4146,54 @@ void freqsummary(char fileres[], int ia
}
else{
fprintf(ficresphtm,"\nIMaCh PHTM_ %s\n %s
%s \
-
\n\
+
\n \
Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\
fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
fprintf(ficresphtm,"Current page is file %s
\n\nFrequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition
\n",fileresphtm, fileresphtm);
-
+
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
fflush(ficlog);
exit(70);
- }
- else{
+ } else{
fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s \
-
\n\
+
\n \
Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\
fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
- fprintf(ficresphtmfr,"Current page is file %s
\n\nFrequencies of all effective transitions by age at begin of transition
Unknown status is -1
\n",fileresphtmfr, fileresphtmfr);
-
+ fprintf(ficresphtmfr,"Current page is file %s
\n\nFrequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate)
Unknown status is -1
\n",fileresphtmfr, fileresphtmfr);
+
freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
j1=0;
/* j=ncoveff; /\* Only fixed dummy covariates *\/ */
j=cptcoveff; /* Only dummy covariates of the model */
if (cptcovn<1) {j=1;ncodemax[1]=1;}
-
+
first=1;
-
+
/* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
reference=low_education V1=0,V2=0
med_educ V1=1 V2=0,
high_educ V1=0 V2=1
Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff
*/
-
+
for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */
posproptt=0.;
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/
for (i=-5; i<=nlstate+ndeath; i++)
for (jk=-5; jk<=nlstate+ndeath; jk++)
- for(m=iagemin; m <= iagemax+3; m++)
- freq[i][jk][m]=0;
-
+ for(m=iagemin; m <= iagemax+3; m++)
+ freq[i][jk][m]=0;
+
for (i=1; i<=nlstate; i++) {
for(m=iagemin; m <= iagemax+3; m++)
- prop[i][m]=0;
+ prop[i][m]=0;
posprop[i]=0;
pospropt[i]=0;
}
@@ -4162,7 +4203,7 @@ Title=%s
Datafile=%s Firstpass=%d La
/* meanqt[m][z1]=0.; */
/* } */
/* } */
-
+
dateintsum=0;
k2cpt=0;
/* For that combination of covariate j1, we count and print the frequencies in one pass */
@@ -4241,35 +4282,41 @@ Title=%s
Datafile=%s Firstpass=%d La
} /* end iind = 1 to imx */
/* prop[s][age] is feeded for any initial and valid live state as well as
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
-
-
+
+
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
pstamp(ficresp);
- /* if (ncoveff>0) { */
- if (cptcoveff>0) {
+ if (cptcoveff>0){
fprintf(ficresp, "\n#********** Variable ");
fprintf(ficresphtm, "\n
********** Variable ");
fprintf(ficresphtmfr, "\n
********** Variable ");
+ fprintf(ficlog, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++){
- fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ if(DummyV[z1]){
+ fprintf(ficresp, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtm, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtmfr, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficlog, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ }else{
+ fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ }
}
fprintf(ficresp, "**********\n#");
fprintf(ficresphtm, "**********
\n");
fprintf(ficresphtmfr, "**********
\n");
- fprintf(ficlog, "\n#********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficlog, "**********\n");
}
fprintf(ficresphtm,"");
for(i=1; i<=nlstate;i++) {
- fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
+ fprintf(ficresp, " Age Prev(%d) N(%d) N ",i,i);
fprintf(ficresphtm, "| Age | Prev(%d) | N(%d) | N | ",i,i);
}
fprintf(ficresp, "\n");
fprintf(ficresphtm, "\n");
-
+
/* Header of frequency table by age */
fprintf(ficresphtmfr,"");
fprintf(ficresphtmfr,"| Age | ");
@@ -4280,115 +4327,115 @@ Title=%s
Datafile=%s Firstpass=%d La
}
}
fprintf(ficresphtmfr, "\n");
-
+
/* For each age */
for(iage=iagemin; iage <= iagemax+3; iage++){
fprintf(ficresphtm,"");
if(iage==iagemax+1){
- fprintf(ficlog,"1");
- fprintf(ficresphtmfr,"
| 0 | ");
+ fprintf(ficlog,"1");
+ fprintf(ficresphtmfr,"
|---|
| 0 | ");
}else if(iage==iagemax+2){
- fprintf(ficlog,"0");
- fprintf(ficresphtmfr,"
|---|
| Unknown | ");
+ fprintf(ficlog,"0");
+ fprintf(ficresphtmfr,"
|---|
| Unknown | ");
}else if(iage==iagemax+3){
- fprintf(ficlog,"Total");
- fprintf(ficresphtmfr,"
|---|
| Total | ");
+ fprintf(ficlog,"Total");
+ fprintf(ficresphtmfr,"
|---|
| Total | ");
}else{
- if(first==1){
- first=0;
- printf("See log file for details...\n");
- }
- fprintf(ficresphtmfr,"
|---|
| %d | ",iage);
- fprintf(ficlog,"Age %d", iage);
+ if(first==1){
+ first=0;
+ printf("See log file for details...\n");
+ }
+ fprintf(ficresphtmfr,"
|---|
| %d | ",iage);
+ fprintf(ficlog,"Age %d", iage);
}
for(jk=1; jk <=nlstate ; jk++){
- for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
- pp[jk] += freq[jk][m][iage];
+ for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+ pp[jk] += freq[jk][m][iage];
}
for(jk=1; jk <=nlstate ; jk++){
- for(m=-1, pos=0; m <=0 ; m++)
- pos += freq[jk][m][iage];
- if(pp[jk]>=1.e-10){
- if(first==1){
- printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
- }
- fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
- }else{
- if(first==1)
- printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
- fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
- }
+ for(m=-1, pos=0; m <=0 ; m++)
+ pos += freq[jk][m][iage];
+ if(pp[jk]>=1.e-10){
+ if(first==1){
+ printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+ }
+ fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+ }else{
+ if(first==1)
+ printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+ fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+ }
}
-
+
for(jk=1; jk <=nlstate ; jk++){
- /* posprop[jk]=0; */
- for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
- pp[jk] += freq[jk][m][iage];
+ /* posprop[jk]=0; */
+ for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
+ pp[jk] += freq[jk][m][iage];
} /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */
-
+
for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
- pos += pp[jk]; /* pos is the total number of transitions until this age */
- posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
- from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
- pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
- from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+ pos += pp[jk]; /* pos is the total number of transitions until this age */
+ posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
+ from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+ pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
+ from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
}
for(jk=1; jk <=nlstate ; jk++){
- if(pos>=1.e-5){
- if(first==1)
- printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
- fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
- }else{
- if(first==1)
- printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
- fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
- }
- if( iage <= iagemax){
- if(pos>=1.e-5){
- fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
- fprintf(ficresphtm,"%d | %.5f | %.0f | %.0f | ",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
- /*probs[iage][jk][j1]= pp[jk]/pos;*/
- /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
- }
- else{
- fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
- fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f | ",iage, prop[jk][iage],pospropta);
- }
- }
- pospropt[jk] +=posprop[jk];
+ if(pos>=1.e-5){
+ if(first==1)
+ printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+ fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+ }else{
+ if(first==1)
+ printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+ fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+ }
+ if( iage <= iagemax){
+ if(pos>=1.e-5){
+ fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+ fprintf(ficresphtm,"%d | %.5f | %.0f | %.0f | ",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+ /*probs[iage][jk][j1]= pp[jk]/pos;*/
+ /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
+ }
+ else{
+ fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
+ fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f | ",iage, prop[jk][iage],pospropta);
+ }
+ }
+ pospropt[jk] +=posprop[jk];
} /* end loop jk */
/* pospropt=0.; */
for(jk=-1; jk <=nlstate+ndeath; jk++){
- for(m=-1; m <=nlstate+ndeath; m++){
- if(freq[jk][m][iage] !=0 ) { /* minimizing output */
- if(first==1){
- printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
- }
- fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
- }
- if(jk!=0 && m!=0)
- fprintf(ficresphtmfr,"%.0f | ",freq[jk][m][iage]);
- }
+ for(m=-1; m <=nlstate+ndeath; m++){
+ if(freq[jk][m][iage] !=0 ) { /* minimizing output */
+ if(first==1){
+ printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
+ }
+ fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
+ }
+ if(jk!=0 && m!=0)
+ fprintf(ficresphtmfr,"%.0f | ",freq[jk][m][iage]);
+ }
} /* end loop jk */
posproptt=0.;
for(jk=1; jk <=nlstate; jk++){
- posproptt += pospropt[jk];
+ posproptt += pospropt[jk];
}
fprintf(ficresphtmfr,"
|---|
\n ");
if(iage <= iagemax){
- fprintf(ficresp,"\n");
- fprintf(ficresphtm,"\n");
+ fprintf(ficresp,"\n");
+ fprintf(ficresphtm,"\n");
}
if(first==1)
- printf("Others in log...\n");
+ printf("Others in log...\n");
fprintf(ficlog,"\n");
} /* end loop age iage */
fprintf(ficresphtm,"| Tot | ");
for(jk=1; jk <=nlstate ; jk++){
if(posproptt < 1.e-5){
- fprintf(ficresphtm,"Nanq | %.0f | %.0f | ",pospropt[jk],posproptt);
+ fprintf(ficresphtm,"Nanq | %.0f | %.0f | ",pospropt[jk],posproptt);
}else{
- fprintf(ficresphtm,"%.5f | %.0f | %.0f | ",pospropt[jk]/posproptt,pospropt[jk],posproptt);
+ fprintf(ficresphtm,"%.5f | %.0f | %.0f | ",pospropt[jk]/posproptt,pospropt[jk],posproptt);
}
}
fprintf(ficresphtm,"
|---|
\n");
@@ -4406,7 +4453,7 @@ Title=%s
Datafile=%s Firstpass=%d La
fprintf(ficresphtmfr,"
\n");
} /* end selected combination of covariate j1 */
dateintmean=dateintsum/k2cpt;
-
+
fclose(ficresp);
fclose(ficresphtm);
fclose(ficresphtmfr);
@@ -6314,10 +6361,7 @@ void printinggnuplot(char fileresu[], ch
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
- fprintf(ficgp,"set xlabel \"Age\" \n\
-set ylabel \"Probability\" \n \
-set ter svg size 640, 480\n \
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
@@ -6431,7 +6475,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
/*3eme*/
for (k1=1; k1<= m ; k1 ++){
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
- if(TKresult[nres]!= k)
+ if(TKresult[nres]!= k1)
continue;
for (cpt=1; cpt<= nlstate ; cpt ++) {
@@ -7791,14 +7835,36 @@ int readdata(char datafile[], int firsto
/*-------- data file ----------*/
FILE *fic;
char dummy[]=" ";
- int i=0, j=0, n=0, iv=0;
+ int i=0, j=0, n=0, iv=0, v;
int lstra;
int linei, month, year,iout;
char line[MAXLINE], linetmp[MAXLINE];
char stra[MAXLINE], strb[MAXLINE];
char *stratrunc;
+ DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
+ FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
+ for(v=1; v <=ncovcol;v++){
+ DummyV[v]=0;
+ FixedV[v]=0;
+ }
+ for(v=ncovcol+1; v <=ncovcol+nqv;v++){
+ DummyV[v]=1;
+ FixedV[v]=0;
+ }
+ for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
+ DummyV[v]=0;
+ FixedV[v]=1;
+ }
+ for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
+ DummyV[v]=1;
+ FixedV[v]=1;
+ }
+ for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
+ printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+ fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+ }
if((fic=fopen(datafile,"r"))==NULL) {
printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
@@ -8419,27 +8485,7 @@ 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\
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);
-
- for(v=1; v <=ncovcol;v++){
- DummyV[v]=0;
- FixedV[v]=0;
- }
- for(v=ncovcol+1; v <=ncovcol+nqv;v++){
- DummyV[v]=1;
- FixedV[v]=0;
- }
- for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
- DummyV[v]=0;
- FixedV[v]=1;
- }
- for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
- DummyV[v]=1;
- FixedV[v]=1;
- }
- for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
- printf("Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
- fprintf(ficlog,"Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
- }
+ 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 */
if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
Fixed[k]= 0;
@@ -8464,7 +8510,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarFind[ncovf]=k;
TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
- }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */
+ }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */
Fixed[k]= 0;
Dummy[k]= 1;
nqfveff++;
@@ -8517,167 +8563,167 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarA[ncova]=Tvar[k];
TvarAind[ncova]=k;
if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
- Fixed[k]= 2;
- Dummy[k]= 2;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APFD;
- /* ncoveff++; */
+ Fixed[k]= 2;
+ Dummy[k]= 2;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APFD;
+ /* ncoveff++; */
}else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
- Fixed[k]= 2;
- Dummy[k]= 3;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APFQ; /* Product age * fixed quantitative */
- /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */
+ Fixed[k]= 2;
+ Dummy[k]= 3;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APFQ; /* Product age * fixed quantitative */
+ /* nqfveff++; /\* Only simple fixed quantitative variable *\/ */
}else if( Tvar[k] <=ncovcol+nqv+ntv ){
- Fixed[k]= 3;
- Dummy[k]= 2;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APVD; /* Product age * varying dummy */
- /* ntveff++; /\* Only simple time varying dummy variable *\/ */
+ Fixed[k]= 3;
+ Dummy[k]= 2;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APVD; /* Product age * varying dummy */
+ /* ntveff++; /\* Only simple time varying dummy variable *\/ */
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 3;
- Dummy[k]= 3;
- modell[k].maintype= ATYPE;
- modell[k].subtype= APVQ; /* Product age * varying quantitative */
- /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
+ Fixed[k]= 3;
+ Dummy[k]= 3;
+ modell[k].maintype= ATYPE;
+ modell[k].subtype= APVQ; /* Product age * varying quantitative */
+ /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
}
}else if (Typevar[k] == 2) { /* product without age */
k1=Tposprod[k];
if(Tvard[k1][1] <=ncovcol){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 1;
- Dummy[k]= 0;
- modell[k].maintype= FTYPE;
- modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */
- ncovf++; /* Fixed variables without age */
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv){
- Fixed[k]= 0; /* or 2 ?*/
- Dummy[k]= 1;
- modell[k].maintype= FTYPE;
- modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */
- ncovf++; /* Varying variables without age */
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 0;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ modell[k].maintype= FTYPE;
+ modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */
+ ncovf++; /* Fixed variables without age */
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 0; /* or 2 ?*/
+ Dummy[k]= 1;
+ modell[k].maintype= FTYPE;
+ modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */
+ ncovf++; /* Varying variables without age */
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else if(Tvard[k1][1] <=ncovcol+nqv){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 0; /* or 2 ?*/
- Dummy[k]= 1;
- modell[k].maintype= FTYPE;
- modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */
- ncovf++; /* Fixed variables without age */
- TvarF[ncovf]=Tvar[k];
- TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 0; /* or 2 ?*/
+ Dummy[k]= 1;
+ modell[k].maintype= FTYPE;
+ modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */
+ ncovf++; /* Fixed variables without age */
+ TvarF[ncovf]=Tvar[k];
+ TvarFind[ncovf]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 0;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 0;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
- Fixed[k]= 1;
- Dummy[k]= 1;
- modell[k].maintype= VTYPE;
- modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */
- ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }
+ if(Tvard[k1][2] <=ncovcol){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ Fixed[k]= 1;
+ Dummy[k]= 1;
+ modell[k].maintype= VTYPE;
+ modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */
+ ncovv++; /* Varying variables without age */
+ TvarV[ncovv]=Tvar[k];
+ TvarVind[ncovv]=k;
+ }
}else{
- printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
- fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
- } /* end k1 */
+ printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+ fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+ } /*end k1*/
}else{
printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
@@ -10009,8 +10055,6 @@ Please run with mle=-1 to get a correct
Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */
Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */
Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */
- DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
- FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
/* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs).
For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4,
Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
@@ -10828,10 +10872,12 @@ Please run with mle=-1 to get a correct
}else
break;
}
+ if (!feof(ficpar))
while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
- if (num_filled == 0)
+ if (num_filled == 0){
resultline[0]='\0';
- else if (num_filled != 1){
+ break;
+ } else if (num_filled != 1){
printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
}
nresult++; /* Sum of resultlines */
@@ -10861,7 +10907,7 @@ Please run with mle=-1 to get a correct
break;
else{ /* Processess output results for this combination of covariate values */
}
- }
+ } /* end while */