--- imach/src/imach.c 2015/05/05 08:51:13 1.190 +++ imach/src/imach.c 2015/07/14 10:00:33 1.191 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.190 2015/05/05 08:51:13 brouard Exp $ +/* $Id: imach.c,v 1.191 2015/07/14 10:00:33 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.191 2015/07/14 10:00:33 brouard + Summary: Some fixes + Revision 1.190 2015/05/05 08:51:13 brouard Summary: Adding digits in output parameters (7 digits instead of 6) @@ -686,11 +689,11 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.190 2015/05/05 08:51:13 brouard Exp $ */ +/* $Id: imach.c,v 1.191 2015/07/14 10:00:33 brouard Exp $ */ /* $State: Exp $ */ char version[]="Imach version 0.98q2, April 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.190 $ $Date: 2015/05/05 08:51:13 $"; +char fullversion[]="$Revision: 1.191 $ $Date: 2015/07/14 10:00:33 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1450,31 +1453,41 @@ values at the three points, fa, fb , and #endif #ifdef MNBRAKORIGINAL #else - 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 */ +/* 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("mnbrak3 fu < fc \n"); - fprintf(ficlog, "mnbrak3 fu < fc\n"); + printf("mnbrak34 fu < or >= fc \n"); + fprintf(ficlog, "mnbrak34 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 @@ -1565,6 +1578,9 @@ void linmin(double p[], double xi[], int } }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]) */ @@ -1574,7 +1590,9 @@ void linmin(double p[], double xi[], int 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) */ @@ -1583,7 +1601,14 @@ void linmin(double p[], double xi[], int 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); } @@ -1780,7 +1805,22 @@ void powell(double p[], double **xi, int } 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 */ @@ -5865,7 +5905,7 @@ BOOL IsWow64() } #endif -void syscompilerinfo() +void syscompilerinfo(int logged) { /* #include "syscompilerinfo.h"*/ /* command line Intel compiler 32bit windows, XP compatible:*/ @@ -5914,52 +5954,52 @@ void syscompilerinfo() int cross = CROSS; if (cross){ printf("Cross-"); - fprintf(ficlog, "Cross-"); + if(logged) fprintf(ficlog, "Cross-"); } #endif #include - 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__ */ @@ -5973,11 +6013,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__) @@ -5990,18 +6030,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 @@ -6009,12 +6049,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(); @@ -6190,11 +6230,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; @@ -6204,18 +6244,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 */ @@ -6276,7 +6316,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: "); @@ -6349,7 +6389,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); @@ -6448,8 +6488,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);