Summary: Finishing move from main to function (hpijx and prevalence_limit)

diff --git a/src/imach.c b/src/imach.c
index 783f2e27cb7f29a713230cc68a067120aab241c1..c20f7e63664b9815b1991f9949605b32858e6ad4 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.179  2015/01/04 09:57:06  brouard
+  Summary: back to OS/X
+
   Revision 1.178  2015/01/04 09:35:48  brouard
   *** empty log message ***
 
@@ -635,7 +638,7 @@ typedef struct {
 /* $Id$ */
 /* $State$ */
 
-char version[]="Imach version 0.98p, December 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char version[]="Imach version 0.98p, Février 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
 char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
@@ -2277,9 +2280,9 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub
 #endif
   free_matrix(xi,1,npar,1,npar);
   fclose(ficrespow);
-  printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
-  fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
-  fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+  printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+  fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+  fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
 
 }
 
@@ -5638,6 +5641,143 @@ void syscompilerinfo()
 
  }
 
+int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){
+  /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
+  int i, j, k, i1 ;
+  double ftolpl = 1.e-10;
+  double age, agebase, agelim;
+
+    strcpy(filerespl,"pl");
+    strcat(filerespl,fileres);
+    if((ficrespl=fopen(filerespl,"w"))==NULL) {
+      printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+      fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+    }
+    printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+    fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+    pstamp(ficrespl);
+    fprintf(ficrespl,"# Period (stable) prevalence \n");
+    fprintf(ficrespl,"#Age ");
+    for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
+    fprintf(ficrespl,"\n");
+  
+    /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
+
+    agebase=ageminpar;
+    agelim=agemaxpar;
+
+    i1=pow(2,cptcoveff);
+    if (cptcovn < 1){i1=1;}
+
+    for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+    /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
+      //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+       k=k+1;
+       /* to clean */
+       //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]);
+       fprintf(ficrespl,"\n#******");
+       printf("\n#******");
+       fprintf(ficlog,"\n#******");
+       for(j=1;j<=cptcoveff;j++) {
+         fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+       }
+       fprintf(ficrespl,"******\n");
+       printf("******\n");
+       fprintf(ficlog,"******\n");
+
+       fprintf(ficrespl,"#Age ");
+       for(j=1;j<=cptcoveff;j++) {
+         fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+       }
+       for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
+       fprintf(ficrespl,"\n");
+       
+       for (age=agebase; age<=agelim; age++){
+       /* for (age=agebase; age<=agebase; age++){ */
+         prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+         fprintf(ficrespl,"%.0f ",age );
+         for(j=1;j<=cptcoveff;j++)
+           fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         for(i=1; i<=nlstate;i++)
+           fprintf(ficrespl," %.5f", prlim[i][i]);
+         fprintf(ficrespl,"\n");
+       } /* Age */
+       /* was end of cptcod */
+    } /* cptcov */
+}
+
+int hPijx(double *p, int bage, int fage){
+    /*------------- h Pij x at various ages ------------*/
+
+  int stepsize;
+  int agelim;
+  int hstepm;
+  int nhstepm;
+  int h, i, i1, j, k;
+
+  double agedeb;
+  double ***p3mat;
+
+    strcpy(filerespij,"pij");  strcat(filerespij,fileres);
+    if((ficrespij=fopen(filerespij,"w"))==NULL) {
+      printf("Problem with Pij resultfile: %s\n", filerespij); return 1;
+      fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1;
+    }
+    printf("Computing pij: result on file '%s' \n", filerespij);
+    fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);
+  
+    stepsize=(int) (stepm+YEARM-1)/YEARM;
+    /*if (stepm<=24) stepsize=2;*/
+
+    agelim=AGESUP;
+    hstepm=stepsize*YEARM; /* Every year of age */
+    hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
+
+    /* hstepm=1;   aff par mois*/
+    pstamp(ficrespij);
+    fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
+    i1= pow(2,cptcoveff);
+   for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+      /*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
+       k=k+1; 
+    /* for (k=1; k <= (int) pow(2,cptcoveff); k++){*/
+       fprintf(ficrespij,"\n#****** ");
+       for(j=1;j<=cptcoveff;j++) 
+         fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+       fprintf(ficrespij,"******\n");
+       
+       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
+         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
+         nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+
+         /*      nhstepm=nhstepm*YEARM; aff par mois*/
+
+         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+         oldm=oldms;savm=savms;
+         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
+         fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
+         for(i=1; i<=nlstate;i++)
+           for(j=1; j<=nlstate+ndeath;j++)
+             fprintf(ficrespij," %1d-%1d",i,j);
+         fprintf(ficrespij,"\n");
+         for (h=0; h<=nhstepm; h++){
+           /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+           fprintf(ficrespij,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm );
+           for(i=1; i<=nlstate;i++)
+             for(j=1; j<=nlstate+ndeath;j++)
+               fprintf(ficrespij," %.5f", p3mat[i][j][h]);
+           fprintf(ficrespij,"\n");
+         }
+         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+         fprintf(ficrespij,"\n");
+       }
+      /*}*/
+    }
+}
+
+
 /***********************************************/
 /**************** Main Program *****************/
 /***********************************************/
@@ -6767,7 +6907,9 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
 
 
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
-#include "prevlim.h"  /* Use ficrespl, ficlog */
+    /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */
+    prlim=matrix(1,nlstate,1,nlstate);
+    prevalence_limit(p, prlim,  ageminpar, agemaxpar);
     fclose(ficrespl);
 
 #ifdef FREEEXIT2
@@ -6775,7 +6917,8 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
 #endif
 
     /*------------- h Pij x at various ages ------------*/
-#include "hpijx.h"
+    /*#include "hpijx.h"*/
+    hPijx(p, bage, fage);
     fclose(ficrespij);
 
   /*-------------- Variance of one-step probabilities---*/