openModeller
Version 1.4.0
|
00001 /* ************************************** 00002 * GARP Modeling Package 00003 * 00004 * ************************************** 00005 * 00006 * Copyright (c), The Center for Research, University of Kansas, 2385 Irving Hill Road, Lawrence, KS 66044-4755, USA. 00007 * Copyright (C), David R.B. Stockwell of Symbiotik Pty. Ltd. 00008 * 00009 * This program is free software; you can redistribute it and/or modify 00010 * it under the terms of the license that is distributed with the software. 00011 * 00012 * This program is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * license.txt file included with this software for more details. 00016 */ 00017 00018 // Algorithm.cpp : Implementation of CGarpAlgorithm 00019 00020 #ifdef WIN32 00021 // avoid warnings caused by problems in VC headers 00022 #define _SCL_SECURE_NO_DEPRECATE 00023 #endif 00024 00025 #include "Rule.h" 00026 #include "GarpAlgorithm.h" 00027 #include "Utilities.h" 00028 #include "EnvCellSet.h" 00029 00030 #include <openmodeller/om.hh> 00031 #include <openmodeller/Random.hh> 00032 #include <openmodeller/ScaleNormalizer.hh> 00033 00034 #include <time.h> 00035 00036 #include <string> 00037 00038 #define NUM_PARAM 4 00039 00040 /****************************************************************/ 00041 /*** Algorithm parameter metadata *******************************/ 00042 00043 static AlgParamMetadata parameters[NUM_PARAM] = 00044 { 00045 // Metadata of the first parameter. 00046 { 00047 "MaxGenerations", // Id. 00048 "Max generations", // Name. 00049 Integer, // Type. 00050 00051 // Overview. 00052 "Maximum number of iterations run by the Genetic Algorithm.", 00053 00054 // Description. 00055 "Maximum number of iterations (generations) run by the Genetic Algorithm.", 00056 00057 1, // Not zero if the parameter has lower limit. 00058 1, // Parameter's lower limit. 00059 0, // Not zero if the parameter has upper limit. 00060 0, // Parameter's upper limit. 00061 "400" // Parameter's typical (default) value. 00062 }, 00063 00064 { 00065 "ConvergenceLimit", // Id. 00066 "Convergence limit", // Name. 00067 Real, // Type. 00068 00069 // Overview. 00070 "Defines the convergence value that makes the algorithm stop (before reaching MaxGenerations).", 00071 00072 // Description. 00073 "Defines the convergence value that makes the algorithm stop (before reaching MaxGenerations).", 00074 00075 1, // Not zero if the parameter has lower limit. 00076 0.0, // Parameter's lower limit. 00077 1, // Not zero if the parameter has upper limit. 00078 1.0, // Parameter's upper limit. 00079 "0.01" // Parameter's typical (default) value. 00080 }, 00081 00082 { 00083 "PopulationSize", // Id. 00084 "Population size", // Name. 00085 Integer, // Type. 00086 00087 "Maximum number of rules to be kept in solution.", // Overview. 00088 "Maximum number of rules to be kept in solution.", // Description 00089 00090 1, // Not zero if the parameter has lower limit. 00091 1, // Parameter's lower limit. 00092 1, // Not zero if the parameter has upper limit. 00093 500, // Parameter's upper limit. 00094 "50" // Parameter's typical (default) value. 00095 }, 00096 00097 { 00098 "Resamples", // Id. 00099 "Resamples", // Name. 00100 Integer, // Type. 00101 00102 // Overview. 00103 "Number of points sampled (with replacement) used to test rules.", 00104 00105 // Description. 00106 "Number of points sampled (with replacement) used to test rules.", 00107 00108 1, // Not zero if the parameter has lower limit. 00109 1, // Parameter's lower limit. 00110 1, // Not zero if the parameter has upper limit. 00111 100000, // Parameter's upper limit. 00112 "2500" // Parameter's typical (default) value. 00113 } 00114 }; 00115 00116 00117 /*****************************************************************/ 00118 /*** Algorithm's general metadata ********************************/ 00119 00120 static AlgMetadata metadata = { 00121 00122 "DG_GARP", // Id. 00123 "GARP (single run) - DesktopGARP implementation", // Name. 00124 "1.1 alpha", // Version. 00125 00126 // Overview. 00127 "This is the 2nd implementation of GARP algorithm, based on the \ 00128 original C code by David Stockwell. This version correspondss to \ 00129 the version in use by the DesktopGarp modeling package, with \ 00130 modifications to use OpenModeller base data access objects.", 00131 00132 // Description. 00133 "GARP is a genetic algorithm that creates ecological niche \ 00134 models for species. The models describe environmental conditions \ 00135 under which the species should be able to maintain populations. \ 00136 For input, GARP uses a set of point localities where the species \ 00137 is known to occur and a set of geographic layers representing \ 00138 the environmental parameters that might limit the species' \ 00139 capabilities to survive.", 00140 00141 // Author 00142 "Stockwell, D. R. B., modified by Ricardo Scachetti Pereira", 00143 00144 // Bibliography. 00145 "Stockwell, D. R. B. 1999. Genetic algorithms II. \ 00146 Pages 123-144 in A. H. Fielding, editor. \ 00147 Machine learning methods for ecological applications. \ 00148 Kluwer Academic Publishers, Boston.\ 00149 \n\ 00150 Stockwell, D. R. B., and D. P. Peters. 1999. \ 00151 The GARP modelling system: Problems and solutions to automated \ 00152 spatial prediction. International Journal of Geographic \ 00153 Information Systems 13:143-158.\ 00154 \n\ 00155 Stockwell, D. R. B., and I. R. Noble. 1992. \ 00156 Induction of sets of rules from animal distribution data: \ 00157 A robust and informative method of analysis. Mathematics and \ 00158 Computers in Simulation 33:385-390.", 00159 00160 "Ricardo Scachetti Pereira", // Code author. 00161 "rpereira [at] ku.edu", // Code author's contact. 00162 00163 0, // Does not accept categorical data. 00164 1, // Does not need (pseudo)absence points. 00165 00166 NUM_PARAM, // Algorithm's parameters. 00167 parameters 00168 }; 00169 00170 00171 00172 /****************************************************************/ 00173 /****************** Algorithm's factory function ****************/ 00174 #ifndef DONT_EXPORT_GARP_FACTORY 00175 OM_ALG_DLL_EXPORT 00176 AlgorithmImpl * 00177 algorithmFactory() 00178 { 00179 return new GarpAlgorithm(); 00180 } 00181 00182 OM_ALG_DLL_EXPORT 00183 AlgMetadata const * 00184 algorithmMetadata() 00185 { 00186 return &metadata; 00187 } 00188 #endif 00189 00190 00191 // ========================================================================== 00192 // ========================================================================== 00193 // GarpAlgorithm implementation 00194 // ========================================================================== 00195 GarpAlgorithm::GarpAlgorithm() 00196 : AlgorithmImpl(& metadata) 00197 { 00198 _normalizerPtr = new ScaleNormalizer( 1.0, 253.0, true ); 00199 00200 GarpUtil::randomize(0); 00201 srand((unsigned)time(NULL)); 00202 00203 // initialize properties to default values 00204 initializeProperties(); 00205 00206 // data points used for training 00207 objTrainSet = NULL; 00208 //Log::instance()->info("GarpAlgorithm::GarpAlgorithm() at %x\n", this ); 00209 } 00210 00211 // ========================================================================== 00212 GarpAlgorithm::~GarpAlgorithm() 00213 { 00214 //objBest.log(); 00215 00216 if (objTrainSet) 00217 delete objTrainSet; 00218 00219 //Log::instance()->info("GarpAlgorithm::~GarpAlgorithm() at %x\n", this ); 00220 } 00221 00222 00223 // ========================================================================== 00224 int GarpAlgorithm::initialize() 00225 { 00226 if (!getParameter("MaxGenerations", &Totalgens)) { 00227 Log::instance()->error("Parameter MaxGenerations not set properly."); 00228 return 0; 00229 } 00230 00231 if (!getParameter("ConvergenceLimit", &Conv_limit)) { 00232 Log::instance()->error("Parameter ConvergenceLimit not set properly."); 00233 return 0; 00234 } 00235 00236 if (!getParameter("PopulationSize", &Popsize)) { 00237 Log::instance()->error("Parameter PopulationSize not set properly."); 00238 return 0; 00239 } 00240 00241 if (!getParameter("Resamples", &Resamples)) { 00242 Log::instance()->error("Parameter Resamples not set properly."); 00243 return 0; 00244 } 00245 00246 //Log::instance()->debug("MaxGenerations set to: %d\n", _max_gen); 00247 //Log::instance()->debug("ConvergenceLimit set to: %.4f\n", _conv_limit); 00248 //Log::instance()->debug("PopulationSize set to: %d\n", _popsize); 00249 //Log::instance()->debug("Resamples set to: %d\n", _resamples); 00250 00251 int dim = _samp->numIndependent(); 00252 iActiveGenes = dim + 1; 00253 objBest.setDimension(dim + 1); 00254 objNew.setDimension(dim + 1); 00255 00256 //printf("Dimensions: %d\n", dim); 00257 00258 //printf("Presences: %d\nAbsences: %d\n", _samp->numPresence(), _samp->numAbsence()); 00259 00260 objTrainSet = new EnvCellSet; 00261 objTrainSet->initialize(Resamples); 00262 for (int i = 0; i < Resamples; i++) 00263 { 00264 OccurrencePtr oc = _samp->getOneSample(); 00265 Sample sample = (*oc).environment(); 00266 Scalar dep = ( (*oc).abundance() > 0.0 ) ? 1.0 : 0.0; 00267 00268 // transfer values to cell 00269 BYTE * values = new BYTE[dim + 2]; 00270 EnvCell * cell = new EnvCell(dim + 2, values); 00271 values[0] = (BYTE) dep; 00272 //printf("%5d]: dep=[%d] ", i, values[0]); 00273 for (int j = 0; j < dim; j++) 00274 { 00275 values[j + 1] = (BYTE) sample[j]; 00276 //printf("%3d (%8.3f) ", values[j + 1], sample[j]); 00277 } 00278 // Always initialize last element (added by rdg 2009-01-09) 00279 values[dim+1] = (BYTE)0; 00280 00281 //printf("\n"); 00282 objTrainSet->add(cell); 00283 } 00284 00285 // create histograms of occurrence of each layer value for bioclim/range rules 00286 objTrainSet->createBioclimHistogram(); 00287 00288 // get initial model 00289 getInitialModel(Popsize, objTrainSet); 00290 00291 return 1; 00292 } 00293 00294 00295 /****************************************************************/ 00296 /****************** getProgress *********************************/ 00297 00298 float GarpAlgorithm::getProgress() const 00299 { 00300 if (done()) 00301 { return 1.0; } 00302 else 00303 { 00304 float byIterations = ( Gen / (float) Totalgens ); 00305 float byConvergence = (float)( Conv_limit / Convergence ); 00306 float progress = (byIterations > byConvergence) ? byIterations : byConvergence; 00307 00308 if (progress > _maxProgress) 00309 { _maxProgress = progress; } 00310 return _maxProgress; 00311 } 00312 } 00313 00314 00315 /****************************************************************/ 00316 /****************** getConvergence ******************************/ 00317 00318 int GarpAlgorithm::getConvergence( Scalar * const val ) const 00319 { 00320 *val = Convergence; 00321 return 0; 00322 } 00323 00324 00325 /****************************************************************/ 00326 /****************** getValue ************************************/ 00327 00328 Scalar GarpAlgorithm::getValue( const Sample& x ) const 00329 { 00330 return objBest.getValue(x); 00331 } 00332 00333 // ========================================================================== 00334 void GarpAlgorithm::initializeProperties() 00335 { 00336 int i, j; 00337 00338 lVersion = 0; 00339 00340 iCPUTime = 0; 00341 00342 Improvements = 0; 00343 00344 Resamples = 2500; 00345 Accuracylimit = 0.0; // The minimium post prob of a rule 00346 MinUsage = 0.00; 00347 Mortality = 0.9; 00348 00349 Experiment = 0; // generations per experiment 00350 Totalgens = 2560; // generations per experiment 00351 Totaltrials = 0; // trials per experiment 00352 Popsize = 50; // population size 00353 C_rate = 0.25; // crossover rate 00354 M_rate = 0.50; // mutation rate 00355 J_rate = 0.25; // join rate 00356 I_rate = 0.0; 00357 Gapsize = 0.8; // fraction of pop generated from archive 00358 Maxspin = 2; // max gens without evals 00359 Resampling_f = 1.0 ; 00360 Significance = 2.70 ; // the minimum level for inc in best 00361 Conv_limit = 0.10 ; // the fractional addition of rules 00362 Cutval = 0.0; 00363 Trials = 0; 00364 00365 // data collection and loop control variables 00366 Ave_current_perf = 0.0; // ave perf in current generation 00367 Best_current_perf = 0.0; // best perf in current generation 00368 Worst_current_perf = 0.0; // worst perf in current generation 00369 Best = 0.0; // best performance seen so far 00370 Worst = 0.0; // worst performance seen so far 00371 Best_guy = 0; // index of best_current_perf 00372 00373 Conv = 0; // number of partially coverged genes 00374 Doneflag = false;// set when termination conditions hold 00375 Gen = 0; // generation counter 00376 Lost = 0; // number of totally coverged positions 00377 Spin = 0; // number of gens since eval occurred 00378 Convergence = 1.0; // the stability of rule set 00379 Improvements = 0;// the number of times the best set has been 00380 // improved by addition or alteration 00381 Resample=1; 00382 00383 for (i = 0; i < 2; i++) 00384 for (j = 0; j < 5; j++) 00385 Heuristic[i][j] = 0; // successes of heuristic operators 00386 00387 // flags 00388 Sigflag = 1; 00389 Postflag = 0; 00390 Compflag = 0; 00391 Adjustflag = 1; 00392 00393 BioclimOnlyFlag = 0; 00394 LogitOnlyFlag = 0; 00395 RangeRuleFlag = 1; 00396 NegatedRuleFlag = 1; 00397 AtomicRuleFlag = 1; 00398 LogitRuleFlag = 1; 00399 00400 00401 // reset gene selection 00402 // gene 0 is always active (it is the species gene) 00403 iActiveGenes = 1; 00404 iGeneIndex[0] = 0; 00405 bGeneIsActive[0] = true; 00406 for (i = 1; i < MAX_ENV_LAYERS; i++) 00407 { 00408 bGeneIsActive[i] = true; 00409 iGeneIndex[i] = i; 00410 } 00411 00412 _maxProgress = 0.0; 00413 } 00414 00415 // ========================================================================== 00416 char * GarpAlgorithm::getParameter2(char * strParamName) 00417 { 00418 const int iStringLen = 512; 00419 char * strResult = new char[iStringLen]; 00420 00421 strcpy(strResult, ""); 00422 00423 // most used parameters 00424 if (strcmp(strParamName, "Gen" ) == 0) sprintf(strResult, "%d", Gen); 00425 00426 // flags 00427 else if (strcmp(strParamName, "Sigflag" ) == 0) sprintf(strResult, "%d", Sigflag); 00428 else if (strcmp(strParamName, "Postflag" ) == 0) sprintf(strResult, "%d", Postflag); 00429 else if (strcmp(strParamName, "Compflag" ) == 0) sprintf(strResult, "%d", Compflag); 00430 else if (strcmp(strParamName, "Adjustflag" ) == 0) sprintf(strResult, "%d", Adjustflag); 00431 else if (strcmp(strParamName, "BioclimOnlyFlag" ) == 0) sprintf(strResult, "%d", BioclimOnlyFlag); 00432 else if (strcmp(strParamName, "LogitOnlyFlag" ) == 0) sprintf(strResult, "%d", LogitOnlyFlag); 00433 else if (strcmp(strParamName, "RangeRuleFlag" ) == 0) sprintf(strResult, "%d", RangeRuleFlag); 00434 else if (strcmp(strParamName, "NegatedRuleFlag" ) == 0) sprintf(strResult, "%d", NegatedRuleFlag); 00435 else if (strcmp(strParamName, "AtomicRuleFlag" ) == 0) sprintf(strResult, "%d", AtomicRuleFlag); 00436 else if (strcmp(strParamName, "LogitRuleFlag" ) == 0) sprintf(strResult, "%d", LogitRuleFlag); 00437 00438 // int parameters 00439 else if (strcmp(strParamName, "Version" ) == 0) sprintf(strResult, "%d", static_cast<int>(lVersion)); 00440 else if (strcmp(strParamName, "CPUTime" ) == 0) sprintf(strResult, "%d", iCPUTime); 00441 else if (strcmp(strParamName, "Resamples" ) == 0) sprintf(strResult, "%d", Resamples); 00442 else if (strcmp(strParamName, "Totalgens" ) == 0) sprintf(strResult, "%d", Totalgens); 00443 else if (strcmp(strParamName, "Totaltrials" ) == 0) sprintf(strResult, "%d", Totaltrials); 00444 else if (strcmp(strParamName, "Popsize" ) == 0) sprintf(strResult, "%d", Popsize); 00445 else if (strcmp(strParamName, "Maxspin" ) == 0) sprintf(strResult, "%d", Maxspin); 00446 else if (strcmp(strParamName, "Best_guy" ) == 0) sprintf(strResult, "%d", Best_guy); 00447 else if (strcmp(strParamName, "Resample" ) == 0) sprintf(strResult, "%d", Resample); 00448 else if (strcmp(strParamName, "Conv" ) == 0) sprintf(strResult, "%d", Conv); 00449 else if (strcmp(strParamName, "Trials" ) == 0) sprintf(strResult, "%d", Trials); 00450 else if (strcmp(strParamName, "Experiment" ) == 0) sprintf(strResult, "%d", Experiment); 00451 else if (strcmp(strParamName, "Lost" ) == 0) sprintf(strResult, "%d", Lost); 00452 else if (strcmp(strParamName, "Spin" ) == 0) sprintf(strResult, "%d", Spin); 00453 else if (strcmp(strParamName, "Improvements" ) == 0) sprintf(strResult, "%d", Improvements); 00454 else if (strcmp(strParamName, "Doneflag" ) == 0) sprintf(strResult, "%d", (int) Doneflag); 00455 else if (strcmp(strParamName, "Heuristic_0" ) == 0) sprintf(strResult, "%d", Heuristic[0][0]); 00456 else if (strcmp(strParamName, "Heuristic_1" ) == 0) sprintf(strResult, "%d", Heuristic[0][1]); 00457 else if (strcmp(strParamName, "Heuristic_2" ) == 0) sprintf(strResult, "%d", Heuristic[0][2]); 00458 else if (strcmp(strParamName, "Heuristic_3" ) == 0) sprintf(strResult, "%d", Heuristic[0][3]); 00459 else if (strcmp(strParamName, "Heuristic_4" ) == 0) sprintf(strResult, "%d", Heuristic[0][4]); 00460 00461 // double parameters 00462 else if (strcmp(strParamName, "Accuracylimit" ) == 0) sprintf(strResult, "%f", Accuracylimit); 00463 else if (strcmp(strParamName, "MinUsage" ) == 0) sprintf(strResult, "%f", MinUsage); 00464 else if (strcmp(strParamName, "Mortality" ) == 0) sprintf(strResult, "%f", Mortality); 00465 else if (strcmp(strParamName, "C_rate" ) == 0) sprintf(strResult, "%f", C_rate); 00466 else if (strcmp(strParamName, "M_rate" ) == 0) sprintf(strResult, "%f", M_rate); 00467 else if (strcmp(strParamName, "J_rate" ) == 0) sprintf(strResult, "%f", J_rate); 00468 else if (strcmp(strParamName, "I_rate" ) == 0) sprintf(strResult, "%f", I_rate); 00469 else if (strcmp(strParamName, "Gapsize" ) == 0) sprintf(strResult, "%f", Gapsize); 00470 else if (strcmp(strParamName, "Resampling_f" ) == 0) sprintf(strResult, "%f", Resampling_f); 00471 else if (strcmp(strParamName, "Significance" ) == 0) sprintf(strResult, "%f", Significance); 00472 else if (strcmp(strParamName, "Conv_limit" ) == 0) sprintf(strResult, "%f", Conv_limit); 00473 else if (strcmp(strParamName, "Cutval" ) == 0) sprintf(strResult, "%f", Cutval); 00474 else if (strcmp(strParamName, "Ave_current_perf" ) == 0) sprintf(strResult, "%f", Ave_current_perf); 00475 else if (strcmp(strParamName, "Best_current_perf" ) == 0) sprintf(strResult, "%f", Best_current_perf); 00476 else if (strcmp(strParamName, "Worst_current_perf" ) == 0) sprintf(strResult, "%f", Worst_current_perf); 00477 else if (strcmp(strParamName, "Best" ) == 0) sprintf(strResult, "%f", Best); 00478 else if (strcmp(strParamName, "Worst" ) == 0) sprintf(strResult, "%f", Worst); 00479 else if (strcmp(strParamName, "Convergence" ) == 0) sprintf(strResult, "%f", Convergence); 00480 00481 // special parameters 00482 else if (strcmp(strParamName, "SelectedLayers" ) == 0) 00483 { 00484 char * strAux = getSelectedLayersAsString(); 00485 sprintf(strResult, "%s", strAux); 00486 delete strAux; 00487 } 00488 00489 if (strlen(strResult) > static_cast<unsigned int>(iStringLen)) 00490 throw GarpException(82, "String size exceeded in getParameter::parametersToXML()"); 00491 00492 return strResult; 00493 } 00494 00495 // ========================================================================== 00496 void GarpAlgorithm::setParameter(char * strParamName, char * strParamValue) 00497 { 00498 // flags 00499 if (strcmp(strParamName, "Sigflag" ) == 0) Sigflag = atoi(strParamValue); 00500 else if (strcmp(strParamName, "Postflag" ) == 0) Postflag = atoi(strParamValue); 00501 else if (strcmp(strParamName, "Compflag" ) == 0) Compflag = atoi(strParamValue); 00502 else if (strcmp(strParamName, "Adjustflag" ) == 0) Adjustflag = atoi(strParamValue); 00503 else if (strcmp(strParamName, "BioclimOnlyFlag" ) == 0) BioclimOnlyFlag = atoi(strParamValue); 00504 else if (strcmp(strParamName, "LogitOnlyFlag" ) == 0) LogitOnlyFlag = atoi(strParamValue); 00505 else if (strcmp(strParamName, "RangeRuleFlag" ) == 0) RangeRuleFlag = atoi(strParamValue); 00506 else if (strcmp(strParamName, "NegatedRuleFlag" ) == 0) NegatedRuleFlag = atoi(strParamValue); 00507 else if (strcmp(strParamName, "AtomicRuleFlag" ) == 0) AtomicRuleFlag = atoi(strParamValue); 00508 else if (strcmp(strParamName, "LogitRuleFlag" ) == 0) LogitRuleFlag = atoi(strParamValue); 00509 00510 // int parameters 00511 else if (strcmp(strParamName, "Version" ) == 0) lVersion = atoi(strParamValue); 00512 else if (strcmp(strParamName, "CPUTime" ) == 0) iCPUTime = atoi(strParamValue); 00513 else if (strcmp(strParamName, "Resamples" ) == 0) Resamples = atoi(strParamValue); 00514 else if (strcmp(strParamName, "Totalgens" ) == 0) Totalgens = atoi(strParamValue); 00515 else if (strcmp(strParamName, "Totaltrials" ) == 0) Totaltrials = atoi(strParamValue); 00516 else if (strcmp(strParamName, "Popsize" ) == 0) Popsize = atoi(strParamValue); 00517 else if (strcmp(strParamName, "Maxspin" ) == 0) Maxspin = atoi(strParamValue); 00518 else if (strcmp(strParamName, "Best_guy" ) == 0) Best_guy = atoi(strParamValue); 00519 else if (strcmp(strParamName, "Resample" ) == 0) Resample = atoi(strParamValue); 00520 else if (strcmp(strParamName, "Conv" ) == 0) Conv = atoi(strParamValue); 00521 else if (strcmp(strParamName, "Trials" ) == 0) Trials = atoi(strParamValue); 00522 else if (strcmp(strParamName, "Experiment" ) == 0) Experiment = atoi(strParamValue); 00523 else if (strcmp(strParamName, "Gen" ) == 0) Gen = atoi(strParamValue); 00524 else if (strcmp(strParamName, "Lost" ) == 0) Lost = atoi(strParamValue); 00525 else if (strcmp(strParamName, "Spin" ) == 0) Spin = atoi(strParamValue); 00526 else if (strcmp(strParamName, "Improvements" ) == 0) Improvements = atoi(strParamValue); 00527 else if (strcmp(strParamName, "Doneflag" ) == 0) Doneflag = ( atoi(strParamValue) ) ? true : false; 00528 else if (strcmp(strParamName, "Heuristic_0" ) == 0) Heuristic[0][0] = atoi(strParamValue); 00529 else if (strcmp(strParamName, "Heuristic_1" ) == 0) Heuristic[0][1] = atoi(strParamValue); 00530 else if (strcmp(strParamName, "Heuristic_2" ) == 0) Heuristic[0][2] = atoi(strParamValue); 00531 else if (strcmp(strParamName, "Heuristic_3" ) == 0) Heuristic[0][3] = atoi(strParamValue); 00532 else if (strcmp(strParamName, "Heuristic_4" ) == 0) Heuristic[0][4] = atoi(strParamValue); 00533 00534 // double parameters 00535 else if (strcmp(strParamName, "Accuracylimit" ) == 0) Accuracylimit = atof(strParamValue); 00536 else if (strcmp(strParamName, "MinUsage" ) == 0) MinUsage = atof(strParamValue); 00537 else if (strcmp(strParamName, "Mortality" ) == 0) Mortality = atof(strParamValue); 00538 else if (strcmp(strParamName, "C_rate" ) == 0) C_rate = atof(strParamValue); 00539 else if (strcmp(strParamName, "M_rate" ) == 0) M_rate = atof(strParamValue); 00540 else if (strcmp(strParamName, "J_rate" ) == 0) J_rate = atof(strParamValue); 00541 else if (strcmp(strParamName, "I_rate" ) == 0) I_rate = atof(strParamValue); 00542 else if (strcmp(strParamName, "Gapsize" ) == 0) Gapsize = atof(strParamValue); 00543 else if (strcmp(strParamName, "Resampling_f" ) == 0) Resampling_f = atof(strParamValue); 00544 else if (strcmp(strParamName, "Significance" ) == 0) Significance = atof(strParamValue); 00545 else if (strcmp(strParamName, "Conv_limit" ) == 0) Conv_limit = atof(strParamValue); 00546 else if (strcmp(strParamName, "Cutval" ) == 0) Cutval = atof(strParamValue); 00547 else if (strcmp(strParamName, "Ave_current_perf" ) == 0) Ave_current_perf = atof(strParamValue); 00548 else if (strcmp(strParamName, "Best_current_perf" ) == 0) Best_current_perf = atof(strParamValue); 00549 else if (strcmp(strParamName, "Worst_current_perf" ) == 0) Worst_current_perf = atof(strParamValue); 00550 else if (strcmp(strParamName, "Best" ) == 0) Best = atof(strParamValue); 00551 else if (strcmp(strParamName, "Worst" ) == 0) Worst = atof(strParamValue); 00552 else if (strcmp(strParamName, "Convergence" ) == 0) Convergence = atof(strParamValue); 00553 00554 // special parameters 00555 else if (strcmp(strParamName, "SelectedLayers" ) == 0) setSelectedLayers(strParamValue); 00556 } 00557 00558 // ========================================================================== 00559 char * GarpAlgorithm::getSelectedLayersAsString() 00560 { 00561 const int iStringSize = 1024; 00562 char strSelectedLayers[iStringSize], strNextGene[16], * strResult; 00563 00564 strcpy(strSelectedLayers, ""); 00565 for (int i = 1; i < iActiveGenes; i++) 00566 { 00567 sprintf(strNextGene, ";%d", iGeneIndex[i] - 1); 00568 strcat(strSelectedLayers, strNextGene); 00569 } 00570 00571 if (strlen(strSelectedLayers) > static_cast<unsigned int>(iStringSize)) 00572 throw GarpException(82, "String size exceeded in getParameter::parametersToXML()"); 00573 00574 strResult = new char[strlen(strSelectedLayers) + 2]; 00575 // get rid of the first colon (;) if there is at least one gene active 00576 if (iActiveGenes > 1) 00577 strcpy(strResult, strSelectedLayers + 1); 00578 else 00579 strcpy(strResult, strSelectedLayers); 00580 00581 return strResult; 00582 } 00583 00584 // ========================================================================== 00585 void GarpAlgorithm::setSelectedLayers(char * strParamValue) 00586 { 00587 // set active genes, given a string with the gene numbers separeted by colons 00588 char strAux[1024], * strToken; 00589 int iCurrentGene; 00590 00591 iActiveGenes = 1; 00592 iGeneIndex[0] = 0; 00593 bGeneIsActive[0] = true; 00594 for (int i = 1; i < MAX_ENV_LAYERS; i++) 00595 bGeneIsActive[i] = false; 00596 00597 strcpy(strAux, strParamValue); 00598 strToken = strtok(strAux, ";"); 00599 while (strToken != NULL) 00600 { 00601 // get the gene number 00602 iCurrentGene = 1 + atoi(strToken); 00603 bGeneIsActive[iCurrentGene] = true; 00604 iGeneIndex[iActiveGenes] = iCurrentGene; 00605 iActiveGenes++; 00606 00607 strToken = strtok(NULL, ";"); 00608 } 00609 } 00610 00611 // ========================================================================== 00612 int GarpAlgorithm::iterate() 00613 { 00614 if (done()) 00615 return 1; 00616 00617 generate(objTrainSet); 00618 00619 return 1; 00620 } 00621 00622 // ========================================================================== 00623 int GarpAlgorithm::done() const 00624 { 00625 return Doneflag = ( Doneflag || 00626 (Gen >= Totalgens) || 00627 (Spin >= Maxspin) || 00628 (Convergence < Conv_limit) ); 00629 } 00630 00631 // ========================================================================== 00632 void GarpAlgorithm::generate(EnvCellSet * objTrainSet) 00633 { 00634 int i, iNewSize; 00635 00636 // create a new population 00637 Spin++; 00638 00639 if (!Gen || Resample) 00640 objTrainSet->resampleInPlace(); 00641 00642 // evaluate the newly formed population 00643 evaluate(&objNew, objTrainSet); 00644 00645 // put rule into archive 00646 iNewSize = objNew.size(); 00647 for(i = 0; i < iNewSize; i++) 00648 Conv += saveRule(i); 00649 00650 objBest.trim(Popsize); 00651 00652 // gather performance statistics 00653 measure(); 00654 00655 // check termination condition for this experiment 00656 Doneflag = ( done() ) ? true : false; 00657 00658 //objBest.gatherRuleSetStats(Gen); 00659 //objNew.gatherRuleSetStats(-Gen); 00660 00661 if (Doneflag) 00662 { 00663 // sort rules by usage 00664 //objBest.sort(9); 00665 00666 // discard rules which performance is worse than the specified significance 00667 objBest.discardRules(8, Significance); 00668 } 00669 else 00670 { 00671 if (objBest.size()) 00672 select(); 00673 00674 // Fill the rest with new rules */ 00675 colonize(&objNew, objTrainSet, Popsize - objNew.size()); 00676 00677 // Kill off proportion of archive 00678 objBest.trim((int) ((double) objBest.size() * Mortality)); 00679 00680 // apply heuristic operators 00681 mutate(); 00682 crossover(); 00683 join(); 00684 } 00685 00686 // one more generation has been computed 00687 Gen++; 00688 } 00689 00690 // ========================================================================== 00691 void GarpAlgorithm::measure() 00692 { 00693 double performance, ch; 00694 int i, iBestSize; 00695 00696 iBestSize = objBest.size(); 00697 00698 // update current statistics 00699 if (iBestSize > 0) 00700 { 00701 performance = objBest.get(0)->dblPerformance[0]; 00702 Best_current_perf = performance; 00703 Worst_current_perf = performance; 00704 Ave_current_perf = performance; 00705 Best_guy = 0; 00706 } 00707 00708 for (i = 1; i < iBestSize; i++) 00709 { 00710 // update current statistics 00711 performance = objBest.get(i)->dblPerformance[0]; 00712 00713 Ave_current_perf += performance; 00714 if (BETTER(performance, Best_current_perf)) 00715 { 00716 Best_current_perf = performance; 00717 Best_guy = i; 00718 } 00719 00720 if (BETTER(Worst_current_perf, performance)) 00721 Worst_current_perf = performance; 00722 } 00723 00724 Ave_current_perf /= Popsize; 00725 00726 // Adjuct heuristic operators 00727 if (Adjustflag) 00728 { 00729 ch = (double) (Heuristic[1][0] - Heuristic[0][0]); 00730 00731 if (ch < 1) 00732 I_rate = Gapsize = M_rate = 0.1; 00733 else 00734 { 00735 I_rate = ch / (double) objBest.size(); 00736 Gapsize = 0.5 * (Heuristic[1][1] - Heuristic[0][1]) / ch + 0.1; 00737 M_rate = 0.5 * (Heuristic[1][2] - Heuristic[0][2]) / ch + 0.1; 00738 C_rate = C_rate ? 0.5 * (Heuristic[1][3] - Heuristic[0][3]) / ch + 0.1 : 0; 00739 J_rate = J_rate ? 0.5 * (Heuristic[1][4] - Heuristic[0][4]) / ch + 0.1 : 0; 00740 } 00741 00742 for (i = 0; i < 5; i++) 00743 Heuristic[0][i] = Heuristic[1][i]; 00744 00745 //printf("gs=%4.2f ", Gapsize); 00746 } 00747 00748 // update overall performance measures 00749 Convergence = converge(); 00750 00751 //DisplayStatus(); 00752 } 00753 00754 // ========================================================================== 00755 void GarpAlgorithm::DisplayStatus() 00756 { 00757 FILE * fp = fopen("status.txt", "a"); 00758 00759 double perf = objBest.getOveralPerformance(8, 5); 00760 fprintf(fp, "%5d %f\n", Gen, perf); 00761 00762 /* 00763 int i; 00764 00765 fprintf(fp,"Max generations %d", Totalgens); 00766 fprintf(fp,"\nGens \tTrials \tConv \t\tData\tResamp\tBest \tAverage"); 00767 fprintf(fp,"\n%4d \t%6d \t%4f \t%4d \t%4d \t%5.2f \t%5.2f", 00768 Gen, Trials, Convergence, objTrainSet->count(), objTrainSet->count(), 00769 Best, Ave_current_perf); 00770 fprintf(fp,"\nHeuristic: Rules Improve Random Mutate Cross Join"); 00771 fprintf(fp,"\nNo\t\t%d", objBest.size()); 00772 for (i=0; i<5; i++) 00773 fprintf(fp,"\t%d", Heuristic[1][i]); 00774 00775 fprintf(fp,"\nRate\t\t\t%4.2f\t%4.2f\t%4.2f\t%4.2f\t%4.2f", I_rate, Gapsize, M_rate, C_rate, J_rate); 00776 00777 if (objBest.size() > 0) 00778 { 00779 //fprintf(fp,"\n<RuleSet>\n"); 00780 for (i = 0; i < objBest.size(); i++) 00781 { 00782 Rule * pRule = objBest.objRules[i]; 00783 00784 if (pRule->dblPerformance[8] > Significance) 00785 fprintf(fp,"%f %5d %f %d %f %f %d\n", pRule->_pXYs, pRule->_no, pRule->_dA, -1, pRule->_dSig, pRule->dblPerformance[8], pRule->Gene[0] == pRule->intConclusion); 00786 //fprintf(fp,"%s", objBest.objRules[i]->headerToXML()); 00787 } 00788 00789 //fprintf(fp,"</RuleSet>"); 00790 } 00791 00792 fprintf(fp,"\nPerformance: %5.3f",objBest.objRules[0]->dblPerformance[0]); 00793 00794 double * Utility; 00795 00796 Utility = objBest.objRules[0]->dblPerformance; 00797 fprintf(fp,"\n\n\t\tStrength\tCertainty\tError\nAll"); 00798 for (i=1; i<4; i++) 00799 fprintf(fp,"\t\t%5.3f",Utility[i]); 00800 fprintf(fp,"\nSelect"); 00801 for (i=4; i<7; i++) 00802 fprintf(fp,"\t\t%5.3f",Utility[i]); 00803 fprintf(fp,"\n\n\t\tno/n\t\tSignificance\tError\n"); 00804 for (i=7; i<10; i++) 00805 fprintf(fp,"\t\t%5.3f",Utility[i]); 00806 */ 00807 00808 fclose(fp); 00809 } 00810 00811 // ========================================================================== 00812 void GarpAlgorithm::updateHeuOpPerformance(char chrType) 00813 { 00814 Improvements++; 00815 00816 Heuristic[1][0] = Improvements; 00817 00818 switch (chrType) 00819 { 00820 case 'r': Heuristic[1][1]++; break; 00821 case 'm': Heuristic[1][2]++; break; 00822 case 'c': Heuristic[1][3]++; break; 00823 case 'j': Heuristic[1][4]++; break; 00824 } 00825 } 00826 00827 // ========================================================================== 00828 void GarpAlgorithm::colonize(RuleSet * objRules, EnvCellSet * objTrainSet, int intNewRules) 00829 { 00830 int i, p, t; 00831 Rule * newRule=0; 00832 00833 // number of available rule types 00834 int iRuleIndex[4]; 00835 int iRuleTypes; 00836 00837 iRuleTypes = 0; 00838 if (RangeRuleFlag) iRuleIndex[iRuleTypes++] = 0; 00839 if (NegatedRuleFlag) iRuleIndex[iRuleTypes++] = 1; 00840 if (AtomicRuleFlag) iRuleIndex[iRuleTypes++] = 2; 00841 if (LogitRuleFlag) iRuleIndex[iRuleTypes++] = 3; 00842 00843 for (i = 0; i < intNewRules; i++) 00844 { 00845 // pick the next rule to be generate 00846 p = i % iRuleTypes; 00847 t = iRuleIndex[p]; 00848 00849 switch (t) 00850 { 00851 case 0: 00852 newRule = new RangeRule(); 00853 break; 00854 00855 case 1: 00856 newRule = new NegatedRangeRule(); 00857 break; 00858 00859 case 2: 00860 newRule = new AtomicRule(); 00861 break; 00862 00863 case 3: 00864 newRule = new LogitRule(); 00865 break; 00866 } 00867 00868 //printf("Active Genes: %3d\n", iActiveGenes); 00869 00870 newRule->initialize(objTrainSet, objRules, bGeneIsActive, iGeneIndex, iActiveGenes); 00871 newRule->lId = Gen * 1000 + i; 00872 newRule->iOrigGen = Gen; 00873 00874 objRules->add(newRule); 00875 } 00876 } 00877 00878 // ========================================================================== 00879 void GarpAlgorithm::evaluate(RuleSet * objRules, EnvCellSet * objTrainSet) 00880 { 00881 Rule * pRule; 00882 register double performance=0.0; 00883 register int i, n; 00884 00885 Conv = 0; 00886 00887 n = objRules->size(); 00888 for (i = 0; i < n; i++) 00889 { 00890 pRule = objRules->get(i); 00891 00892 if (pRule->needsEvaluation()) 00893 { 00894 performance = pRule->testWithData(objTrainSet); 00895 00896 /* Make not rule if significance is decreased */ 00897 /*if (Speciesflag && performance<0) 00898 { 00899 New[i]->Perf[0] = -New[i]->Perf[0]; 00900 New[i]->Perf[8] = -New[i]->Perf[8]; 00901 New[i]->Type = '!'; 00902 }*/ 00903 00904 pRule->intGens += 1; 00905 pRule->intTrials += 1; 00906 pRule->blnNeedsEvaluation = false; 00907 Trials++; 00908 } 00909 00910 Spin = 0; /* we're making progress */ 00911 00912 if (Trials == 1) 00913 Best = performance; 00914 00915 if (BETTER(performance, Best)) 00916 Best = performance; 00917 } 00918 } 00919 00920 // ========================================================================== 00921 int GarpAlgorithm::saveRule(int iIndex) 00922 { 00923 // Save the ith structure in current population 00924 // if it is one of the Popsize best seen so far 00925 00926 int j, l; // loop control var 00927 int found, ind, iBestSize; 00928 char cNewType; 00929 Rule * Temp, * oNewRule; 00930 00931 // Set criteria for optimizing 00932 ind = 0; 00933 00934 // get best size 00935 iBestSize = objBest.size(); 00936 00937 // get rule from objNew being inserted into objBest 00938 oNewRule = objNew.get(iIndex); 00939 cNewType = oNewRule->type(); 00940 00941 // get Gene length 00942 oNewRule->intGenes * 2; 00943 00944 // Check if an identical or more general structure is already there 00945 for (j = 0, found = 0; j < iBestSize && (!found); j++) 00946 { 00947 if (cNewType == objBest.get(j)->type()) 00948 found = oNewRule->similar(objBest.get(j)); 00949 } 00950 00951 if (found) 00952 { 00953 j--; 00954 00955 //Rule * oRuleToBeReplaced = objBest.get(j); 00956 00957 if (oNewRule->dblPerformance[ind] > objBest.objRules[j]->dblPerformance[ind]) 00958 // Replace if better 00959 { 00960 delete objBest.objRules[j]; 00961 objBest.objRules[j] = oNewRule->clone(); 00962 00963 /* 00964 int k; 00965 for (k = 0; k < iLength; k++) 00966 oRuleToBeReplaced->Gene[k] = oNewRule->Gene[k]; 00967 00968 for (k = 0; k < 10; k++) 00969 oRuleToBeReplaced->dblPerformance[k] = oNewRule->dblPerformance[k]; 00970 00971 oRuleToBeReplaced->intGens = oNewRule->intGens; 00972 oRuleToBeReplaced->chrOrigin = oNewRule->chrOrigin; 00973 oRuleToBeReplaced->intTrials++; 00974 */ 00975 00976 objBest.objRules[j]->intTrials++; 00977 updateHeuOpPerformance(oNewRule->chrOrigin); 00978 } 00979 00980 return 1; 00981 } 00982 00983 00984 // No similar structure found 00985 // allocate new structure at the end if room 00986 00987 Rule * oRuleBeingInserted; 00988 if (iBestSize < Popsize ) 00989 { 00990 // create dummy rule and insert it into the best rule set 00991 oRuleBeingInserted = oNewRule->clone(); 00992 oRuleBeingInserted->dblPerformance[ind] = 0.0; 00993 objBest.add(oRuleBeingInserted); 00994 } 00995 00996 // get best size again 00997 iBestSize = objBest.size(); 00998 00999 // find insertion point and sort j 01000 for (j = 0; (iBestSize > 1 && j < iBestSize - 2) && 01001 (objBest.get(j)->dblPerformance[ind] > oNewRule->dblPerformance[ind]); j++) 01002 { 01003 if (objBest.get(j)->dblPerformance[ind] < objBest.get(j + 1)->dblPerformance[ind]) 01004 { 01005 Temp = objBest.objRules[j]; 01006 objBest.objRules[j] = objBest.objRules[j + 1]; 01007 objBest.objRules[j + 1] = Temp; 01008 } 01009 } 01010 01011 if (j >= Popsize) 01012 return 0; 01013 01014 // Inserting new rule in j 01015 // Shift rules down 01016 for (l = iBestSize - 1; (l >= j); l--) 01017 { 01018 objBest.objRules[l + 1] = objBest.objRules[l]; 01019 } 01020 01021 objBest.objRules[j] = oNewRule->clone(); 01022 objBest.intRules++; 01023 01024 updateHeuOpPerformance(oNewRule->chrOrigin); 01025 01026 return 1; 01027 } 01028 01029 // ========================================================================== 01030 void GarpAlgorithm::saveBestRules(RuleSet * toRuleSet, RuleSet * fromRuleSet) 01031 { 01032 // deprecated!!! 01033 throw GarpException(1, "GarpAlgorithm::saveBestRules() method is deprecated"); 01034 01035 /* 01036 // put marks on both sets so it is easy to tell where each rule came from 01037 toRuleSet->setPad(' '); 01038 fromRuleSet->setPad('n'); 01039 01040 // concatenate the rulesets, replacing the similar rules and adding the different ones 01041 concatenateRuleSets(toRuleSet, fromRuleSet); 01042 01043 // sort the result set 01044 toRuleSet->sort(0); 01045 01046 // trim it to the max number of rules 01047 toRuleSet->trim(Popsize); 01048 01049 // count how many of the remaining rules came from the new set 01050 updateHeuOpPerformance(toRuleSet, 'n'); 01051 */ 01052 } 01053 01054 // ========================================================================== 01055 void GarpAlgorithm::concatenateRuleSets(RuleSet * toRuleSet, RuleSet * fromRuleSet) 01056 { 01057 // deprecated!!! 01058 throw GarpException(1, "GarpAlgorithm::concatenateRuleSets() method is deprecated"); 01059 01060 01061 Rule * fromRule; 01062 int i, j, nf; 01063 bool found; 01064 01065 //nt = toRuleSet->size(); 01066 nf = fromRuleSet->size(); 01067 01068 for (i = 0; i < nf; i++) 01069 { 01070 fromRule = fromRuleSet->get(i); 01071 01072 // Check if an identical or more general structure is already there 01073 for (j = 0, found = false; j < toRuleSet->size() && (!found); j++) 01074 { 01075 if (fromRule->type() == toRuleSet->get(j)->type()) 01076 found = fromRule->similar(toRuleSet->get(j)); 01077 } 01078 01079 j--; 01080 01081 if (found) 01082 toRuleSet->get(j)->copy(fromRule); 01083 else 01084 toRuleSet->add(fromRule->clone()); 01085 } 01086 } 01087 01088 // ========================================================================== 01089 void GarpAlgorithm::select() 01090 { 01091 double expected; // expected number of offspring 01092 double factor; // normalizer for expected value 01093 double ptr; // determines fractional selection 01094 double sum; // control for selection loop 01095 01096 int i, j, k, n, temp; 01097 01098 int Sample[MAX_RULES]; 01099 01100 n = objBest.size(); 01101 for (i = 0; i < MAX_RULES; i++) 01102 Sample[i] = i % n; 01103 01104 // normalizer for proportional selection probabilities 01105 if (Ave_current_perf - Worst) 01106 factor = Maxflag ? 1.0 / (Ave_current_perf - Worst) : 1.0/(Worst - Ave_current_perf); 01107 else 01108 factor=1.0; 01109 01110 01111 // Stochastic universal sampling algorithm by James E. Baker 01112 k=0; // index of next Selected structure 01113 ptr = GarpUtil::random(); // spin the wheel one time 01114 01115 for (sum = i = 0; i < objBest.size(); i++) 01116 { 01117 if (Maxflag) 01118 { 01119 if (objBest.get(i)->dblPerformance[0] > Worst) 01120 expected = (objBest.get(i)->dblPerformance[0] - Worst) * factor; 01121 else 01122 expected = 0.0; 01123 } 01124 01125 else 01126 { 01127 if (objBest.get(i)->dblPerformance[0] < Worst) 01128 expected = (Worst - objBest.get(i)->dblPerformance[0]) * factor; 01129 else 01130 expected = 0.0; 01131 } 01132 01133 for (sum += expected; (sum > ptr) && (k<=Popsize); ptr++) 01134 Sample[k++] = i; 01135 } 01136 01137 // randomly shuffle pointers to new structures 01138 for (i = 0; i < Popsize; i++) 01139 { 01140 j = GarpUtil::randint (i, Popsize - 1); 01141 temp = Sample[j]; 01142 Sample[j] = Sample[i]; 01143 Sample[i] = temp; 01144 } 01145 01146 // finally, form the new population 01147 // Gapsize giving the proportion contribution 01148 // to the new population from the objBest archive set 01149 Rule * oRuleBeingInserted; 01150 double dSize; 01151 01152 objNew.clear(); 01153 dSize = ((double) Popsize) * Gapsize; 01154 01155 for (i = 0; i < dSize; i++) 01156 { 01157 oRuleBeingInserted = objBest.objRules[Sample[i]]->clone(); 01158 oRuleBeingInserted->blnNeedsEvaluation = true; 01159 objNew.add(oRuleBeingInserted); 01160 } 01161 } 01162 01163 // ========================================================================== 01164 double GarpAlgorithm::converge() 01165 { 01166 if (Heuristic[1][0] != 0) 01167 Convergence = (Convergence + ((double)Conv)/Improvements)/2; 01168 else 01169 Convergence = 1.0; 01170 01171 return Convergence; 01172 } 01173 01174 // ========================================================================== 01175 void GarpAlgorithm::join() 01176 { 01177 return; 01178 01179 register int mom, dad; /* participants in the crossover */ 01180 register int xpoint1; /* first crossover point w.r.t. structure */ 01181 register int xpoint2; /* second crossover point w.r.t. structure */ 01182 register int i; /* loop control variable */ 01183 register BYTE temp; /* used for swapping alleles */ 01184 static double last; /* last element to undergo Crossover */ 01185 int diff; /* set if parents differ from offspring */ 01186 BYTE *kid1; /* pointers to the offspring */ 01187 BYTE *kid2; 01188 01189 last = Popsize * J_rate; 01190 01191 if (J_rate > 0.0) 01192 { 01193 for (mom = Popsize - 1; mom > Popsize-last ; mom -= 2) 01194 { 01195 /* diff wasn't beeing initialized in original Stockwell's code */ 01196 diff = 0; 01197 01198 dad = mom - 1; 01199 01200 /* kids start as identical copies of parents */ 01201 kid1 = objNew.get(mom)->Gene; 01202 kid2 = objNew.get(dad)->Gene; 01203 01204 /* choose two Crossover points */ 01205 xpoint1 = GarpUtil::randint(0, objNew.get(mom)->intLength); 01206 xpoint2 = GarpUtil::randint(0, objNew.get(mom)->intLength); 01207 01208 01209 /* perform crossover */ 01210 for (i=xpoint1 ; i % objNew.get(mom)->intLength == xpoint2; i++) 01211 { 01212 temp = kid1[i]; 01213 if (kid1[i] == MAX_BYTE) kid1[i] = kid2[i]; 01214 if (temp != MAX_BYTE) kid2[i] = temp; 01215 diff += (kid1[i] != kid2[i]); 01216 } 01217 01218 if (diff) /* kids differ from parents */ 01219 { 01220 /* set evaluation flags */ 01221 objNew.get(mom)->blnNeedsEvaluation = true; 01222 objNew.get(mom)->intGens = 0; 01223 objNew.get(mom)->chrOrigin = 'j'; 01224 } 01225 } 01226 } 01227 } 01228 01229 // ========================================================================== 01230 void GarpAlgorithm::mutate() 01231 { 01232 int i; 01233 int Temperature = MAX_MUTATION_TEMPERATURE; 01234 01235 if (Gen) Temperature = (int)(MAX_MUTATION_TEMPERATURE/(double)Gen); 01236 01237 for (i = 0; i < Popsize; i++) 01238 objNew.get(i)->mutate(Temperature); 01239 01240 return; 01241 } 01242 01243 // ========================================================================== 01244 void GarpAlgorithm::crossover() 01245 { 01246 return; 01247 01248 register int mom, dad; /* participants in the crossover */ 01249 register int xpoint1; /* first crossover point w.r.t. structure */ 01250 register int xpoint2; /* second crossover point w.r.t. structure */ 01251 register int i; /* loop control variable */ 01252 register char temp; /* used for swapping alleles */ 01253 static double last; /* last element to undergo Crossover */ 01254 int diff; /* set if parents differ from offspring */ 01255 BYTE *kid1; /* pointers to the offspring */ 01256 BYTE *kid2; 01257 01258 last = C_rate * Popsize; 01259 01260 if (C_rate > 0.0) 01261 { 01262 for (mom=0; mom < last ; mom += 2) 01263 { 01264 /* diff wasn't beeing initialized in original Stockwell's code */ 01265 diff = 0; 01266 01267 dad = mom + 1; 01268 01269 /* kids start as identical copies of parents */ 01270 kid1 = objNew.get(mom)->Gene; 01271 kid2 = objNew.get(dad)->Gene; 01272 01273 /* choose two Crossover points */ 01274 xpoint1 = GarpUtil::randint(0,objNew.get(mom)->intLength); 01275 xpoint2 = GarpUtil::randint(0,objNew.get(mom)->intLength); 01276 01277 01278 /* perform crossover */ 01279 for (i=xpoint1 ; i%(objNew.get(mom)->intLength)==xpoint2; i++) 01280 { 01281 temp = kid1[i]; 01282 kid1[i] = kid2[i]; 01283 kid2[i] = temp; 01284 diff += (kid1[i] != kid2[i]); 01285 } 01286 01287 if (diff) /* kids differ from parents */ 01288 { 01289 /* set evaluation flags */ 01290 objNew.get(mom)->blnNeedsEvaluation = true; 01291 objNew.get(dad)->blnNeedsEvaluation = true; 01292 objNew.get(mom)->intGens = 0; 01293 objNew.get(dad)->intGens = 0; 01294 objNew.get(mom)->chrOrigin = 'c'; 01295 objNew.get(dad)->chrOrigin = 'c'; 01296 } 01297 } 01298 } 01299 } 01300 01301 // ========================================================================== 01302 void GarpAlgorithm::getInitialModel(int intSize, EnvCellSet * objTrainSet) 01303 { 01304 Popsize = intSize; 01305 colonize(&objNew, objTrainSet, Popsize); 01306 } 01307 01308 // ========================================================================== 01309 01310 /****************************************************************/ 01311 /****************** configuration *******************************/ 01312 void 01313 GarpAlgorithm::_getConfiguration( ConfigurationPtr& config ) const 01314 { 01315 01316 // Only get the Model portion. Since the parameter portion is 01317 // handled by the base class AlgorithmImpl. 01318 if ( !done() ) 01319 return; 01320 01321 // Const hack. 01322 // To accomodate other's const-incorrectness. 01323 RuleSet* rs = const_cast<RuleSet*>( & objBest ); 01324 01325 ConfigurationPtr model_config ( new ConfigurationImpl("Garp") ); 01326 config->addSubsection( model_config ); 01327 01328 const int rule_count = rs->size(); 01329 01330 model_config->addNameValue( "ActiveGenes", iActiveGenes ); 01331 model_config->addNameValue( "RuleCount", rule_count ); 01332 01333 for ( int i=0; i<rule_count; i++ ) { 01334 01335 //Log::instance()->debug( "Gettting Rule %d\n",i ); 01336 01337 Rule *rule = rs->get(i); 01338 01339 ConfigurationPtr rule_config ( new ConfigurationImpl("Rule") ); 01340 model_config->addSubsection( rule_config ); 01341 01342 char type[2]; 01343 sprintf(type, "%c", rule->type() ); 01344 rule_config->addNameValue( "Type", type ); 01345 01346 rule_config->addNameValue( "Performance", rule->dblPerformance, 10 ); 01347 01348 rule_config->addNameValue( "Genes", rule->Gene, rule->intLength ); 01349 01350 } 01351 01352 01353 } 01354 01355 void 01356 GarpAlgorithm::_setConfiguration( const ConstConfigurationPtr& config ) 01357 { 01358 01359 ConstConfigurationPtr model_config = config->getSubsection( "Garp", false ); 01360 01361 if (!model_config) 01362 return; 01363 01364 Doneflag = true; 01365 // 01366 // Clear out the rule set. 01367 objBest.clear(); 01368 01369 iActiveGenes = model_config->getAttributeAsInt( "ActiveGenes", 0 ); 01370 01371 objBest.setDimension( iActiveGenes ); 01372 01373 Configuration::subsection_list::const_iterator ss; 01374 for( ss = model_config->getAllSubsections().begin(); 01375 ss != model_config->getAllSubsections().end(); 01376 ++ss ) { 01377 01378 const ConstConfigurationPtr& c(*ss); 01379 01380 std::string type = c->getAttribute( "Type" ); 01381 01382 double *perf; 01383 c->getAttributeAsDoubleArray( "Performance", &perf, 0 ); 01384 01385 int n_genes; 01386 unsigned char *p_genes; 01387 c->getAttributeAsByteArray( "Genes", &p_genes, &n_genes ); 01388 01389 Rule * newRule=0; 01390 switch( type [0] ) { 01391 01392 case 'a': 01393 newRule = new AtomicRule(); 01394 newRule->RestoreRule( perf, p_genes, n_genes, iGeneIndex ); 01395 break; 01396 01397 case 'd': 01398 newRule = new RangeRule(); 01399 newRule->RestoreRule( perf, p_genes, n_genes, iGeneIndex ); 01400 break; 01401 01402 case 'r': 01403 newRule = new LogitRule(); 01404 newRule->RestoreRule( perf, p_genes, n_genes, iGeneIndex ); 01405 break; 01406 01407 case '!': 01408 newRule = new NegatedRangeRule(); 01409 newRule->RestoreRule( perf, p_genes, n_genes, iGeneIndex ); 01410 break; 01411 01412 } 01413 01414 objBest.add(newRule); 01415 delete [] p_genes; 01416 delete [] perf; 01417 } 01418 01419 return; 01420 01421 }