[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.246 and 1.251

version 1.246, 2016/09/02 08:49:22 version 1.251, 2016/09/15 15:01:13
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.251  2016/09/15 15:01:13  brouard
     Summary: not working
   
     Revision 1.250  2016/09/08 16:07:27  brouard
     Summary: continue
   
     Revision 1.249  2016/09/07 17:14:18  brouard
     Summary: Starting values from frequencies
   
     Revision 1.248  2016/09/07 14:10:18  brouard
     *** empty log message ***
   
     Revision 1.247  2016/09/02 11:11:21  brouard
     *** empty log message ***
   
   Revision 1.246  2016/09/02 08:49:22  brouard    Revision 1.246  2016/09/02 08:49:22  brouard
   *** empty log message ***    *** empty log message ***
   
Line 903  typedef struct { Line 918  typedef struct {
 /* #include <libintl.h> */  /* #include <libintl.h> */
 /* #define _(String) gettext (String) */  /* #define _(String) gettext (String) */
   
 #define MAXLINE 1024 /* Was 256. Overflow with 312 with 2 states and 4 covariates. Should be ok */  #define MAXLINE 2048 /* Was 256 and 1024. Overflow with 312 with 2 states and 4 covariates. Should be ok */
   
 #define GNUPLOTPROGRAM "gnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
Line 2584  Earliest age to start was %d-%d=%d, ncvl Line 2599  Earliest age to start was %d-%d=%d, ncvl
   /* If we start from prlim again, prlim tends to a constant matrix */    /* If we start from prlim again, prlim tends to a constant matrix */
   
   int i, ii,j,k;    int i, ii,j,k;
     int first=0;
   double *min, *max, *meandiff, maxmax,sumnew=0.;    double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **bmij();    double **out, cov[NCOVMAX+1], **bmij();
Line 2709  Earliest age to start was %d-%d=%d, ncvl Line 2725  Earliest age to start was %d-%d=%d, ncvl
     }      }
   } /* age loop */    } /* age loop */
     /* After some age loop it doesn't converge */      /* After some age loop it doesn't converge */
   printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\    if(first){
       first=1;
       printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\
   Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
     }
     fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
 Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);  Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
   /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */    /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
   free_vector(min,1,nlstate);    free_vector(min,1,nlstate);
Line 3596  double funcone( double *x) Line 3617  double funcone( double *x)
       agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */        agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
       ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */        ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
       for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */        for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */
         /* for(d=0; d<=0; d++){  /\* Delay between two effective waves Only one matrix to speed up*\/ */
         /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]          /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
           and mw[mi+1][i]. dh depends on stepm.*/            and mw[mi+1][i]. dh depends on stepm.*/
         newm=savm;          newm=savm;
         agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;  /* Here d is needed */
         cov[2]=agexact;          cov[2]=agexact;
         if(nagesqr==1)          if(nagesqr==1)
           cov[3]= agexact*agexact;            cov[3]= agexact*agexact;
Line 4167  void pstamp(FILE *fichier) Line 4189  void pstamp(FILE *fichier)
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \  void  freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
                   int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \                    int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \
                   int firstpass,  int lastpass, int stepm, int weightopt, char model[])                    int firstpass,  int lastpass, int stepm, int weightopt, char model[])
 {  /* Some frequencies */  {  /* Some frequencies as well as proposing some starting values */
       
   int i, m, jk, j1, bool, z1,j, k, iv;    int i, m, jk, j1, bool, z1,j, k, iv, jj=0;
   int iind=0, iage=0;    int iind=0, iage=0;
   int mi; /* Effective wave */    int mi; /* Effective wave */
   int first;    int first;
Line 4185  void  freqsummary(char fileres[], int ia Line 4207  void  freqsummary(char fileres[], int ia
   double agebegin, ageend;    double agebegin, ageend;
           
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
   prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);     prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); 
   posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */     posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ 
   pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */     pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ 
   /* prop=matrix(1,nlstate,iagemin,iagemax+3); */    /* prop=matrix(1,nlstate,iagemin,iagemax+3); */
Line 4229  Title=%s <br>Datafile=%s Firstpass=%d La Line 4251  Title=%s <br>Datafile=%s Firstpass=%d La
   }    }
   fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies 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) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);    fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies 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) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
       
   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);    freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+4+AGEMARGE);
   j1=0;    j1=0;
       
   /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */    /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */
   j=cptcoveff;  /* Only dummy covariates of the model */    j=cptcoveff;  /* Only dummy covariates of the model */
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    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:    /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
      reference=low_education V1=0,V2=0       reference=low_education V1=0,V2=0
Line 4244  Title=%s <br>Datafile=%s Firstpass=%d La Line 4265  Title=%s <br>Datafile=%s Firstpass=%d La
      high_educ               V1=0 V2=1       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        Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff 
   */    */
       dateintsum=0;
   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 */    k2cpt=0;
     posproptt=0.;  
     /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);    for (j = 0; j <= cptcoveff; j+=cptcoveff){   /* j= 0 constant model */
       scanf("%d", i);*/      first=1;
     for (i=-5; i<=nlstate+ndeath; i++)        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 */
       for (jk=-5; jk<=nlstate+ndeath; jk++)          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 (i=1; i<=nlstate; i++)  {
         for(m=iagemin; m <= iagemax+3; m++)          for(m=iagemin; m <= iagemax+3; m++)
           freq[i][jk][m]=0;            prop[i][m]=0;
               posprop[i]=0;
     for (i=1; i<=nlstate; i++)  {          pospropt[i]=0;
       for(m=iagemin; m <= iagemax+3; m++)        }
         prop[i][m]=0;        /* for (z1=1; z1<= nqfveff; z1++) {   */
       posprop[i]=0;        /*   meanq[z1]+=0.; */
       pospropt[i]=0;        /*   for(m=1;m<=lastpass;m++){ */
     }        /*        meanqt[m][z1]=0.; */
     /* for (z1=1; z1<= nqfveff; z1++) {   */        /*   } */
     /*   meanq[z1]+=0.; */        /* } */
     /*   for(m=1;m<=lastpass;m++){ */        
     /*  meanqt[m][z1]=0.; */        /* dateintsum=0; */
     /*   } */        /* k2cpt=0; */
     /* } */        
             /* For that combination of covariate j1, we count and print the frequencies in one pass */
     dateintsum=0;        for (iind=1; iind<=imx; iind++) { /* For each individual iind */
     k2cpt=0;          bool=1;
     /* For that combination of covariate j1, we count and print the frequencies in one pass */          if(j !=0){
     for (iind=1; iind<=imx; iind++) { /* For each individual iind */            if(anyvaryingduminmodel==0){ /* If All fixed covariates */
       bool=1;              if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
       if(anyvaryingduminmodel==0){ /* If All fixed covariates */                /* for (z1=1; z1<= nqfveff; z1++) {   */
         if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */                /*   meanq[z1]+=coqvar[Tvar[z1]][iind];  /\* Computes mean of quantitative with selected filter *\/ */
           /* for (z1=1; z1<= nqfveff; z1++) {   */                /* } */
           /*   meanq[z1]+=coqvar[Tvar[z1]][iind];  /\* Computes mean of quantitative with selected filter *\/ */                for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */
           /* } */                  /* if(Tvaraff[z1] ==-20){ */
           for (z1=1; z1<=cptcoveff; z1++) {                    /*       /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
             /* if(Tvaraff[z1] ==-20){ */                  /* }else  if(Tvaraff[z1] ==-10){ */
             /*   /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */                  /*       /\* sumnew+=coqvar[z1][iind]; *\/ */
             /* }else  if(Tvaraff[z1] ==-10){ */                  /* }else  */
             /*   /\* sumnew+=coqvar[z1][iind]; *\/ */                  if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */
             /* }else  */                    /* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */
             if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){                    bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
               /* Tests if this individual iind responded to j1 (V4=1 V3=0) */                    /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
               bool=0;                       bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
               /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",                        j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
                  bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),                    /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
                  j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/                  } /* Onlyf fixed */
               /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/                } /* end z1 */
             } /* Onlyf fixed */              } /* cptcovn > 0 */
           } /* end z1 */            } /* end any */
         } /* cptcovn > 0 */          }/* end j==0 */
       } /* end any */          if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */
       if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */            /* for(m=firstpass; m<=lastpass; m++){ */
         /* for(m=firstpass; m<=lastpass; m++){ */            for(mi=1; mi<wav[iind];mi++){ /* For that wave */
         for(mi=1; mi<wav[iind];mi++){ /* For that wave */              m=mw[mi][iind];
           m=mw[mi][iind];              if(j!=0){
           if(anyvaryingduminmodel==1){ /* Some are varying covariates */                if(anyvaryingduminmodel==1){ /* Some are varying covariates */
             for (z1=1; z1<=cptcoveff; z1++) {                  for (z1=1; z1<=cptcoveff; z1++) {
               if( Fixed[Tmodelind[z1]]==1){                    if( Fixed[Tmodelind[z1]]==1){
                 iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;                      iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
                 if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */                      if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality. If covariate's 
                   bool=0;                                                                                        value is -1, we don't select. It differs from the 
               }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */                                                                                        constant and age model which counts them. */
                 if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {                        bool=0; /* not selected */
                   bool=0;                    }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
                       if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
                         bool=0;
                       }
                     }
                 }                  }
                 }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */
               } /* end j==0 */
               /* bool =0 we keep that guy which corresponds to the combination of dummy values */
               if(bool==1){
                 /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
                    and mw[mi+1][iind]. dh depends on stepm. */
                 agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
                 ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
                 if(m >=firstpass && m <=lastpass){
                   k2=anint[m][iind]+(mint[m][iind]/12.);
                   /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
                   if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
                   if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
                   if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */
                     prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */
                   if (m<lastpass) {
                     /* if(s[m][iind]==4 && s[m+1][iind]==4) */
                     /*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
                     if(s[m][iind]==-1)
                       printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
                     freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
                     /* if((int)agev[m][iind] == 55) */
                     /*   printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */
                     /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
                     freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
                   }
                 } /* end if between passes */  
                 if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99) && (j==0)) {
                   dateintsum=dateintsum+k2; /* on all covariates ?*/
                   k2cpt++;
                   /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
               }                }
             }              }else{
           }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */                bool=1;
           /* bool =0 we keep that guy which corresponds to the combination of dummy values */              }/* end bool 2 */
           if(bool==1){            } /* end m */
             /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]          } /* end bool */
                and mw[mi+1][iind]. dh depends on stepm. */        } /* end iind = 1 to imx */
             agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/        /* prop[s][age] is feeded for any initial and valid live state as well as
             ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */           freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
             if(m >=firstpass && m <=lastpass){        
               k2=anint[m][iind]+(mint[m][iind]/12.);        
               /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/        /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
               if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */        pstamp(ficresp);
               if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */        if  (cptcoveff>0 && j!=0){
               if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */          printf( "\n#********** Variable "); 
                 prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */          fprintf(ficresp, "\n#********** Variable "); 
               if (m<lastpass) {          fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
                 /* if(s[m][iind]==4 && s[m+1][iind]==4) */          fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
                 /*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */          fprintf(ficlog, "\n#********** Variable "); 
                 if(s[m][iind]==-1)          for (z1=1; z1<=cptcoveff; z1++){
                   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));            if(!FixedV[Tvaraff[z1]]){
                 freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */              printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
                 /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */              fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
                 freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */              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)]);
             } /* end if between passes */                fprintf(ficlog, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
             if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {            }else{
               dateintsum=dateintsum+k2;              printf( "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
               k2cpt++;              fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
               /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */              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)]);
           } /* end bool 2 */              fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
         } /* end m */            }
       } /* end bool */          }
     } /* end iind = 1 to imx */          printf( "**********\n#");
     /* prop[s][age] is feeded for any initial and valid live state as well as          fprintf(ficresp, "**********\n#");
        freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */          fprintf(ficresphtm, "**********</h3>\n");
               fprintf(ficresphtmfr, "**********</h3>\n");
               fprintf(ficlog, "**********\n");
     /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/        }
     pstamp(ficresp);        fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
     if  (cptcoveff>0){        for(i=1; i<=nlstate;i++) {
       fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, " Age Prev(%d)  N(%d)  N  ",i,i);
       fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");           fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
       fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");         }
       fprintf(ficlog, "\n#********** Variable ");         fprintf(ficresp, "\n");
       for (z1=1; z1<=cptcoveff; z1++){        fprintf(ficresphtm, "\n");
         if(DummyV[z1]){        
           fprintf(ficresp, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);        /* Header of frequency table by age */
           fprintf(ficresphtm, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);        fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
           fprintf(ficresphtmfr, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);        fprintf(ficresphtmfr,"<th>Age</th> ");
           fprintf(ficlog, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);        for(jk=-1; jk <=nlstate+ndeath; jk++){
         }else{          for(m=-1; m <=nlstate+ndeath; m++){
           fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);            if(jk!=0 && m!=0)
           fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);              fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
           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, "**********</h3>\n");  
       fprintf(ficresphtmfr, "**********</h3>\n");  
       fprintf(ficlog, "**********\n");  
     }  
     fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");  
     for(i=1; i<=nlstate;i++) {  
       fprintf(ficresp, " Age Prev(%d)  N(%d)  N  ",i,i);  
       fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);  
     }  
     fprintf(ficresp, "\n");  
     fprintf(ficresphtm, "\n");  
       
     /* Header of frequency table by age */  
     fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");  
     fprintf(ficresphtmfr,"<th>Age</th> ");  
     for(jk=-1; jk <=nlstate+ndeath; jk++){  
       for(m=-1; m <=nlstate+ndeath; m++){  
         if(jk!=0 && m!=0)  
           fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);  
       }  
     }  
     fprintf(ficresphtmfr, "\n");  
       
     /* For each age */  
     for(iage=iagemin; iage <= iagemax+3; iage++){  
       fprintf(ficresphtm,"<tr>");  
       if(iage==iagemax+1){  
         fprintf(ficlog,"1");  
         fprintf(ficresphtmfr,"<tr><th>0</th> ");  
       }else if(iage==iagemax+2){  
         fprintf(ficlog,"0");  
         fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");  
       }else if(iage==iagemax+3){  
         fprintf(ficlog,"Total");  
         fprintf(ficresphtmfr,"<tr><th>Total</th> ");  
       }else{  
         if(first==1){  
           first=0;  
           printf("See log file for details...\n");  
         }          }
         fprintf(ficresphtmfr,"<tr><th>%d</th> ",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(jk=1; jk <=nlstate ; jk++){        fprintf(ficresphtmfr, "\n");
         for(m=-1, pos=0; m <=0 ; m++)      
           pos += freq[jk][m][iage];        /* For each age */
         if(pp[jk]>=1.e-10){        for(iage=iagemin; iage <= iagemax+3; iage++){
           fprintf(ficresphtm,"<tr>");
           if(iage==iagemax+1){
             fprintf(ficlog,"1");
             fprintf(ficresphtmfr,"<tr><th>0</th> ");
           }else if(iage==iagemax+2){
             fprintf(ficlog,"0");
             fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
           }else if(iage==iagemax+3){
             fprintf(ficlog,"Total");
             fprintf(ficresphtmfr,"<tr><th>Total</th> ");
           }else{
           if(first==1){            if(first==1){
             printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);              first=0;
               printf("See log file for details...\n");
             }
             fprintf(ficresphtmfr,"<tr><th>%d</th> ",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(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);
           }            }
           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];  
       } /* 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++){          for(jk=1; jk <=nlstate ; jk++){ 
         pos += pp[jk]; /* pos is the total number of transitions until this age */            /* posprop[jk]=0; */
         posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state            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] */                                            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){          for(jk=1; jk <=nlstate ; jk++){
           if(pos>=1.e-5){            if(pos>=1.e-5){
             fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);              if(first==1)
             fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);                printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
             /*probs[iage][jk][j1]= pp[jk]/pos;*/              fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*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{
           }              if(first==1)
           else{                printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
             fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);              fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
             fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);            }
             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,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",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,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
               }
           }            }
         }            pospropt[jk] +=posprop[jk];
         pospropt[jk] +=posprop[jk];          } /* end loop jk */
       } /* end loop jk */          /* pospropt=0.; */
       /* pospropt=0.; */          for(jk=-1; jk <=nlstate+ndeath; jk++){
       for(jk=-1; jk <=nlstate+ndeath; jk++){            for(m=-1; m <=nlstate+ndeath; m++){
         for(m=-1; m <=nlstate+ndeath; m++){              if(freq[jk][m][iage] !=0 ) { /* minimizing output */
           if(freq[jk][m][iage] !=0 ) { /* minimizing output */                if(first==1){
             if(first==1){                  printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
               printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);                }
                 fprintf(ficlog," %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,"<td>%.0f</td> ",freq[jk][m][iage]);
           }            }
           if(jk!=0 && m!=0)          } /* end loop jk */
             fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);          posproptt=0.; 
           for(jk=1; jk <=nlstate; jk++){
             posproptt += pospropt[jk];
           }
           fprintf(ficresphtmfr,"</tr>\n ");
           if(iage <= iagemax){
             fprintf(ficresp,"\n");
             fprintf(ficresphtm,"</tr>\n");
           }
           if(first==1)
             printf("Others in log...\n");
           fprintf(ficlog,"\n");
         } /* end loop age iage */
         fprintf(ficresphtm,"<tr><th>Tot</th>");
         for(jk=1; jk <=nlstate ; jk++){
           if(posproptt < 1.e-5){
             fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt); 
           }else{
             fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);  
         }          }
       } /* end loop jk */  
       posproptt=0.;   
       for(jk=1; jk <=nlstate; jk++){  
         posproptt += pospropt[jk];  
       }  
       fprintf(ficresphtmfr,"</tr>\n ");  
       if(iage <= iagemax){  
         fprintf(ficresp,"\n");  
         fprintf(ficresphtm,"</tr>\n");  
       }        }
       if(first==1)        fprintf(ficresphtm,"</tr>\n");
         printf("Others in log...\n");        fprintf(ficresphtm,"</table>\n");
       fprintf(ficlog,"\n");        fprintf(ficresphtmfr,"</table>\n");
     } /* end loop age iage */  
     fprintf(ficresphtm,"<tr><th>Tot</th>");  
     for(jk=1; jk <=nlstate ; jk++){  
       if(posproptt < 1.e-5){        if(posproptt < 1.e-5){
         fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);             fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
           fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
           fprintf(ficres,"\n  This combination (%d) is not valid and no result will be produced\n\n",j1);
           invalidvarcomb[j1]=1;
       }else{        }else{
         fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);              fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);
           invalidvarcomb[j1]=0;
       }        }
     }        fprintf(ficresphtmfr,"</table>\n");
     fprintf(ficresphtm,"</tr>\n");        fprintf(ficlog,"\n");
     fprintf(ficresphtm,"</table>\n");        if(j!=0){
     fprintf(ficresphtmfr,"</table>\n");          printf("#Freqsummary: Starting values for combination j1=%d:\n", j1);
     if(posproptt < 1.e-5){          for(i=1,jk=1; i <=nlstate; i++){
       fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);            for(k=1; k <=(nlstate+ndeath); k++){
       fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);              if (k != i) {
       fprintf(ficres,"\n  This combination (%d) is not valid and no result will be produced\n\n",j1);                printf("%d%d ",i,k);
       invalidvarcomb[j1]=1;                fprintf(ficlog,"%d%d ",i,k);
     }else{                for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */
       fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);                  if(jj==1){  /* Constant case */
       invalidvarcomb[j1]=0;                    if(j1==1){ /* All dummy covariates to zero */
     }                      freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */
     fprintf(ficresphtmfr,"</table>\n");                      freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */
   } /* end selected combination of covariate j1 */                    }
                     printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]));
                     fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
                     pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
                   }else if( (log(j1-1)/log(2)+1 == jj -2 -nagesqr)  && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ 
                     pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]));
                     printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4]));
                   }else if(jj==2 || nagesqr==1){ /* age or age*age parameter */
                     ;
                   }else{ /* Other cases, like quantitative fixed or varying covariates */
                     ;
                   }
                   /* printf("%12.7f )", param[i][jj][k]); */
                   /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
                   jk++; 
                 } /* end jj */
                 printf("\n");
                 fprintf(ficlog,"\n");
               } /* end k!= i */
             } /* end k */
           } /* end i, jk */
         } /* end j !=0 */
       } /* end selected combination of covariate j1 */
       if(j==0){ /* We can estimate starting values from the occurences in each case */
         printf("#Freqsummary: Starting values for the constants:\n");
         fprintf(ficlog,"\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);
               for(jj=1; jj <=ncovmodel; jj++){
                 if(jj==1){
                   printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
                   fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
                 }
                 /* printf("%12.7f )", param[i][jj][k]); */
                 /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
                 jk++; 
               }
               printf("\n");
               fprintf(ficlog,"\n");
             }
           }
         }
         printf("#Freqsummary\n");
         fprintf(ficlog,"\n");
         for(jk=-1; jk <=nlstate+ndeath; jk++){
           for(m=-1; m <=nlstate+ndeath; m++){
             /* param[i]|j][k]= freq[jk][m][iagemax+3] */
             printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
             fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
             /* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */
             /*   printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
             /*   fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
             /* } */
           }
         } /* end loop jk */
         
         printf("\n");
         fprintf(ficlog,"\n");
       } /* end j=0 */
     } /* end j */
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
       
   fclose(ficresp);    fclose(ficresp);
Line 4522  Title=%s <br>Datafile=%s Firstpass=%d La Line 4637  Title=%s <br>Datafile=%s Firstpass=%d La
   fclose(ficresphtmfr);    fclose(ficresphtmfr);
   free_vector(meanq,1,nqfveff);    free_vector(meanq,1,nqfveff);
   free_matrix(meanqt,1,lastpass,1,nqtveff);    free_matrix(meanqt,1,lastpass,1,nqtveff);
   free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE);    free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_vector(pospropt,1,nlstate);    free_vector(pospropt,1,nlstate);
   free_vector(posprop,1,nlstate);    free_vector(posprop,1,nlstate);
   free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);    free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
   /* End of freqsummary */    /* End of freqsummary */
 }  }
Line 4552  void prevalence(double ***probs, double Line 4667  void prevalence(double ***probs, double
   iagemin= (int) agemin;    iagemin= (int) agemin;
   iagemax= (int) agemax;    iagemax= (int) agemax;
   /*pp=vector(1,nlstate);*/    /*pp=vector(1,nlstate);*/
   prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);     prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); 
   /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/    /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
   j1=0;    j1=0;
       
Line 4562  void prevalence(double ***probs, double Line 4677  void prevalence(double ***probs, double
   first=1;    first=1;
   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */    for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
     for (i=1; i<=nlstate; i++)        for (i=1; i<=nlstate; i++)  
       for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)        for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++)
         prop[i][iage]=0.0;          prop[i][iage]=0.0;
     printf("Prevalence combination of varying and fixed dummies %d\n",j1);      printf("Prevalence combination of varying and fixed dummies %d\n",j1);
     /* fprintf(ficlog," V%d=%d ",Tvaraff[j1],nbcode[Tvaraff[j1]][codtabm(k,j1)]); */      /* fprintf(ficlog," V%d=%d ",Tvaraff[j1],nbcode[Tvaraff[j1]][codtabm(k,j1)]); */
Line 4593  void prevalence(double ***probs, double Line 4708  void prevalence(double ***probs, double
             if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */              if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
               if(agev[m][i]==0) agev[m][i]=iagemax+1;                if(agev[m][i]==0) agev[m][i]=iagemax+1;
               if(agev[m][i]==1) agev[m][i]=iagemax+2;                if(agev[m][i]==1) agev[m][i]=iagemax+2;
               if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){                if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+4+AGEMARGE){
                 printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d  m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m);                   printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d  m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); 
                 exit(1);                  exit(1);
               }                }
Line 4630  void prevalence(double ***probs, double Line 4745  void prevalence(double ***probs, double
       
   /*  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/    /*  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
   /*free_vector(pp,1,nlstate);*/    /*free_vector(pp,1,nlstate);*/
   free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE);    free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+4+AGEMARGE);
 }  /* End of prevalence */  }  /* End of prevalence */
   
 /************* Waves Concatenation ***************/  /************* Waves Concatenation ***************/
Line 8833  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8948  Dummy[k] 0=dummy (0 1), 1 quantitative (
 }  }
   
 int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )  int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
 {  {/* Check ages at death */
   int i, m;    int i, m;
   int firstone=0;    int firstone=0;
       
Line 9566  int main(int argc, char *argv[]) Line 9681  int main(int argc, char *argv[])
   double **prlim;    double **prlim;
   double **bprlim;    double **bprlim;
   double ***param; /* Matrix of parameters */    double ***param; /* Matrix of parameters */
   double  *p;    double ***paramstart; /* Matrix of starting parameter values */
     double  *p, *pstart; /* p=param[1][1] pstart is for starting values guessed by freqsummary */
   double **matcov; /* Matrix of covariance */    double **matcov; /* Matrix of covariance */
   double **hess; /* Hessian matrix */    double **hess; /* Hessian matrix */
   double ***delti3; /* Scale */    double ***delti3; /* Scale */
Line 9877  int main(int argc, char *argv[]) Line 9993  int main(int argc, char *argv[])
   delti=delti3[1][1];    delti=delti3[1][1];
   /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/    /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
   if(mle==-1){ /* Print a wizard for help writing covariance matrix */    if(mle==-1){ /* Print a wizard for help writing covariance matrix */
   /* We could also provide initial parameters values giving by simple logistic regression 
    * only one way, that is without matrix product. We will have nlstate maximizations */
         /* for(i=1;i<nlstate;i++){ */
         /*        /\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
         /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
         /* } */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);      prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);      printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
     fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);      fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
Line 9885  int main(int argc, char *argv[]) Line 10007  int main(int argc, char *argv[])
     fclose (ficlog);      fclose (ficlog);
     goto end;      goto end;
     exit(0);      exit(0);
     } else if(mle==-2) { /* Guessing from means */
       prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
       printf(" You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
       fprintf(ficlog," You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
      
   }  else if(mle==-5) { /* Main Wizard */    }  else if(mle==-5) { /* Main Wizard */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);      prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);      printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
Line 9906  int main(int argc, char *argv[]) Line 10033  int main(int argc, char *argv[])
     ungetc(c,ficpar);      ungetc(c,ficpar);
           
     param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);      param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
       paramstart= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
     for(i=1; i <=nlstate; i++){      for(i=1; i <=nlstate; i++){
       j=0;        j=0;
       for(jj=1; jj <=nlstate+ndeath; jj++){        for(jj=1; jj <=nlstate+ndeath; jj++){
Line 9942  run imach with mle=-1 to get a correct t Line 10070  run imach with mle=-1 to get a correct t
     }        }  
     fflush(ficlog);      fflush(ficlog);
           
     /* Reads scales values */      /* Reads parameters values */
     p=param[1][1];      p=param[1][1];
       pstart=paramstart[1][1];
           
     /* Reads comments: lines beginning with '#' */      /* Reads comments: lines beginning with '#' */
     while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
Line 10208  Please run with mle=-1 to get a correct Line 10337  Please run with mle=-1 to get a correct
   */    */
   
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);    concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   /* */    /* Concatenates waves */
     
   free_vector(moisdc,1,n);    free_vector(moisdc,1,n);
   free_vector(andc,1,n);    free_vector(andc,1,n);
Line 10389  Title=%s <br>Datafile=%s Firstpass=%d La Line 10518  Title=%s <br>Datafile=%s Firstpass=%d La
   /* Calculates basic frequencies. Computes observed prevalence at single age     /* Calculates basic frequencies. Computes observed prevalence at single age 
                  and for any valid combination of covariates                   and for any valid combination of covariates
      and prints on file fileres'p'. */       and prints on file fileres'p'. */
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \    freqsummary(fileres, p, pstart, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
               firstpass, lastpass,  stepm,  weightopt, model);                firstpass, lastpass,  stepm,  weightopt, model);
   
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
Line 10411  Interval (in months) between two waves: Line 10540  Interval (in months) between two waves:
   /* For mortality only */    /* For mortality only */
   if (mle==-3){    if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM);       ximort=matrix(1,NDIM,1,NDIM); 
                 for(i=1;i<=NDIM;i++)      for(i=1;i<=NDIM;i++)
                         for(j=1;j<=NDIM;j++)        for(j=1;j<=NDIM;j++)
                                 ximort[i][j]=0.;          ximort[i][j]=0.;
     /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */      /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
     cens=ivector(1,n);      cens=ivector(1,n);
     ageexmed=vector(1,n);      ageexmed=vector(1,n);
Line 10649  Please run with mle=-1 to get a correct Line 10778  Please run with mle=-1 to get a correct
     printf("\n");      printf("\n");
     if(mle>=1){ /* Could be 1 or 2, Real Maximization */      if(mle>=1){ /* Could be 1 or 2, Real Maximization */
       /* mlikeli uses func not funcone */        /* mlikeli uses func not funcone */
         /* for(i=1;i<nlstate;i++){ */
         /*        /\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
         /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
         /* } */
       mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);        mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
     }      }
     if(mle==0) {/* No optimization, will print the likelihoods for the datafile */      if(mle==0) {/* No optimization, will print the likelihoods for the datafile */
Removed from v.1.246  
changed lines
  Added in v.1.251

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>