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

Diff for /imach/src/imach.c between versions 1.11 and 1.12

version 1.11, 2001/05/17 16:07:14 version 1.12, 2002/02/20 16:57:00
Line 22 Line 22
   delay between waves is not identical for each individual, or if some    delay between waves is not identical for each individual, or if some
   individual missed an interview, the information is not rounded or lost, but    individual missed an interview, the information is not rounded or lost, but
   taken into account using an interpolation or extrapolation.    taken into account using an interpolation or extrapolation.
   hPijx is the probability to be    hPijx is the probability to be
   observed in state i at age x+h conditional to the observed state i at age    observed in state i at age x+h conditional to the observed state i at age
   x. The delay 'h' can be split into an exact number (nh*stepm) of    x. The delay 'h' can be split into an exact number (nh*stepm) of
   unobserved intermediate  states. This elementary transition (by month or    unobserved intermediate  states. This elementary transition (by month or
Line 693  double **prevalim(double **prlim, int nl Line 693  double **prevalim(double **prlim, int nl
   }    }
 }  }
   
 /*************** transition probabilities **********/  /*************** transition probabilities ***************/
   
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {  {
Line 716  double **pmij(double **ps, double *cov, Line 716  double **pmij(double **ps, double *cov,
         s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];          s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
         /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/          /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
       }        }
       ps[i][j]=s2;        ps[i][j]=(s2);
     }      }
   }    }
       /*ps[3][2]=1;*/
   
   for(i=1; i<= nlstate; i++){    for(i=1; i<= nlstate; i++){
      s1=0;       s1=0;
     for(j=1; j<i; j++)      for(j=1; j<i; j++)
Line 740  double **pmij(double **ps, double *cov, Line 742  double **pmij(double **ps, double *cov,
     }      }
   }    }
   
   
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){    /*   for(ii=1; ii<= nlstate+ndeath; ii++){
     for(jj=1; jj<= nlstate+ndeath; jj++){      for(jj=1; jj<= nlstate+ndeath; jj++){
      printf("%lf ",ps[ii][jj]);       printf("%lf ",ps[ii][jj]);
Line 804  double ***hpxij(double ***po, int nhstep Line 807  double ***hpxij(double ***po, int nhstep
       cov[1]=1.;        cov[1]=1.;
       cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;        cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
       for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];        for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
 for (k=1; k<=cptcovage;k++)        for (k=1; k<=cptcovage;k++)
         cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];          cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
    for (k=1; k<=cptcovprod;k++)        for (k=1; k<=cptcovprod;k++)
         cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];          cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
   
   
Line 913  void hesscov(double **matcov, double p[] Line 916  void hesscov(double **matcov, double p[]
   void lubksb(double **a, int npar, int *indx, double b[]) ;    void lubksb(double **a, int npar, int *indx, double b[]) ;
   void ludcmp(double **a, int npar, int *indx, double *d) ;    void ludcmp(double **a, int npar, int *indx, double *d) ;
   
   
   hess=matrix(1,npar,1,npar);    hess=matrix(1,npar,1,npar);
   
   printf("\nCalculation of the hessian matrix. Wait...\n");    printf("\nCalculation of the hessian matrix. Wait...\n");
Line 921  void hesscov(double **matcov, double p[] Line 923  void hesscov(double **matcov, double p[]
     printf("%d",i);fflush(stdout);      printf("%d",i);fflush(stdout);
     hess[i][i]=hessii(p,ftolhess,i,delti);      hess[i][i]=hessii(p,ftolhess,i,delti);
     /*printf(" %f ",p[i]);*/      /*printf(" %f ",p[i]);*/
       /*printf(" %lf ",hess[i][i]);*/
   }    }
    
   for (i=1;i<=npar;i++) {    for (i=1;i<=npar;i++) {
     for (j=1;j<=npar;j++)  {      for (j=1;j<=npar;j++)  {
       if (j>i) {        if (j>i) {
         printf(".%d%d",i,j);fflush(stdout);          printf(".%d%d",i,j);fflush(stdout);
         hess[i][j]=hessij(p,delti,i,j);          hess[i][j]=hessij(p,delti,i,j);
         hess[j][i]=hess[i][j];          hess[j][i]=hess[i][j];    
           /*printf(" %lf ",hess[i][j]);*/
       }        }
     }      }
   }    }
Line 1032  double hessii( double x[], double delta, Line 1036  double hessii( double x[], double delta,
     }      }
   }    }
   delti[theta]=delts;    delti[theta]=delts;
   return res;    return res;
     
 }  }
   
Line 1311  void  concatwav(int wav[], int **dh, int Line 1315  void  concatwav(int wav[], int **dh, int
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           if (j <= jmin) jmin=j;            if (j <= jmin) jmin=j;
           sum=sum+j;            sum=sum+j;
           if (j<0) printf("j=%d num=%d ",j,i);            /* if (j<10) printf("j=%d num=%d ",j,i); */
           }            }
         }          }
         else{          else{
Line 1319  void  concatwav(int wav[], int **dh, int Line 1323  void  concatwav(int wav[], int **dh, int
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            else if (j <= jmin)jmin=j;
             /*   if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           sum=sum+j;            sum=sum+j;
         }          }
         jk= j/stepm;          jk= j/stepm;
Line 1335  void  concatwav(int wav[], int **dh, int Line 1340  void  concatwav(int wav[], int **dh, int
   }    }
   jmean=sum/k;    jmean=sum/k;
   printf("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 Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
 }   }
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)  void tricode(int *Tvar, int **nbcode, int imx)
 {  {
Line 1373  void tricode(int *Tvar, int **nbcode, in Line 1378  void tricode(int *Tvar, int **nbcode, in
   
  for (k=0; k<19; k++) Ndum[k]=0;   for (k=0; k<19; k++) Ndum[k]=0;
   
  for (i=1; i<=ncovmodel; i++) {   for (i=1; i<=ncovmodel-2; i++) {
       ij=Tvar[i];        ij=Tvar[i];
       Ndum[ij]++;        Ndum[ij]++;
     }      }
Line 1448  void varevsij(char fileres[], double *** Line 1453  void varevsij(char fileres[], double ***
   double **dnewm,**doldm;    double **dnewm,**doldm;
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h;
   int k, cptcode;    int k, cptcode;
    double *xp;    double *xp;
   double **gp, **gm;    double **gp, **gm;
   double ***gradg, ***trgradg;    double ***gradg, ***trgradg;
   double ***p3mat;    double ***p3mat;
Line 1663  int main() Line 1668  int main()
   int ju,jl, mi;    int ju,jl, mi;
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;    int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;    int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
     
   int hstepm, nhstepm;    int hstepm, nhstepm;
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 1772  split(pathtot, path,optionfile); Line 1777  split(pathtot, path,optionfile);
       fprintf(ficparo,"\n");        fprintf(ficparo,"\n");
     }      }
     
   npar= (nlstate+ndeath-1)*nlstate*ncovmodel;      npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
   
   p=param[1][1];    p=param[1][1];
     
   /* Reads comments: lines beginning with '#' */    /* Reads comments: lines beginning with '#' */
Line 1862  split(pathtot, path,optionfile); Line 1868  split(pathtot, path,optionfile);
     tab=ivector(1,NCOVMAX);      tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,8);      ncodemax=ivector(1,8);
   
     i=1;      i=1;
     while (fgets(line, MAXLINE, fic) != NULL)    {      while (fgets(line, MAXLINE, fic) != NULL)    {
       if ((i >= firstobs) && (i <=lastobs)) {        if ((i >= firstobs) && (i <=lastobs)) {
                 
Line 1884  split(pathtot, path,optionfile); Line 1890  split(pathtot, path,optionfile);
           cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);            cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
         }          }
         num[i]=atol(stra);          num[i]=atol(stra);
          
         /*printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/          /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
             printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
   
         i=i+1;          i=i+1;
       }        }
     }      }
       /* printf("ii=%d", ij);
     /*scanf("%d",i);*/         scanf("%d",i);*/
   imx=i-1; /* Number of individuals */    imx=i-1; /* Number of individuals */
   
     /* for (i=1; i<=imx; i++){
       if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
       if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
       if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
     }
     for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
   
   /* Calculation of the number of parameter from char model*/    /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15);    Tvar=ivector(1,15);
   Tprod=ivector(1,15);    Tprod=ivector(1,15);
Line 1975  split(pathtot, path,optionfile); Line 1989  split(pathtot, path,optionfile);
     }      }
     /*-calculation of age at interview from date of interview and age at death -*/      /*-calculation of age at interview from date of interview and age at death -*/
     agev=matrix(1,maxwav,1,imx);      agev=matrix(1,maxwav,1,imx);
   
      for (i=1; i<=imx; i++)
        for(m=2; (m<= maxwav); m++)
          if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
            anint[m][i]=9999;
            s[m][i]=-1;
          }
         
     for (i=1; i<=imx; i++)  {      for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);        agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
Line 2083  printf("Total number of individuals= %d, Line 2104  printf("Total number of individuals= %d,
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */      oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
           
     /* For Powell, parameters are in a vector p[] starting at p[1]      /* For Powell, parameters are in a vector p[] starting at p[1]
        so we point p on param[1][1] so that p[1] maps on param[1][1][1] */         so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
     p=param[1][1]; /* *(*(*(param +1)+1)+0) */      p=param[1][1]; /* *(*(*(param +1)+1)+0) */
Line 2666  strcpy(fileresvpl,"vpl"); Line 2687  strcpy(fileresvpl,"vpl");
   /*printf("Total time was %d uSec.\n", total_usecs);*/    /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/    /*------ End -----------*/
   
   
  end:   end:
 #ifdef windows  #ifdef windows
  chdir(pathcd);   chdir(pathcd);
Removed from v.1.11  
changed lines
  Added in v.1.12

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