#include #include #include #include #include #include /* for command line option stuff */ #include "fly_opt.h" #include "global.h" #include #include /* DistP.variables and prototypes */ /*** Constants *************************************************************/ /* command line option string */ /* #define OPTS ":a:b:Bc:C:d:De:E:f:g:hi:lL:nopQr:s:Stvw:W:y:" */ #define OPTS ":A:a:b:Bc:C:De:E:f:g:G:hi:I:lL:nNopQ:R:r:s:tvw:V:W:x:X:y:z:" /* command line option string */ /* D will be debug, like scramble, score */ /* must start with :, option with argument must have a : following */ /*** STATIC VARIABLES ******************************************************/ /* Help, usage and version messages */ #ifdef MPI static const char usage[] = "Usage: FLY_**.mpi [-b ] [-B] [-C ] \n" " [-D] [-e ][-E ] [-f ] [-g ]\n" " [-h] [-i ] [-i ][-l] [-L ] [-n] [-N] [-p] \n" " [-r rmsflag] [-s ] [-S] [-t] [-T ] [-v] [-w ]\n" " [-W ] [-y ]\n" " [-j] [-r ]\n" " [-c ]\n" " \n"; #else static const char usage[] = "Usage: FLY_** [-a ] [-b ] [-B] [-e ] [-E ]\n" " [-f ] [-g ] [-h] [-i ][-I ]\n" " [-l] [-L ] [-n] [-N] [-p] [-Q] [-r rmsflag] [-s ] [-v]\n" " [-w ] [-y ]\n" " [-j] [-r ]\n" " [-c ]\n" " \n"; #endif static const char help[] = #ifdef MPI "Usage: fly_***.mpi [options] \n\n" #else "Usage: fly_*** [options] \n\n" #endif "Argument:\n" " input data file\n\n" "Options:\n" " -a solver accuracy for adaptive stepsize ODE solvers\n" " -b write state file every * tau moves\n" " -B run in benchmark mode (only do fixed initial steps)\n" #ifdef MPI " -C set covar sample interval to * tau\n" #endif " -D debugging mode, prints all kinds of debugging info\n" " -e set annealing freeze criterion to \n" " -E define ES_TYPE \n" " -f float precision of parameters is \n" " -g chooses g(u): e = exp, h = hvs, s = sqrt, t = tanh\n" " -h prints this help message\n" " -i sets ODE solver stepsize (in minutes)\n" " -I <0/1> set individual Or integrated data type to use\n" " -l echo log to the terminal\n" #ifdef MPI " -L define LM_TYPE \n" #endif " -n nofile: don't print .log or .state files\n" " -N generates landscape to .landscape file in equilibrate mode \n" " -o use oldstyle cell division times (3 div only)\n" " -p prints move acceptance stats to .prolix file\n" #ifndef MPI " -Q quenchit mode, T is lowered immediately to zero\n" #endif " -s choose ODE solver\n" #ifdef MPI " -S disable tuning stop flag\n" #endif " -t write timing information to .times file\n" " -v print version and compilation date\n" " -w write output to instead of \n" #ifdef MPI " -W tuning stats written times per interval\n" #endif " -y write log every * tau moves\n" " -j Allow gnuplot visualization of the simulation\n" " -r write residual computation in file\n" " -c by default it is the LSE: \n" " - lse or LSE \n " " - mad or MAD\n\n " "Please report bugs to . Thank you!\n"; static char version[MAX_RECORD]; /* version gets set below */ static const char verstring[] = "compiled by: %s\n" " on: %s\n" " using: %s\n" " flags: %s\n" " date: %s at %s\n"; /* other static variables */ static char *inname; /* filename of input file */ static char *outDir,*outname; /* filename of output file */ static char *ga_config; /*filename where ga is configured*/ static int ga_type; static int ls_type; static char *argvsave; /* static string for saving command line */ static double energy; static double stepsize = 1.; /* stepsize for solver */ static double accuracy = 0.001; /* accuracy for solver (not used yet) */ static int precision = 8; /* precision for eqparms */ static int prolix_flag = 0; /* to prolix or not to prolix */ static int landscape_flag = 1; /* generate energy landscape data */ /* set the landscape flag (and the landscape filename) in lsa.c */ /* the following lines define a pointers to: */ /* - pd: dvdt function, currently only DvdtOrig in zygotic.c */ /* - pj: Jacobian function, in zygotic.c */ /* This pointer needs to be static since both InitialMove and Restore- */ /* State need it; it gets set in ParseCommandLine() */ /* */ /* NOTE: ps (solver) is declared as global in integrate.h */ //static void (*pd)(double *, double, double *, int); //static void (*pj)(double, double *, double *, double **, int); /* log and input/output file related variables *****************************/ static char *inputfile; /* name of the input file */ static char *outputfile; /* name of the output file */ static char *statefile; /* name of the state file */ static char *logfile; /* name of the global .log file */ static char *landscapefile; /* filename of landscape file */ #ifdef MPI static char *l_logfile; /* name of the local .llog file */ /* Files for tuning */ static char *lbfile; /* name of .lb file (for lower_bound) */ static char *ubfile; /* name of .ub file (for upper_bound) */ static char *mbfile; /* name of .mb file (for mix_bound) */ #endif /* flag used by Landscape generation ****************************************/ /* Set by InitLandscape called from xxx_sa.c****************************/ static int landscape = 0; static int prolix; /* flag for printing stat info (prolix) */ static char *prolixfile; /* filename of prolix file */ StopStyle stop_flag; /* type of stop criterion (see above) */ static int ndp; static double penalty; int rmsflag = 0 ; /* flag for displaying log to stdout */ int time_flag; /* flag for timing the code */ int log_flag; /* flag for displaying log to stdout */ int nofile_flag; /* flag for not writing .state or .log files */ long state_write; /* frequency for writing state files (in tau) */ long print_freq; /* frequency for printing to log (in tau) */ long captions; /* option for printing freqency of log captions */ /* Special annealing modes: ************************************************ * * * the following are flags that should be set by cmd line opts for some * * special annealing modes that we use for testing or gathering stats: * * * * - bench: does no loop, i.e. only runs through initial steps; this * * can be useful to do solver benchmarks in -B mode * * - equil: runs at a fixed temperature to collect equilibrium stats * * - quenchit: means lower the temperature to zero *instantly* after fini- * * shing initialize; this is not implemented for parallel code * * since it doesn't take a long time to finish * * * ***************************************************************************/ int bench; int equil; int quenchit; int IndivOrIntegrated; int restart=0; static ParamList *ptab; /* array of pointers to parameters and ranges */ char *section_title; /* parameter section name */ void WriteVersion(char *filename) { char *temp; /* temporary file name */ char *record; /* record to be read and written */ char *record_ptr; /* pointer used to remember record for 'free' */ char *convline; /* temporarily saves conv version line */ char *shell_cmd; /* used by 'system' below */ FILE *outfile; /* name of output file */ FILE *tmpfile; /* name of temporary file */ int versionflag = 0; /* version section already present or not? */ temp = (char *)calloc(MAX_RECORD, sizeof(char)); record = (char *)calloc(MAX_RECORD, sizeof(char)); convline = (char *)calloc(MAX_RECORD, sizeof(char)); shell_cmd = (char *)calloc(MAX_RECORD, sizeof(char)); record_ptr = record; /* this is to remember record for 'free' */ outfile = fopen(filename, "r"); /* open outfile for reading */ if ( !outfile ) /* sorry for the confusion! */ error("WriteVersion: error opening %s", filename); if (FindSection(outfile, "version")) /* version section already there? */ versionflag = 1; rewind(outfile); temp = strcpy(temp,"versionXXXXXX"); /* required by mkstemp() */ if ( mkstemp(temp) == -1 ) /* get unique name for temp file */ error("WriteVersion: error creating temporary file"); tmpfile = fopen(temp, "w"); /* ... and open it for writing */ if ( !tmpfile ) error("WriteVersion: error opening temporary file %s", temp); if ( !versionflag ) { /* CASE 1: no version section yet */ fprintf(tmpfile, "$version\n"); /* -> write new one at beginning */ fprintf(tmpfile, "%s\n", version); /* of the data file */ fprintf(tmpfile, "%s", argvsave); fprintf(tmpfile, "$$\n\n"); while ( (record = fgets(record, MAX_RECORD, outfile)) ) fputs(record, tmpfile); /* CASE 2: version section already present */ } else { /* -> append new lines to existing section */ while ( strncmp(record=fgets(record, MAX_RECORD, outfile), "$version", 8) ) /* first write everything before */ fputs(record, tmpfile); /* version section */ fputs(record, tmpfile); while ( strncmp(record=fgets(record, MAX_RECORD, outfile), "$$", 2) ) /* save existing converter line */ if ( !strncmp(record, "converted", 9) ) /* skip all the rest */ convline = strcpy(convline, record); fprintf(tmpfile, "%s\n", version); /* new lines are written here */ fprintf(tmpfile, "%s", argvsave); if ( strlen(convline) > 0 ) /* conversion line is appended at end */ fprintf(tmpfile, "%s", convline); fprintf(tmpfile, "$$\n"); while ( (record=fgets(record, MAX_RECORD, outfile)) ) fputs(record, tmpfile); /* ... and then write all the rest */ } fclose(outfile); fclose(tmpfile); /* rename tmpfile into new file */ sprintf(shell_cmd, "cp -f %s %s", temp, filename); if ( -1 == system(shell_cmd) ) error("WriteVersion: error renaming temp file %s", temp); if ( remove(temp) ) warning("WriteVersion: temp file %s could not be deleted", temp); /* clean up */ free(temp); free(record_ptr); free(convline); free(shell_cmd); } /*** WriteTimes: writes the time-structure to a .times file **************** ***************************************************************************/ void WriteTimes(double *times) { char *timefile; /* name of the .times file */ FILE *timeptr; /* file pointer for .times file */ /* create time file name by appending .times to input file name */ timefile = (char *)calloc(MAX_RECORD, sizeof(char)); timefile = strcpy(timefile, outname); timefile = strcat(timefile, ".times"); timeptr = fopen(timefile, "w"); if ( !timeptr ) file_error("main"); PrintTimes(timeptr, times); /* write times to .times file */ fclose(timeptr); /* clean up */ free(timefile); } /*** PrintTimes: writes two (parallel: three) times sections *************** ***************************************************************************/ void PrintTimes(FILE *fp, double *times) { fprintf(fp, "wallclock: %.3f\n", times[0]); fprintf(fp, "user: %.3f\n", times[1]); } /*** FUNCTIONS USED TO COMMUNICATE WITH OTHER FILES ***********************/ /**************************************************************************** *** InitLandscape: sets flag for printing landscape output and acceptance**** * landscape and initializes the landscape file names *** * called from xxx_sa.c to make filenames static and *** * set landscape flag *** ****************************************************************************/ void InitLandscape(int value, char *file) { const char *suffix = ".landscape"; /* landscape file in equilibrate */ /* sets the landscape and acceptance landscape file names static to lsa.c */ landscapefile = (char *)calloc(MAX_RECORD, sizeof(char)); landscapefile = strcpy(landscapefile, file); landscapefile = strcat(landscapefile, suffix); landscape = value; } /*** FUNCTIONS THAT COMMUNICATE WITH OTHER SOURCE FILES ********************/ /*** SetProlix: sets flag for printing prolix output on acceptance stats *** * and initializes the prolix file if required * ***************************************************************************/ void SetProlix(int value, char *file, int init_flag) { FILE *prolixptr; const char *suffix = ".prolix"; /* suffix for prolix file */ prolixfile = (char *)calloc(MAX_RECORD, sizeof(char)); prolixfile = strcpy(prolixfile, file); prolixfile = strcat(prolixfile, suffix); if ( init_flag ) { /* this deletes an old .prolix file if present */ prolixptr = fopen(prolixfile, "w"); if ( !prolixptr ) file_error("SetProlix"); fclose(prolixptr); } prolix = value; } void InitDistribution(FILE *fp) { /* begin init distribution */ /* read in the data from file */ /* read in header info then output it*/ /* At the end of file exit from the loop */ printf("init distribution_parameters\n"); fp = FindSection(fp,"distribution_parameters"); /* find distribution section */ if( !fp ) /* to be backward compatible set default to exp */ { printf("ReadTune: could not locate distribution parameters:using exponential.\n"); DistP.distribution = 1; DistP.q == 1.0; } /* end of default to exp */ else { /* input user selected distribution parameters */ fscanf(fp,"%*s\n"); /* advance past title line no blanks allowed! */ /* read distribution stuff */ if ( 2!= (fscanf(fp,"%d %lf\n", &(DistP.distribution), &(DistP.q) )) ) error ("ReadTune: error reading distribution stuff"); if ((DistP.distribution > 11) || (DistP.distribution < 1)) {error ("fly_***: distribution must be int from 1 to 11 \n"); } else if ((DistP.distribution == 4)||(DistP.distribution == 3)) {error ("fly_***: PLEASE use 5 for Lorentz or 10 for normal distribution \n"); } else if ((DistP.distribution == 6)||(DistP.distribution == 9)) {error ("fly_***: 6=poisson or 9=pareto distribution returns positive values--do not use for fly \n"); } else if (DistP.distribution == 7 ) { /* general distribution */ if ((DistP.q >= 3.0) || (DistP.q <= 1.0)) {error ("tsp_sa: q must be between 1 and 3 \n");} else if (DistP.q == 2.0) {DistP.distribution = 5; /* fly needs lorentz, tsp use abs lorentz(4)*/ printf ("fly_***: q=2 is lorentz--setting distribution to 5\n");} else if (DistP.q > 2.0) {qgt2_init(); } else {qlt2_init(); } } /* end of general distribution */ /***************LG 05-02***************/ /* calculate q dependent factors that */ /* do not change for entire run. */ /* these live in distributions.h */ /***************LG 05-02***************/ } /* end of input user selected distribution parameters */ } /* end InitDistribution */ /*** Initial: initializes the following stuff: ************************* * - various static things for filenames and such * ***************************************************************************/ double InitialMove(int argc, char **argv, int opt_index, NucStatePtr state_ptr, double *p_chisq) { FILE *infile; int i; char *p; char * shell_cmd; char *basename; int mkDirOut; double i_temp; double energy; GAType in_tune; /* temporary struct to read in tune parameters */ LocalParam in_lm; PArrPtr pl; /* local copy of PArrPtr */ /* save some things in static storage for FinalMove() and StateWrite() */ inname = (char *)calloc(MAX_RECORD + 1, sizeof(char)); /* input file */ inname = strcpy(inname, argv[opt_index]); if ( !outDir) { /* in case there's no -w, outname = inname */ outDir = (char *)calloc(MAX_RECORD + 1, sizeof(char)); outDir = strcpy(outDir, inname); outDir=strcat(outDir,"_out"); }else{ //outDir=strcat(outDir,"_esds"); } basename = (char *)calloc(MAX_RECORD, sizeof(char)); outname = (char *)calloc(MAX_RECORD, sizeof(char)); shell_cmd = (char *)calloc(MAX_RECORD, sizeof(char)); // sprintf(shell_cmd, "mkdir %s ",outDir); mkDirOut= system(shell_cmd); printf ("The value returned was: %d.\n",mkDirOut); switch (mkDirOut){ case -1: error("fly_ea: error creating the directory file %s", outDir); break; case 256: sprintf(shell_cmd, "rm %s/* ",outDir); system(shell_cmd); sprintf(shell_cmd, "rm %s/*.*",outDir); system(shell_cmd); break; case 0: break; } basename=strrchr(outDir, '/') +1; outname=strncat(outDir,"/",strlen(basename)-1); outname=strncat(outname,basename,strlen(basename)-1); // outname=strcpy(outname,outname); sprintf(shell_cmd, "cp -f %s %s", inname, outname); if ( -1 == system(shell_cmd) ) error("fly_ea: error creating the output file %s",outname); free(shell_cmd); argvsave = (char *)calloc(MAX_RECORD, sizeof(char)); for(i=0; i < argc; i++){ if ( i>0 ) argvsave = strcat(argvsave, " "); argvsave = strcat(argvsave, argv[i]); } argvsave = strcat(argvsave, "\n"); /* set the prolix flag (and the .prolix filename) in moves.c, if necessary */ if ( prolix_flag ) SetProlix(prolix_flag, outname, 1); /* 1 means: new .prolix file */ if ( landscape_flag ) /* 1 means gen landscape files */ InitLandscape(landscape_flag, outname); /*lives in lsa.c */ /* initialize some Lam/Greening structures */ p = state_ptr->tune.progname; /* tune.progname contains program name */ p = strcpy(p, version); state_ptr->tune.debuglevel = 0; /* following stuff not used now */ p = state_ptr->tune.tunename; p = strcpy(p, "The Other One"); /* Grateful Dead tune, what else? */ state_ptr->tune.tunefile = NULL; /* read data file and initialize different things for scoring and move */ /* generation (in zygotic.c, score.c and moves.c */ printf("attempt to open: %s\n",argv[opt_index]); infile = fopen(argv[opt_index], "r"); if ( !infile ) file_error("fly_*** but"); InitZygote(infile, pd, pj, section_title); // InitZygote2(infile, "input"); /* init problem, parms, bias, bcd, static structs for deriv and deriv */ InitScoring(infile); /* initializes facts and limits */ InitStepsize(stepsize, accuracy, NULL, argv[opt_index]); InitTweak(infile); in_tune = ReadTune(infile); /* read tune_parameter section */ in_lm=ReadLocalParam(infile); pl = Translate(); /* read the list of parameters to be tweaked */ // nparams = pl.size; /* nparams is static to moves.c */ ptab = pl.array; /* make parameter array static to moves.c */ in_tune.len_chrom=pl.size; // i_temp = InitMoves(infile); /* set initial temperature and initialize */ /* random number generator and annealing parameters */ InitDistribution(infile); /* initialize distribution stuff */ fclose(infile); /* initialize Lam parameters (see sa.h for further detail) */ state_ptr->tune.restart=restart; state_ptr->tune.inname=inname; state_ptr->tune.outname=outname; state_ptr->tune.outNamePop = (char *)calloc(MAX_RECORD, sizeof(char)); state_ptr->tune.outNamePop=strcpy( state_ptr->tune.outNamePop,outname); state_ptr->tune.outNamePop=strcat( state_ptr->tune.outNamePop,"_pop"); state_ptr->tune.outNameLandscapeFirst = (char *)calloc(MAX_RECORD, sizeof(char)); state_ptr->tune.outNameLandscapeFirst=strcpy( state_ptr->tune.outNameLandscapeFirst,outname); state_ptr->tune.outNameLandscapeFirst=strcat( state_ptr->tune.outNameLandscapeFirst,"_landscapeFirst"); state_ptr->tune.outNameLandscapeSec = (char *)calloc(MAX_RECORD, sizeof(char)); state_ptr->tune.outNameLandscapeSec=strcpy( state_ptr->tune.outNameLandscapeSec,outname); state_ptr->tune.outNameLandscapeSec=strcat( state_ptr->tune.outNameLandscapeSec,"_landscapeSec"); state_ptr->tune.configOut = (char *)calloc(MAX_RECORD, sizeof(char)); state_ptr->tune.configOut=strcpy( state_ptr->tune.configOut,outname); state_ptr->tune.configOut=strcat( state_ptr->tune.configOut,"_config"); state_ptr->tune.outDir=outDir; state_ptr->tune.evoltype= in_tune.evoltype; state_ptr->tune.seed=in_tune.seed; state_ptr->tune.nbPopulation= in_tune.nbPopulation; state_ptr->tune.populationSize= in_tune.populationSize; state_ptr->tune.num_chrom= in_tune.num_chrom; state_ptr->tune.len_chrom= in_tune.len_chrom; state_ptr->tune.crossover= in_tune.crossover; state_ptr->tune.mutation= in_tune.mutation; state_ptr->tune.migration= in_tune.migration; state_ptr->tune.replace= in_tune.replace; state_ptr->tune.nbGeneration=in_tune.nbGeneration; state_ptr->tune.minGeneration=in_tune.minGeneration; state_ptr->tune.steadyGen=in_tune.steadyGen; state_ptr->tune.maxFuncEval =in_tune.maxFuncEval; state_ptr->tune.convergence= in_tune.convergence; state_ptr->tune.popConvergence= in_tune.popConvergence; state_ptr->tune.monOut=in_tune.monOut; state_ptr->tune.fileOut=in_tune.fileOut; state_ptr->tune.plotOut=in_tune.plotOut; state_ptr->tune.lambda_rate= in_tune.lambda_rate; state_ptr->tune.mu= in_tune.mu; switch (ga_type){ case 0: state_ptr->tune.evoltype=MuLambda; break; case 1: state_ptr->tune.evoltype=MuPlusLambda; break; default: error("unknow type of ES optimizatiom method"); break; } switch (ls_type){ case 0: state_ptr->tune.lstype=SHHSCH; break; case 1: state_ptr->tune.lstype=SMDSCH; break; case 2: state_ptr->tune.lstype=NMSCH; break; case 3: state_ptr->tune.lstype= EDHJ; break; case 4: state_ptr->tune.lstype= HJ; break; case 5: state_ptr->tune.lstype=NLESS; break; case 6: state_ptr->tune.lstype=COMPSS; break; case 7: state_ptr->tune.lstype=COORD; break; case 8: state_ptr->tune.lstype= HNM; break; default: error("unknow type of local search optimizatiom method"); break; } state_ptr->tune.mix_interval= in_tune.mix_interval; state_ptr->tune.energy_threshold = in_tune.energy_threshold; state_ptr->tune.count_threshold = in_tune.count_threshold; state_ptr->tune.check_interval = in_tune.check_interval; state_ptr->lmParam.nbIterations=in_lm.nbIterations; state_ptr->lmParam.mu=in_lm.mu; state_ptr->lmParam.epsilon1=in_lm.epsilon1; state_ptr->lmParam.epsilon2=in_lm.epsilon2; state_ptr->lmParam.epsilon3=in_lm.epsilon3; state_ptr->lmParam.delta=in_lm.delta; energy = Score(); /* score first time to check if: */ /* - parameters within range? */ if ( energy == FORBIDDEN_MOVE ) /* - set initial energy (p_chisq) */ error("fly_***: the initial state was forbidden, cannot proceed"); ndp = GetNDatapoints(); penalty = GetPenalty(); *p_chisq = energy; /* set initial score */ return(i_temp); /* initial temperature */ } /*** ReadTune: reads the tune_parameters section in a data file and ******** * turns a GAType structure to the caller * ***************************************************************************/ GAType ReadTune(FILE *fp) { GAType in_tune= InitDefaultGA(); /* this is the GAType struct that we return */ unsigned ui; long lbuf; int intbuf1; /* following four are used as temp. buffer for ints */ int intbuf2; int intbuf3; int intbuf4; double dbuf1; /* following four are used as temp. buffer for doubles */ double dbuf2; double dbuf3; double dbuf4; fp = FindSection(fp, "ga_tune_parameters"); /* find tune section */ if( !fp ) error("ReadTune: could not locate ga_tune_parameters"); fscanf(fp,"%*s\n"); /* advance past title line 1 */ /* read seed */ if ( 1 != fscanf(fp, "%u\n",&ui)) error("ReadTune: error reading seed interval"); in_tune.seed=ui; fscanf(fp,"%*s\n"); /* advance past title line 2 */ /* read number and size of populations, and number of chromosomes */ if ( 3 !=(fscanf(fp,"%d %d %d\n",&intbuf1,&intbuf2,&intbuf3))) error ("ReadTune: error reading GA stuff"); in_tune.nbPopulation=intbuf1; in_tune.populationSize=intbuf2; in_tune.num_chrom=intbuf3; fscanf(fp,"%*s\n"); /* advance past title line 3 */ /* read recombination and mutation rate, replacement and migration proportions */ if ( 4!= (fscanf(fp,"%lg %lg %lg %lg\n", &dbuf1, &dbuf2, &dbuf3, &dbuf4)) ) error ("ReadTune: error reading GA stuff"); in_tune.crossover = dbuf1; in_tune.mutation = dbuf2; in_tune.replace = dbuf3; in_tune.migration = dbuf4; fscanf(fp,"%*s\n"); /* advance past title line 4*/ /* read nbGeneration_minGeneration_steadyGen_maxFuncEval */ if ( 4!= (fscanf(fp,"%d %d %d %d\n", &intbuf1, &intbuf2, &intbuf3, &intbuf4)) ) error ("ReadTune: error reading GA stuff"); in_tune.nbGeneration = intbuf1; in_tune.minGeneration=intbuf2; in_tune.steadyGen=intbuf3; in_tune.maxFuncEval=intbuf4; fscanf(fp,"%*s\n"); /* advance past title line 5 */ /* read convergence criteria */ if ( 2!= (fscanf(fp,"%lg %lg \n", &dbuf1, &dbuf2)) ) error ("ReadTune: error reading GA stuff"); in_tune.convergence = dbuf1; in_tune.popConvergence = dbuf2; fscanf(fp,"%*s\n"); /* advance past title line 6 */ /* read monOut,fileOut and plotOut */ if ( 3!= (fscanf(fp,"%d %d %d\n", &intbuf1, &intbuf2, &intbuf3)) ) error ("ReadTune: error reading GA stuff"); in_tune.monOut = intbuf1; in_tune.fileOut=intbuf2; in_tune.plotOut=intbuf3; fscanf(fp,"%*s\n"); /* advance past title line 7 */ /* read lambda and mu */ if ( 2!= (fscanf(fp,"%lg %d \n", &dbuf1, &intbuf2)) ) error ("ReadTune: error reading GA stuff"); in_tune.lambda_rate = dbuf1; in_tune.mu=intbuf2; fscanf(fp,"%*s\n"); /* advance past title line 8 */ /* read migration interval */ if ( 1 != fscanf(fp, "%d\n", &intbuf1) ) error("ReadTune: error reading mixing interval"); in_tune.mix_interval = intbuf1; fscanf(fp,"%*s\n"); /* advance past title line 9 */ /* read stop criteria stuff */ if ( 3!= (fscanf(fp,"%lg %d %d \n", &dbuf1, &intbuf2,&intbuf3)) ) error ("ReadTune: error reading GA stuff"); in_tune.energy_threshold = dbuf1; in_tune.count_threshold = intbuf2; in_tune.check_interval=intbuf3; in_tune.tunefile = NULL; /* We're not using this stuff */ in_tune.debuglevel = 0; return (in_tune); } LocalParam ReadLocalParam(FILE *fp){ LocalParam in_lm= InitDefaultLM(); /* this is the GAType struct that we return */ int intbuf1; /* following four are used as temp. buffer for ints */ double dbuf1; /* following four are used as temp. buffer for doubles */ double dbuf2; double dbuf3; double dbuf4; double dbuf5; fp = FindSection(fp, "local_search_parameters"); /* find tune section */ if( !fp ) error("ReadTune: could not locatelocal_search_parameters"); fscanf(fp,"%*s\n"); /* advance past title line 1 */ if ( 1 !=(fscanf(fp,"%d\n",&intbuf1))) error ("ReadTune: error reading LM stuff"); in_lm.nbIterations=intbuf1; fscanf(fp,"%*s\n"); /* advance past title line 1 */ /* read Lam lambda stuff */ if ( 5!= (fscanf(fp,"%lg %lg %lg %lg %lg\n", &dbuf1, &dbuf2, &dbuf3, &dbuf4, &dbuf5)) ) error ("ReadTune: error reading LM stuff"); in_lm.mu = dbuf1; in_lm.epsilon1 = dbuf2; in_lm.epsilon2 = dbuf3; in_lm.epsilon3 = dbuf4; in_lm.delta=dbuf5; return (in_lm); } /* FUNCTIONS *************************************************************/ /*** COMMAND LINE OPTS ARE ALL PARSED HERE *********************************/ /*** ParseCommandLine: well, parses the command line and returns an index ** * to the 1st argument after the command line options * ***************************************************************************/ GAType InitDefaultGA(){ GAType in_tune; /* in_tune.scheme=GA_SCHEME_DARWIN; in_tune.elit =GA_ELITISM_PARENTS_DIE; in_tune.strategy=GA_DE_STRATEGY_RANDTOBEST; */ in_tune.evoltype=0 ; in_tune.lstype=3; in_tune.nbPopulation=1; in_tune.populationSize=10; in_tune.num_chrom=1; in_tune.len_chrom=4; in_tune.crossover=0.9; in_tune.mutation=0.01; in_tune.migration=0.0; in_tune.replace=0.0; in_tune.nbGeneration=100; in_tune.convergence=0; in_tune.popConvergence=0; return in_tune; } LocalParam InitDefaultLM(){ LocalParam in_lm; in_lm.nbIterations=100; in_lm.mu=1E-03; in_lm.epsilon1=1E-17 ; in_lm.epsilon2=1E-17 ; in_lm.epsilon3=1E-17 ; in_lm.delta=0; return in_lm; } int ParseCommandLine(int argc, char **argv) { int c; /* used to parse command line options */ /* external declarations for command line option parsing (unistd.h) */ extern char *optarg; /* command line option argument */ extern int optind; /* pointer to current element of argv */ extern int optopt; /* contain option character upon error */ /* set the version string */ /* following part sets default values for command line options */ IndivOrIntegrated =INTEGRATED; ga_type=0; pd = DvdtOrig; /* default derivative function */ pj = JacobnOrig; /* default Jacobian function */ ps = Rk4; /* default solver */ captions = 100000000; /* default freq for writing captions (off) */ print_freq = 100; /* default freq for writing to log file */ state_write = 1000; /* default freq for writing to state file */ stop_flag = absolute_freeze; /* type of stop criterion */ time_flag = 0; /* flag for timing */ log_flag = 0; /* flag for writing logs to the screen */ nofile_flag = 0; /* flog for not writing .log and .state files */ #ifdef MPI tuning = 0; /* tuning mode is off by default */ covar_index = 1; /* covariance sample will be covar_index * tau */ write_tune_stat = 1; /* how many times do we write tuning stats? */ auto_stop_tune = 1; /* auto stop tuning runs? default: on */ write_llog = 0; /* write local llog files when tuning; default: off */ #endif /* following part parses command line for options and their arguments */ optarg = NULL; while( (c = getopt(argc, argv, OPTS)) != -1 ) { switch(c) { case 'a': accuracy = atof(optarg); if ( accuracy <= 0 ) error("fly_sa: accuracy (%g) is too small", accuracy); break; case 'b': /* -b sets backup frequency (to write state file) */ state_write = strtol(optarg, NULL, 0); if ( state_write < 1 ) error("fly_sa: max. backup frequency is every tau steps i.e. -b 1"); if ( state_write == LONG_MAX ) error("fly_sa: argument for -b too large"); break; case 'B': /* -B runs in benchmark mode -> quit after intial steps */ bench = 1; time_flag = 1; break; case 'c': /* choice of objective function */ setCostType(optarg); break; case 'C': /* parallel code: -C sets the covar sample index */ #ifdef MPI /* -C is currently disabled since I (Yogi) could not get it to work and I * * got tired of debugging and I didn't really need it at that time and so * * on and so forth; try running on a smaller number of nodes or get it to * * work yourself; that's what I say. Grmbl. */ /* covar_index = atoi(optarg); if ( covar_index < 1 ) error("fly_sa: covariation sample index must be >= 1"); */ error("fly_sa: -C does not work yet, try to run on fewer nodes"); #else error("fly_sa: can't use -C in serial, tuning only in parallel"); #endif break; case 'e': /* -e sets the stopping criterion */ if( !(strcmp(optarg,"pfreeze")) ) stop_flag = proportional_freeze; else if( !(strcmp(optarg,"afreeze")) ) stop_flag = absolute_freeze; else if( !(strcmp(optarg,"abs")) ) stop_flag = absolute_energy; else error("fly_sa: valid stopping criteria are pfreeze, afreeze, abs"); break; case 'E': /* define GA type */ ga_type=atoi(optarg); if (ga_type<0 || ga_type>1) error("fly_**: ga type is between [0,1]"); break; case 'f': precision = atoi(optarg); /* -f determines float precision */ if ( precision < 0 ) error("fly_sa: what exactly would a negative precision be???"); if ( precision > MAX_PRECISION ) error("fly_sa: max. float precision is %d!", MAX_PRECISION); break; case 'g': /* -g choose g(u) function */ if( !(strcmp(optarg, "s")) ) gofu = Sqrt; else if ( !(strcmp(optarg, "t"))) gofu = Tanh; else if ( !(strcmp(optarg, "e"))) gofu = Exp; else if ( !(strcmp(optarg, "h"))) gofu = Hvs; else error("fly_sa: %s is an invalid g(u), should be e, h, s or t", optarg); break; case 'h': /* -h help option */ PrintMsg(help, 0); break; case 'i': /* -i sets the stepsize */ stepsize = atof(optarg); if ( stepsize < 0 ) error("fly_sa: going backwards? (hint: check your -i)"); if ( stepsize == 0 ) error("fly_sa: going nowhere? (hint: check your -i)"); if ( stepsize > MAX_STEPSIZE ) error("fly_sa: stepsize %g too large (max. is %g)", stepsize, MAX_STEPSIZE); break; case 'I': IndivOrIntegrated=atoi(optarg); if (IndivOrIntegrated < 0 || IndivOrIntegrated >1) warning("fly_ used integrated data (0) or individual data (1). Default is 0"); break; case 'j': break; case 'l': /* -l displays the log to the screen */ log_flag = 1; break; case 'L': /* -L writes local .llog files when tuning */ ls_type=atoi(optarg); if (ls_type<0 || ls_type>8) error("fly_**: local search type is between [0,8]"); break; case 'n': /* -n suppresses .state and .log files */ nofile_flag = 1; break; case 'N': /* -N sets laNdscape flag and Equilibrate mode */ equil = 1; landscape_flag = 1; break; case 'o': /* -o sets old division style (ndivs = 3 only! ) */ olddivstyle = 1; break; case 'p': /* -p sets prolix printing */ prolix_flag = 1; break; case 'Q': /* -Q sets quenchit mode (serial code only) */ #ifdef MPI error("fly_sa: can't use -Q in parallel, quenchit only in serial"); #else quenchit = 1; #endif break; case 'r': rmsflag = atoi(optarg); /* flag for printing root mean square */ if (rmsflag==0 || rmsflag==1) rmsflag; else error("fly_**: rmsflag=1 or rmsflag=0"); break; case 'R': restart=1; break; case 's': /* -s sets solver to be used */ if ( !(strcmp(optarg, "a")) ) ps = Adams; else if ( !(strcmp(optarg, "bd")) ) ps = BaDe; else if ( !(strcmp(optarg, "bs")) ) ps = BuSt; else if ( !(strcmp(optarg, "e")) ) ps = Euler; else if ( !(strcmp(optarg, "h")) ) ps = Heun; else if ( !(strcmp(optarg, "mi")) || !(strcmp(optarg, "m")) ) ps = Milne; else if ( !(strcmp(optarg, "me")) ) ps = Meuler; else if ( !(strcmp(optarg, "r4")) || !(strcmp(optarg, "r")) ) ps = Rk4; else if ( !(strcmp(optarg, "r2")) ) ps = Rk2; else if ( !(strcmp(optarg, "rck")) ) ps = Rkck; else if ( !(strcmp(optarg, "rf")) ) ps = Rkf; else error("fly_sa: invalid solver (%s), use: a,bs,e,h,mi,me,r{2,4,ck,f} imex", optarg); break; case 'S': /* -S unsets the auto_stop_tune flag */ #ifdef MPI auto_stop_tune = 0; #else error("fly_sa: can't use -S in serial, tuning only in parallel"); #endif break; case 't': /* -t saves times in data and .times files */ time_flag = 1; break; #ifdef MPI tuning = 1; #else error("fly_sa: can't use -T in serial, tuning only in parallel"); #endif break; case 'v': /* -v prints version message */ fprintf(stderr, "%s\n", version); // fprintf(stderr, verstring, USR, MACHINE, COMPILER, FLAGS, __DATE__, __TIME__); exit(0); case 'w': /* -w sets output file */ outDir=(char *)calloc(MAX_RECORD, sizeof(char)); outDir=strcpy(outDir, optarg); break; case 'W': /* -W determines the frequency of writing tuning stats */ #ifdef MPI write_tune_stat = atoi(optarg); if ( write_tune_stat < 1 ) error("fly_sa: frequency of writing tune stats must be >= 1"); #else error("fly_sa: can't use -W in serial, tuning only in parallel"); #endif break; case 'x': section_title = strcpy(section_title, optarg); break; case 'X': ga_config = (char *)calloc(MAX_RECORD, sizeof(char)); ga_config = strcpy(ga_config, optarg); break; case 'y': /* -y set frequency to print log */ print_freq = strtol(optarg, NULL, 0); if ( print_freq < 1 ) error("fly_sa: can't print status more than every tau steps (-y 1)"); if ( print_freq == LONG_MAX ) error("fly_sa: argument for -y too large"); break; case 'z': custom_gast = atof(optarg); if ( custom_gast < 0. ) error("unfold: gastrulation time must be positive"); if ( custom_gast > 10000000. ) error("unfold: gastrulation time must be smaller than 10'000'000"); break; case ':': error("fly_sa: need an argument for option -%c", optopt); break; case '?': default: error("fly_sa: unrecognized option -%c", optopt); } } /* error checking here */ #ifdef MPI if ( (tuning == 1) && (equil == 1) ) error("fly_sa: can't combine -E with -T"); if ( write_llog && !tuning ) error("fly_sa: -L only makes sense when tuning"); #else if ( (quenchit == 1) && (equil == 1) ) error("fly_sa: can't combine -E with -Q"); #endif if ( ((argc - (optind - 1)) != 2) ) PrintMsg(usage, 1); initResCompPar(argv[argc-1]); return optind; } NucStatePtr fly_main(int argc, char **argv){ int opt_index; /* pointer to current argument of command line */ int stateflag = 0; /* state file or not? */ NucStatePtr ga_state; /* global GA parameter struct */ char sep; /* allocate memory for static file names */ inputfile = (char *)calloc(MAX_RECORD, sizeof(char)); statefile = (char *)calloc(MAX_RECORD, sizeof(char)); /* allocate memory for static Lam parameters and dance partners */ ga_state = (NucStateType *)malloc(sizeof(NucStateType)); section_title = (char *)calloc(MAX_RECORD, sizeof(char)); section_title = strcpy(section_title, "input"); /* default is eqparms */ opt_index = ParseCommandLine(argc, argv); inputfile = strcpy(inputfile, argv[opt_index]); #ifdef MPI if ( nnodes>1 ) { if ( nnodes <= 10 ) sprintf(statefile, "%s_%d.state", inputfile, myid); else if ( (nnodes > 10) && (nnodes <= 100) ) sprintf(statefile, "%s_%02d.state", inputfile, myid); else if ( (nnodes > 100) && (nnodes <= 1000) ) sprintf(statefile, "%s_%03d.state", inputfile, myid); else if ( (nnodes > 1000) && (nnodes <= 10000) ) sprintf(statefile, "%s_%04d.state", inputfile, myid); else if ( (nnodes > 10000) && (nnodes <= 100000) ) sprintf(statefile, "%s_%05d.state", inputfile, myid); else error("Initialize: can't open more than 100'000 state files"); } else #endif sprintf(statefile, "%s.state", inputfile); /* check if a state file exists (access() is in unistd.h) */ if ( 0 == access(statefile, F_OK) ){ stateflag = 1; printf("continue with a previous run\n"); } #ifdef MPI /* parallel code: make sure that all state files are present */ MPI_Allreduce(&stateflag, &flagsum, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); if ( (flagsum > 0) && (flagsum != nnodes) && (stateflag == 0) ) error("Initialize: state file for process %d is missing"); #endif /* first get Lam parameters, initial temp and energy and initialize S_0 */ /* if we restore a run from a state file: call RestoreState() */ if ( !stateflag ) { InitialMove(argc, argv, opt_index, ga_state, &energy ); } else{ // RestoreState(statefile, state, &energy); // reInitDistribution(); printf("restore old statefile\n"); } ga_state->tune.programName= (char *)calloc(MAX_RECORD, sizeof(char)); ga_state->tune.programName=argv[0]; ga_state->tune.configName = (char *)calloc(MAX_RECORD, sizeof(char)); if (ga_config!=NULL){ ga_state->tune.configName=ga_config; printf(" config file is : %s\n", ga_state->tune.configName); } else{ ga_state->tune.configName= ""; } return ga_state; } double checkConstraint(){ return checkBound(); } double getFlyScore(){ double chisq ; double ms, rms; chisq =Score(); if (chisq ==FORBIDDEN_MOVE) return FORBIDDEN_MOVE; if (chisq ==OUT_OF_BOUND) // return OUT_OF_BOUND; return FORBIDDEN_MOVE; if ( rmsflag==1 ) { /* print rms */ ms = (chisq - GetCurPenalty()) / (double)ndp; rms = sqrt(ms); return rms; }else return chisq; } double getFlyScoreNoCheck(){ double chisq ; double ms, rms; chisq =ScoreNoCheck(); if (chisq ==FORBIDDEN_MOVE) return FORBIDDEN_MOVE; if (chisq ==OUT_OF_BOUND) return OUT_OF_BOUND; if ( rmsflag==1 ) { /* print rms */ ms = (chisq - GetCurPenalty()) / (double)ndp; rms = sqrt(ms); return rms; }else return chisq; } double getRMSScore(double lse){ double ms, rms; if ( rmsflag ) { /* print rms */ ms = (lse - GetPenalty()) / (double)ndp; rms = sqrt(ms); return lse; }else return lse; } double *getFlyScoreLM(int diffOrNoDiff, int DIForJAC){ double aa; // aa=getFlyScoreNoCheck(); // printf("aa :%f\n",aa); return ScoreNoCheckLM(diffOrNoDiff,DIForJAC); } Range getRange(int ivar){ return *(ptab[ivar].param_range); } double getParameter(int ivar){ SearchSpace *limits; /* local copy of the SearchSpace */ limits = GetLimits(); /* get pointer to SearchSpace in score.c */ if ( limits->pen_vec == 0 ) return *(ptab[ivar].param); else { if ((ptab[ivar].scale)==NULL) return *(ptab[ivar].param); else return *(ptab[ivar].param)**(ptab[ivar].scale); } } double *getParameters(int nbParameters){ double *array=AllocDbl1D(nbParameters); int ivar; SearchSpace *limits; /* local copy of the SearchSpace */ limits = GetLimits(); /* get pointer to SearchSpace in score.c */ for (ivar=0;ivar pen_vec == 0 ) array[ivar] =*(ptab[ivar].param); else { if ((ptab[ivar].scale)==NULL) array[ivar] =*(ptab[ivar].param); else array[ivar] = *(ptab[ivar].param)**(ptab[ivar].scale); } } return array; } double pamInLim(double param, int index){ static SearchSpace *limits; /* local copy of the SearchSpace */ limits = GetLimits(); /* get pointer to SearchSpace in score.c */ if ( limits->pen_vec == 0 ) return param; else { if ((ptab[index].scale)==NULL) return param; else return param**(ptab[index].scale); } } void setParameters(int i, double val){ if ( GetLimits()->pen_vec == 0 ) *(ptab[i].param)= val ; else{ if ((ptab[i].scale)==NULL) *(ptab[i].param)= val ; else *(ptab[i].param)= val/ *(ptab[i].scale); } } void setAllParameter(EqParms newparm){ int i,j,k; int index; index=0; for ( i=0; i < defs.ngenes; i++ ) /* pointers to R stuff */ *(ptab[index+1].param)= newparm.R[i]; for ( i=0; i < defs.ngenes; i++ ) /* pointers to T stuff */ for ( j=0; j < defs.ngenes; j++ ) *(ptab[index++].param) = newparm.T[j+i*defs.ngenes]; for ( i=0; i < defs.ngenes; i++ ) /* pointers to m stuff */ *(ptab[index++].param)= newparm.m[i]; for ( i=0; i < defs.ngenes; i++ ) /* pointers to h stuff */ *(ptab[index++].param) = newparm.h[i]; for ( i=0; i < defs.ngenes; i++ ) /* pointers to d stuff */ *(ptab[index++].param) = newparm.d[i]; if ( (defs.diff_schedule == 'A') || (defs.diff_schedule == 'C') ) i=defs.ngenes+1; for ( i=0; i < defs.ngenes; i++ ) /* pointers to lambda stuff */ *(ptab[index++].param) = newparm.lambda[i]; } void plot(DArrPtr tabtimes){ plotFly(tabtimes); }