--- imach/src/imach.c 2014/08/27 16:26:55 1.157 +++ imach/src/imach.c 2014/09/15 20:41:41 1.161 @@ -1,6 +1,19 @@ -/* $Id: imach.c,v 1.157 2014/08/27 16:26:55 brouard Exp $ +/* $Id: imach.c,v 1.161 2014/09/15 20:41:41 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.161 2014/09/15 20:41:41 brouard + Summary: Problem with macro SQR on Intel compiler + + Revision 1.160 2014/09/02 09:24:05 brouard + *** empty log message *** + + Revision 1.159 2014/09/01 10:34:10 brouard + Summary: WIN32 + Author: Brouard + + Revision 1.158 2014/08/27 17:11:51 brouard + *** empty log message *** + Revision 1.157 2014/08/27 16:26:55 brouard Summary: Preparing windows Visual studio version Author: Brouard @@ -475,13 +488,18 @@ #include #include #include + +#ifdef _WIN32 +#include +#else #include +#endif #include #include #include #include -extern int errno; +/* extern int errno; */ /* #ifdef LINUX */ /* #include */ @@ -532,11 +550,11 @@ extern int errno; #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.157 2014/08/27 16:26:55 brouard Exp $ */ +/* $Id: imach.c,v 1.161 2014/09/15 20:41:41 brouard Exp $ */ /* $State: Exp $ */ char version[]="Imach version 0.98nX, August 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)"; -char fullversion[]="$Revision: 1.157 $ $Date: 2014/08/27 16:26:55 $"; +char fullversion[]="$Revision: 1.161 $ $Date: 2014/09/15 20:41:41 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -774,7 +792,7 @@ char *cutl(char *blocc, char *alocc, cha gives blocc="abcdef2ghi" and alocc="j". If occ is not found blocc is null and alocc is equal to in. Returns blocc */ - char *s, *t, *bl; + char *s, *t; t=in;s=in; while ((*in != occ) && (*in != '\0')){ *alocc++ = *in++; @@ -858,6 +876,24 @@ int nbocc(char *s, char occ) /* } */ /* } */ +#ifdef _WIN32 +char * strsep(char **pp, const char *delim) +{ + char *p, *q; + + if ((p = *pp) == NULL) + return 0; + if ((q = strpbrk (p, delim)) != NULL) + { + *pp = q + 1; + *q = '\0'; + } + else + *pp = 0; + return p; +} +#endif + /********************** nrerror ********************/ void nrerror(char error_text[]) @@ -1080,7 +1116,7 @@ double brent(double ax, double bx, doubl { int iter; double a,b,d,etemp; - double fu,fv,fw,fx; + double fu=0,fv,fw,fx; double ftemp; double p,q,r,tol1,tol2,u,v,w,x,xm; double e=0.0; @@ -1367,23 +1403,43 @@ void powell(double p[], double **xi, int return; } if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); - for (j=1;j<=n;j++) { + for (j=1;j<=n;j++) { /* Computes an extrapolated point */ ptt[j]=2.0*p[j]-pt[j]; xit[j]=p[j]-pt[j]; pt[j]=p[j]; } fptt=(*func)(ptt); - if (fptt < fp) { - t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); - if (t < 0.0) { - linmin(p,xit,n,fret,func); + if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */ + /* x1 f1=fp x2 f2=*fret x3 f3=fptt, xm fm */ + /* From x1 (P0) distance of x2 is at h and x3 is 2h */ + /* Let f"(x2) be the 2nd derivative equal everywhere. Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) + will reach at f3 = fm + h^2/2 f''m ; f" = (f1 -2f2 +f3 ) / h**2 */ + /* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */ + /* Thus we compare delta(2h) with observed f1-f3 */ + /* or best gain on one ancient line 'del' with total gain f1-f2 = f1 - f2 - 'del' with del */ + /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */ + t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); + t= t- del*SQR(fp-fptt); + printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t); + fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t); +#ifdef DEBUG + printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), + (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt)); + fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), + (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt)); + printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); + fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); +#endif + if (t < 0.0) { /* Then we use it for last direction */ + linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/ for (j=1;j<=n;j++) { - xi[j][ibig]=xi[j][n]; - xi[j][n]=xit[j]; + xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */ + xi[j][n]=xit[j]; /* and nth direction by the extrapolated */ } + printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction :\n",n,ibig); + #ifdef DEBUG - printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); - fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); for(j=1;j<=n;j++){ printf(" %.12e",xit[j]); fprintf(ficlog," %.12e",xit[j]); @@ -6536,7 +6592,8 @@ Interval (in months) between two waves: for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ oldm=oldms;savm=savms; /* Segmentation fault */ - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); + cptcod= 0; /* To be deleted */ + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); if(vpopbased==1) fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); @@ -6696,15 +6753,15 @@ Interval (in months) between two waves: sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); #endif if(!stat(plotcmd,&info)){ - printf("Error or gnuplot program not found: %s\n",plotcmd);fflush(stdout); + printf("Error or gnuplot program not found: '%s'\n",plotcmd);fflush(stdout); if(!stat(getenv("GNUPLOTBIN"),&info)){ - printf("Error or gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout); + printf("Error or gnuplot program not found: '%s' Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout); }else strcpy(pplotcmd,plotcmd); #ifdef __unix strcpy(plotcmd,GNUPLOTPROGRAM); if(!stat(plotcmd,&info)){ - printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); + printf("Error gnuplot program not found: '%s'\n",plotcmd);fflush(stdout); }else strcpy(pplotcmd,plotcmd); #endif @@ -6712,23 +6769,23 @@ Interval (in months) between two waves: strcpy(pplotcmd,plotcmd); sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); - printf("Starting graphs with: %s\n",plotcmd);fflush(stdout); + printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); if((outcmd=system(plotcmd)) != 0){ - printf("gnuplot command might not be in your path: %s, err=%d\n", plotcmd, outcmd); + printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); if((outcmd=system(plotcmd)) != 0) printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); } - printf(" Successul, please wait..."); + printf(" Successful, please wait..."); while (z[0] != 'q') { /* chdir(path); */ printf("\nType e to edit results with your browser, g to graph again and q for exit: "); scanf("%s",z); /* if (z[0] == 'c') system("./imach"); */ if (z[0] == 'e') { -#ifdef _APPLE_ +#ifdef __APPLE__ sprintf(pplotcmd, "open %s", optionfilehtm); #elif __linux sprintf(pplotcmd, "xdg-open %s", optionfilehtm); @@ -6747,6 +6804,3 @@ Interval (in months) between two waves: scanf("%s",z); } } - - -