--- imach/src/imach.c	2015/01/04 09:57:06	1.179
+++ imach/src/imach.c	2015/02/11 17:33:45	1.180
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.179 2015/01/04 09:57:06 brouard Exp $
+/* $Id: imach.c,v 1.180 2015/02/11 17:33:45 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.180  2015/02/11 17:33:45  brouard
+  Summary: Finishing move from main to function (hpijx and prevalence_limit)
+
   Revision 1.179  2015/01/04 09:57:06  brouard
   Summary: back to OS/X
 
@@ -635,11 +638,11 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.179 2015/01/04 09:57:06 brouard Exp $ */
+/* $Id: imach.c,v 1.180 2015/02/11 17:33:45 brouard Exp $ */
 /* $State: Exp $ */
 
-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 fullversion[]="$Revision: 1.179 $ $Date: 2015/01/04 09:57:06 $"; 
+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: 1.180 $ $Date: 2015/02/11 17:33:45 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -2280,9 +2283,9 @@ void mlikeli(FILE *ficres,double p[], in
 #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));
 
 }
 
@@ -5641,6 +5644,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 *****************/
 /***********************************************/
@@ -6770,7 +6910,9 @@ Interval (in months) between two waves:
 
 
     /*--------------- 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
@@ -6778,7 +6920,8 @@ Interval (in months) between two waves:
 #endif
 
     /*------------- h Pij x at various ages ------------*/
-#include "hpijx.h"
+    /*#include "hpijx.h"*/
+    hPijx(p, bage, fage);
     fclose(ficrespij);
 
   /*-------------- Variance of one-step probabilities---*/