Summary: Some fixes

diff --git a/src/imach.c b/src/imach.c
index 4d16e5bec68b221640c2839eb7d03dfe73b26ec1..edd8bad1d202a8c22b81eb1e5fa45f6f15ec40f6 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,11 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.190  2015/05/05 08:51:13  brouard
+  Summary: Adding digits in output parameters (7 digits instead of 6)
+
+  Fix 1+age+.
+
   Revision 1.189  2015/04/30 14:45:16  brouard
   Summary: 0.98q2
 
@@ -1445,31 +1450,41 @@ values at the three points, fa, fb , and fc such that fa > fb and fb < fc.
 #endif 
 #ifdef MNBRAKORIGINAL
 #else
-      if (fu > *fc) {
+/*       if (fu > *fc) { */
+/* #ifdef DEBUG */
+/*       printf("mnbrak4  fu > fc \n"); */
+/*       fprintf(ficlog, "mnbrak4 fu > fc\n"); */
+/* #endif */
+/*     /\* SHFT(u,*cx,*cx,u) /\\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\\/  *\/ */
+/*     /\* SHFT(*fa,*fc,fu,*fc) /\\* (b, u, c) is a bracket while test fb > fc will be fu > fc  will exit *\\/ *\/ */
+/*     dum=u; /\* Shifting c and u *\/ */
+/*     u = *cx; */
+/*     *cx = dum; */
+/*     dum = fu; */
+/*     fu = *fc; */
+/*     *fc =dum; */
+/*       } else { /\* end *\/ */
+/* #ifdef DEBUG */
+/*       printf("mnbrak3  fu < fc \n"); */
+/*       fprintf(ficlog, "mnbrak3 fu < fc\n"); */
+/* #endif */
+/*     dum=u; /\* Shifting c and u *\/ */
+/*     u = *cx; */
+/*     *cx = dum; */
+/*     dum = fu; */
+/*     fu = *fc; */
+/*     *fc =dum; */
+/*       } */
 #ifdef DEBUG
-      printf("mnbrak4  fu > fc \n");
-      fprintf(ficlog, "mnbrak4 fu > fc\n");
+      printf("mnbrak34  fu < or >= fc \n");
+      fprintf(ficlog, "mnbrak34 fu < fc\n");
 #endif
-       /* SHFT(u,*cx,*cx,u) /\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\/  */
-       /* SHFT(*fa,*fc,fu,*fc) /\* (b, u, c) is a bracket while test fb > fc will be fu > fc  will exit *\/ */
-       dum=u; /* Shifting c and u */
-       u = *cx;
-       *cx = dum;
-       dum = fu;
-       fu = *fc;
-       *fc =dum;
-      } else { /* end */
-#ifdef DEBUG
-      printf("mnbrak3  fu < fc \n");
-      fprintf(ficlog, "mnbrak3 fu < fc\n");
-#endif
-       dum=u; /* Shifting c and u */
-       u = *cx;
-       *cx = dum;
-       dum = fu;
-       fu = *fc;
-       *fc =dum;
-      }
+      dum=u; /* Shifting c and u */
+      u = *cx;
+      *cx = dum;
+      dum = fu;
+      fu = *fc;
+      *fc =dum;
 #endif
     } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */
 #ifdef DEBUG
@@ -1560,6 +1575,9 @@ void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [
     }
   }while(fx != fx);
 
+#ifdef DEBUGLINMIN
+  printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
+#endif
   *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/
   /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */
   /* fmin = f(p[j] + xmin * xi[j]) */
@@ -1569,7 +1587,9 @@ void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [
   printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
   fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
 #endif
-  /* printf("linmin end "); */
+#ifdef DEBUGLINMIN
+  printf("linmin end ");
+#endif
   for (j=1;j<=n;j++) { 
     /* printf(" before xi[%d]=%12.8f", j,xi[j]); */
     xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */
@@ -1578,7 +1598,14 @@ void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [
     p[j] += xi[j]; /* Parameters values are updated accordingly */
   } 
   /* printf("\n"); */
-  /* printf("Comparing last *frec(xmin)=%12.8f from Brent and frec(0.)=%12.8f \n", *fret, (*func)(p)); */
+#ifdef DEBUGLINMIN
+  printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
+  for (j=1;j<=n;j++) { 
+    printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]);
+    if(j % ncovmodel == 0)
+      printf("\n");
+  }
+#endif
   free_vector(xicom,1,n); 
   free_vector(pcom,1,n); 
 } 
@@ -1774,8 +1801,23 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
     } 
       if (directest < 0.0) { /* Then we use it for new direction */
+#endif
+#ifdef DEBUGLINMIN
+       printf("Before linmin in direction P%d-P0\n",n);
+       for (j=1;j<=n;j++) { 
+         printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+         if(j % ncovmodel == 0)
+           printf("\n");
+       }
 #endif
        linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+#ifdef DEBUGLINMIN
+       for (j=1;j<=n;j++) { 
+         printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+         if(j % ncovmodel == 0)
+           printf("\n");
+       }
+#endif
        for (j=1;j<=n;j++) { 
          xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
          xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
@@ -5860,7 +5902,7 @@ BOOL IsWow64()
 }
 #endif
 
-void syscompilerinfo()
+void syscompilerinfo(int logged)
  {
    /* #include "syscompilerinfo.h"*/
    /* command line Intel compiler 32bit windows, XP compatible:*/
@@ -5909,52 +5951,52 @@ void syscompilerinfo()
    int cross = CROSS;
    if (cross){
           printf("Cross-");
-          fprintf(ficlog, "Cross-");
+          if(logged) fprintf(ficlog, "Cross-");
    }
 #endif
 
 #include <stdint.h>
 
-   printf("Compiled with:");fprintf(ficlog,"Compiled with:");
+   printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:");
 #if defined(__clang__)
-   printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM");        /* Clang/LLVM. ---------------------------------------------- */
+   printf(" Clang/LLVM");if(logged)fprintf(ficlog," Clang/LLVM");      /* Clang/LLVM. ---------------------------------------------- */
 #endif
 #if defined(__ICC) || defined(__INTEL_COMPILER)
-   printf(" Intel ICC/ICPC");fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */
+   printf(" Intel ICC/ICPC");if(logged)fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */
 #endif
 #if defined(__GNUC__) || defined(__GNUG__)
-   printf(" GNU GCC/G++");fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */
+   printf(" GNU GCC/G++");if(logged)fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */
 #endif
 #if defined(__HP_cc) || defined(__HP_aCC)
-   printf(" Hewlett-Packard C/aC++");fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */
+   printf(" Hewlett-Packard C/aC++");if(logged)fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */
 #endif
 #if defined(__IBMC__) || defined(__IBMCPP__)
-   printf(" IBM XL C/C++"); fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */
+   printf(" IBM XL C/C++"); if(logged) fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */
 #endif
 #if defined(_MSC_VER)
-   printf(" Microsoft Visual Studio");fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */
+   printf(" Microsoft Visual Studio");if(logged)fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */
 #endif
 #if defined(__PGI)
-   printf(" Portland Group PGCC/PGCPP");fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */
+   printf(" Portland Group PGCC/PGCPP");if(logged) fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */
 #endif
 #if defined(__SUNPRO_C) || defined(__SUNPRO_CC)
-   printf(" Oracle Solaris Studio");fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */
+   printf(" Oracle Solaris Studio");if(logged)fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */
 #endif
-   printf(" for ");fprintf(ficlog," for ");
+   printf(" for "); if (logged) fprintf(ficlog, " for ");
    
 // http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros
 #ifdef _WIN32 // note the underscore: without it, it's not msdn official!
     // Windows (x64 and x86)
-   printf("Windows (x64 and x86) ");fprintf(ficlog,"Windows (x64 and x86) ");
+   printf("Windows (x64 and x86) ");if(logged) fprintf(ficlog,"Windows (x64 and x86) ");
 #elif __unix__ // all unices, not all compilers
     // Unix
-   printf("Unix ");fprintf(ficlog,"Unix ");
+   printf("Unix ");if(logged) fprintf(ficlog,"Unix ");
 #elif __linux__
     // linux
-   printf("linux ");fprintf(ficlog,"linux ");
+   printf("linux ");if(logged) fprintf(ficlog,"linux ");
 #elif __APPLE__
     // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though..
-   printf("Mac OS ");fprintf(ficlog,"Mac OS ");
+   printf("Mac OS ");if(logged) fprintf(ficlog,"Mac OS ");
 #endif
 
 /*  __MINGW32__          */
@@ -5968,11 +6010,11 @@ void syscompilerinfo()
 /* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */
 
 #if UINTPTR_MAX == 0xffffffff
-   printf(" 32-bit"); fprintf(ficlog," 32-bit");/* 32-bit */
+   printf(" 32-bit"); if(logged) fprintf(ficlog," 32-bit");/* 32-bit */
 #elif UINTPTR_MAX == 0xffffffffffffffff
-   printf(" 64-bit"); fprintf(ficlog," 64-bit");/* 64-bit */
+   printf(" 64-bit"); if(logged) fprintf(ficlog," 64-bit");/* 64-bit */
 #else
-   printf(" wtf-bit"); fprintf(ficlog," wtf-bit");/* wtf */
+   printf(" wtf-bit"); if(logged) fprintf(ficlog," wtf-bit");/* wtf */
 #endif
 
 #if defined(__GNUC__)
@@ -5985,18 +6027,18 @@ void syscompilerinfo()
                             + __GNUC_MINOR__ * 100)
 # endif
    printf(" using GNU C version %d.\n", __GNUC_VERSION__);
-   fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__);
+   if(logged) fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__);
 
    if (uname(&sysInfo) != -1) {
      printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
-     fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
+        if(logged) fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
    }
    else
       perror("uname() error");
    //#ifndef __INTEL_COMPILER 
 #if !defined (__INTEL_COMPILER) && !defined(__APPLE__)
    printf("GNU libc version: %s\n", gnu_get_libc_version()); 
-   fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version());
+   if(logged) fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version());
 #endif
 #endif
 
@@ -6004,12 +6046,12 @@ void syscompilerinfo()
    //   {
 #if defined(_MSC_VER)
    if (IsWow64()){
-          printf("The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
-          fprintf(ficlog, "The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
+          printf("\nThe program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
+          if (logged) fprintf(ficlog, "\nThe program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
    }
    else{
-          printf("The process is not running under WOW64 (i.e probably on a 64bit Windows).\n");
-          fprintf(ficlog,"The programm is not running under WOW64 (i.e probably on a 64bit Windows).\n");
+          printf("\nThe program is not running under WOW64 (i.e probably on a 64bit Windows).\n");
+          if (logged) fprintf(ficlog, "\nThe programm is not running under WOW64 (i.e probably on a 64bit Windows).\n");
    }
    //     printf("\nPress Enter to continue...");
    //     getchar();
@@ -6185,11 +6227,11 @@ int main(int argc, char *argv[])
   /*  FILE *fichtm; *//* Html File */
   /* FILE *ficgp;*/ /*Gnuplot File */
   struct stat info;
-  double agedeb;
+  double agedeb=0.;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
 
   double fret;
-  double dum; /* Dummy variable */
+  double dum=0.; /* Dummy variable */
   double ***p3mat;
   double ***mobaverage;
 
@@ -6199,18 +6241,18 @@ int main(int argc, char *argv[])
   char *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10;
   int c,  h , cpt;
-  int jl;
-  int i1, j1, jk, stepsize;
+  int jl=0;
+  int i1, j1, jk, stepsize=0;
   int *tab; 
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;
-  int hstepm, nhstepm;
+  int hstepm=0, nhstepm=0;
   int agemortsup;
   float  sumlpop=0.;
   double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
   double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;
 
-  double bage=0, fage=110, age, agelim, agebase;
+  double bage=0, fage=110., age, agelim=0., agebase=0.;
   double ftolpl=FTOL;
   double **prlim;
   double ***param; /* Matrix of parameters */
@@ -6271,7 +6313,7 @@ int main(int argc, char *argv[])
 #else
   getcwd(pathcd, size);
 #endif
-
+  syscompilerinfo(0);
   printf("\n%s\n%s",version,fullversion);
   if(argc <=1){
     printf("\nEnter the parameter file name: ");
@@ -6344,7 +6386,7 @@ int main(int argc, char *argv[])
  optionfilext=%s\n\
  optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
 
-  syscompilerinfo();
+  syscompilerinfo(0);
 
   printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);
@@ -6443,8 +6485,8 @@ int main(int argc, char *argv[])
   /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
   if(mle==-1){ /* Print a wizard for help writing covariance matrix */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
-    printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
-    fprintf(ficlog," You choose 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);
     free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
     fclose (ficparo);
     fclose (ficlog);