openModeller
Version 1.4.0
|
00001 00034 #include "garp.hh" 00035 #include "rules_base.hh" 00036 #include "rules_range.hh" 00037 #include "rules_negrange.hh" 00038 #include "rules_logit.hh" 00039 #include "ruleset.hh" 00040 #include "bioclim_histogram.hh" 00041 #include "regression.hh" 00042 #include <openmodeller/Configuration.hh> 00043 #include <openmodeller/Random.hh> 00044 #include <openmodeller/Exceptions.hh> 00045 #include <openmodeller/ScaleNormalizer.hh> 00046 00047 #include <string> 00048 using std::string; 00049 00050 #define NUM_PARAM 4 00051 00052 /****************************************************************/ 00053 /*** Algorithm parameter metadata *******************************/ 00054 00055 AlgParamMetadata parameters[NUM_PARAM] = 00056 { 00057 // Metadata of the first parameter. 00058 { 00059 "MaxGenerations", // Id. 00060 "Max generations", // Name. 00061 Integer, // Type. 00062 00063 // Overview. 00064 "Maximum number of iterations run by the Genetic Algorithm.", 00065 00066 // Description. 00067 "Maximum number of iterations (generations) run by the Genetic\ 00068 Algorithm.", 00069 00070 1, // Not zero if the parameter has lower limit. 00071 1, // Parameter's lower limit. 00072 0, // Not zero if the parameter has upper limit. 00073 0, // Parameter's upper limit. 00074 "400" // Parameter's typical (default) value. 00075 }, 00076 00077 { 00078 "ConvergenceLimit", // Id. 00079 "Convergence limit", // Name. 00080 Real, // Type. 00081 00082 // Overview. 00083 "Defines the convergence value that makes the algorithm stop\ 00084 (before reaching MaxGenerations).", 00085 00086 // Description. 00087 "", 00088 00089 1, // Not zero if the parameter has lower limit. 00090 0.0, // Parameter's lower limit. 00091 1, // Not zero if the parameter has upper limit. 00092 1.0, // Parameter's upper limit. 00093 "0.01" // Parameter's typical (default) value. 00094 }, 00095 00096 { 00097 "PopulationSize", // Id. 00098 "Population size", // Name. 00099 Integer, // Type. 00100 00101 "Maximum number of rules to be kept in solution.", // Overview. 00102 "Maximum number of rules to be kept in solution.", // Description 00103 00104 1, // Not zero if the parameter has lower limit. 00105 1, // Parameter's lower limit. 00106 1, // Not zero if the parameter has upper limit. 00107 500, // Parameter's upper limit. 00108 "50" // Parameter's typical (default) value. 00109 }, 00110 00111 { 00112 "Resamples", // Id. 00113 "Resamples", // Name. 00114 Integer, // Type. 00115 00116 // Overview. 00117 "Number of points sampled (with replacement) used to test rules.", 00118 00119 // Description. 00120 "Number of points sampled (with replacement) used to test rules.", 00121 00122 1, // Not zero if the parameter has lower limit. 00123 1, // Parameter's lower limit. 00124 1, // Not zero if the parameter has upper limit. 00125 100000, // Parameter's upper limit. 00126 "2500" // Parameter's typical (default) value. 00127 } 00128 }; 00129 00130 00131 /*****************************************************************/ 00132 /*** Algorithm's general metadata ********************************/ 00133 00134 AlgMetadata metadata = { 00135 00136 "GARP", // Id. 00137 "GARP (single run) - new openModeller implementation", // Name. 00138 "3.3", // Version. 00139 00140 // Overview. 00141 "GARP is a genetic algorithm that creates ecological niche \ 00142 models for species. The models describe environmental conditions \ 00143 under which the species should be able to maintain populations. \ 00144 For input, GARP uses a set of point localities where the species \ 00145 is known to occur and a set of geographic layers representing \ 00146 the environmental parameters that might limit the species' \ 00147 capabilities to survive. Please refer to algorithm description for \ 00148 more information about the differences between this new GARP \ 00149 implementation and the Desktop GARP implementation.", 00150 00151 // Description. 00152 "GARP is a genetic algorithm that creates ecological niche \ 00153 models for species. The models describe environmental conditions \ 00154 under which the species should be able to maintain populations. \ 00155 For input, GARP uses a set of point localities where the species \ 00156 is known to occur and a set of geographic layers representing \ 00157 the environmental parameters that might limit the species' \ 00158 capabilities to survive. This implementation is a complete rewrite \ 00159 of the DesktopGarp code, and it also contains the following \ 00160 changes/improvements: (1) Gene values changed from integers (between 1 \ 00161 and 253) to floating point numbers (between -1.0 and 1.0). This avoids \ 00162 precision problems in environment values during projection (for example, \ 00163 if an environment variable has the value 2.56 in some raster cell and \ 00164 2.76 in another one, DesktopGarp rounds them off to 3). (2) Atomic rules \ 00165 were removed since they seem to have little significance compared to the \ 00166 other rules. (3) Heuristic operator parameters (percentage of mutation \ 00167 and crossover per iteration) are now static since they used to converge \ 00168 to fixed values during the very first iterations. This implementation \ 00169 simply keeps the converged values. (4) A bug was fixed in the procedure \ 00170 responsible for ordering the rules. When a rule was only replacing \ 00171 another, it was being included in the wrong position.", 00172 00173 // Author 00174 "Stockwell, D. R. B., modified by Ricardo Scachetti Pereira", 00175 00176 // Bibliography. 00177 "Stockwell, D. R. B. 1999. Genetic algorithms II. \ 00178 Pages 123-144 in A. H. Fielding, editor. \ 00179 Machine learning methods for ecological applications. \ 00180 Kluwer Academic Publishers, Boston.\ 00181 \n\ 00182 Stockwell, D. R. B., and D. P. Peters. 1999. \ 00183 The GARP modelling system: Problems and solutions to automated \ 00184 spatial prediction. International Journal of Geographic \ 00185 Information Systems 13:143-158.\ 00186 \n\ 00187 Stockwell, D. R. B., and I. R. Noble. 1992. \ 00188 Induction of sets of rules from animal distribution data: \ 00189 A robust and informative method of analysis. Mathematics and \ 00190 Computers in Simulation 33:385-390.", 00191 00192 "Ricardo Scachetti Pereira", // Code author. 00193 "ricardo [at] tdwg . org", // Code author's contact. 00194 00195 0, // Does not accept categorical data. 00196 1, // Does not need (pseudo)absence points. 00197 00198 NUM_PARAM, // Algorithm's parameters. 00199 parameters 00200 }; 00201 00202 00203 00204 /****************************************************************/ 00205 /****************** Algorithm's factory function ****************/ 00206 #ifndef DONT_EXPORT_GARP_FACTORY 00207 OM_ALG_DLL_EXPORT 00208 AlgorithmImpl * 00209 algorithmFactory() 00210 { 00211 return new Garp(); 00212 } 00213 00214 OM_ALG_DLL_EXPORT 00215 AlgMetadata const * 00216 algorithmMetadata() 00217 { 00218 return &metadata; 00219 } 00220 #endif 00221 00222 00223 /****************************************************************/ 00224 /****************** Garp class **********************************/ 00225 00226 const PerfIndex defaultPerfIndex = PerfSig; 00227 00228 /****************************************************************/ 00229 /****************** Garp constructor ****************************/ 00230 00231 Garp::Garp() 00232 : AlgorithmImpl(& metadata) 00233 { 00234 // fill in default values for parameters 00235 _popsize = 0; 00236 _resamples = 0; 00237 _max_gen = 0; 00238 _conv_limit = 0.0; 00239 00240 _mortality = 0.9; 00241 _gapsize = 0.1; 00242 _acc_limit = 0.0; 00243 00244 _crossover_rate = 0.1; 00245 _mutation_rate = 0.6; 00246 00247 _significance = 2.70; 00248 00249 _maxProgress = 0.0; 00250 00251 // reset private attributes 00252 _fittest = _offspring = NULL; 00253 00254 _gen = 0; 00255 _convergence = 1.0; 00256 _improvements = 0; 00257 00258 int i; 00259 for (i = 0; i < 5; i++) 00260 { 00261 _curr_heur_count[i] = 0; 00262 _prev_heur_count[i] = 0; 00263 } 00264 00265 _normalizerPtr = new ScaleNormalizer( -1.0, +1.0, true ); 00266 } 00267 00268 /****************************************************************/ 00269 /****************** Garp destructor *****************************/ 00270 Garp::~Garp() 00271 { 00272 // debug 00273 if ( _fittest ) 00274 { 00275 //Log::instance()->debug( "Resulting rules:\n"); 00276 //_fittest->log(); 00277 } 00278 00279 if (_offspring) 00280 delete _offspring; 00281 00282 if (_fittest) 00283 delete _fittest; 00284 } 00285 00286 00287 // **************************************************************** 00288 // ************* initialize *************************************** 00289 00290 int Garp::initialize() 00291 { 00292 if (!getParameter("MaxGenerations", &_max_gen)) { 00293 Log::instance()->error("Parameter MaxGenerations not set properly."); 00294 return 0; 00295 } 00296 00297 if (!getParameter("ConvergenceLimit", &_conv_limit)) { 00298 Log::instance()->error("Parameter ConvergenceLimit not set properly."); 00299 return 0; 00300 } 00301 00302 if (!getParameter("PopulationSize", &_popsize)) { 00303 Log::instance()->error("Parameter PopulationSize not set properly."); 00304 return 0; 00305 } 00306 00307 if (!getParameter("Resamples", &_resamples)) { 00308 Log::instance()->error("Parameter Resamples not set properly."); 00309 return 0; 00310 } 00311 00312 //Log::instance()->debug("MaxGenerations set to: %d\n", _max_gen); 00313 //Log::instance()->debug("ConvergenceLimit set to: %.4f\n", _conv_limit); 00314 //Log::instance()->debug("PopulationSize set to: %d\n", _popsize); 00315 //Log::instance()->debug("Resamples set to: %d\n", _resamples); 00316 00317 _offspring = new GarpRuleSet(2 * _popsize); 00318 _fittest = new GarpRuleSet(2 * _popsize); 00319 00320 cacheSamples(_samp, _cachedOccs, _resamples); 00321 _bioclimHistogram.initialize(_cachedOccs); 00322 _regression.calculateParameters(_cachedOccs); 00323 00324 colonize(_offspring, _popsize); 00325 00326 return 1; 00327 } 00328 00329 /****************************************************************/ 00330 /****************************************************************/ 00331 void Garp::cacheSamples(const SamplerPtr& sampler, 00332 OccurrencesPtr& cachedOccs, 00333 int resamples) 00334 { 00335 OccurrencesImpl * occs = new OccurrencesImpl( "", 00336 GeoTransform::getDefaultCS() ); 00337 occs->reserve( resamples ); 00338 cachedOccs = ReferenceCountedPointer<OccurrencesImpl>( occs ); 00339 00340 for (int i = 0; i < resamples; ++i) 00341 { 00342 OccurrencePtr oc = sampler->getOneSample(); 00343 cachedOccs->insert(oc); 00344 } 00345 } 00346 00347 00348 /****************************************************************/ 00349 /****************** iterate *************************************/ 00350 00351 int Garp::iterate() 00352 { 00353 double perfBest, perfWorst, perfAvg; 00354 00355 if (done()) 00356 return 1; 00357 00358 _gen++; 00359 00360 evaluate(_offspring); 00361 keepFittest(_offspring, _fittest, defaultPerfIndex); 00362 _fittest->trim(_popsize); 00363 00364 //_fittest->gatherRuleSetStats(_gen); 00365 //_offspring->gatherRuleSetStats(-_gen); 00366 00367 00368 _fittest->performanceSummary(defaultPerfIndex, 00369 &perfBest, &perfWorst, &perfAvg); 00370 /* 00371 // log info about current iteration 00372 Log::instance()->debug( "%4d] ", _gen ); 00373 Log::instance()->debug( "[%2d] conv=%+7.4f | perfs=%+8.3f, %+8.3f, %+8.3f\n", _fittest->numRules(), 00374 _convergence, perfBest, perfWorst, perfAvg ); 00375 */ 00376 00377 if (done()) 00378 { 00379 // finalize processing of model 00380 // by filtering out rules that have low performance 00381 _fittest->filter(defaultPerfIndex, _significance); 00382 return 1; 00383 } 00384 00385 // algorithm is not done yet 00386 // select fittest individuals 00387 select(_fittest, _offspring, _gapsize); 00388 00389 // create new offspring 00390 colonize(_offspring, _popsize); 00391 _offspring->trim(_popsize); 00392 mutate(_offspring); 00393 crossover(_offspring); 00394 00395 return 1; 00396 } 00397 00398 /****************************************************************/ 00399 /****************** getProgress *********************************/ 00400 00401 float Garp::getProgress() const 00402 { 00403 if (done()) 00404 { return 1.0; } 00405 else 00406 { 00407 float byIterations = ( _gen / (float) _max_gen ); 00408 float byConvergence = (float)( _conv_limit / _convergence ); 00409 float progress = (byIterations > byConvergence) ? byIterations : byConvergence; 00410 if (progress > _maxProgress) 00411 { _maxProgress = progress; } 00412 return _maxProgress; 00413 } 00414 } 00415 00416 /****************************************************************/ 00417 /****************** done ****************************************/ 00418 00419 int Garp::done() const 00420 { 00421 return ( (_gen >= _max_gen) || (_convergence < _conv_limit) ); 00422 } 00423 00424 /****************************************************************/ 00425 /****************** getValue ************************************/ 00426 00427 Scalar Garp::getValue( const Sample& x ) const 00428 { 00429 return _fittest->getValue(x); 00430 } 00431 00432 /****************************************************************/ 00433 /****************** getConvergence ******************************/ 00434 00435 int Garp::getConvergence( Scalar * const val ) const 00436 { 00437 *val = _convergence; 00438 return 0; 00439 } 00440 00441 /****************************************************************/ 00442 /****************** configuration *******************************/ 00443 void 00444 Garp::_getConfiguration( ConfigurationPtr& config ) const 00445 { 00446 if ( !done() ) 00447 return; 00448 00449 ConfigurationPtr model_config ( new ConfigurationImpl("Garp") ); 00450 config->addSubsection( model_config ); 00451 00452 model_config->addNameValue( "Generations", _gen ); 00453 model_config->addNameValue( "AccuracyLimit", _acc_limit ); 00454 model_config->addNameValue( "Mortality", _mortality ); 00455 model_config->addNameValue( "Significance", _significance ); 00456 model_config->addNameValue( "FinalCrossoverRate", _crossover_rate ); 00457 model_config->addNameValue( "FinalMutationRate", _mutation_rate ); 00458 model_config->addNameValue( "FinalGapSize", _gapsize ); 00459 00460 if ( _fittest ) { 00461 00462 int nrules = _fittest->numRules(); 00463 00464 ConfigurationPtr rules_config( new ConfigurationImpl("FittestRules") ); 00465 model_config->addSubsection( rules_config ); 00466 00467 rules_config->addNameValue( "Count", nrules ); 00468 00469 for( int i=0; i<nrules; i++ ) { 00470 00471 GarpRule *rule = _fittest->get(i); 00472 char type[16]; 00473 sprintf(type, "%c", rule->type() ); 00474 00475 ConfigurationPtr rule_config( new ConfigurationImpl("Rule") ); 00476 rules_config->addSubsection( rule_config ); 00477 00478 rule_config->addNameValue( "Type", type ); 00479 rule_config->addNameValue( "Prediction", rule->getPrediction() ); 00480 rule_config->addNameValue( "Chromosome1", rule->getChrom1()); 00481 rule_config->addNameValue( "Chromosome2", rule->getChrom2()); 00482 rule_config->addNameValue( "Performance", rule->getPerformanceArray(), 10 ); 00483 } 00484 } 00485 00486 } 00487 00488 void 00489 Garp::_setConfiguration( const ConstConfigurationPtr& config ) 00490 { 00491 ConstConfigurationPtr model_config = config->getSubsection( "Garp", false ); 00492 00493 if (!model_config) { 00494 00495 return; 00496 } 00497 00498 _gen = model_config->getAttributeAsInt( "Generations", 0 ); 00499 _acc_limit = model_config->getAttributeAsDouble( "AccuracyLimit", 0.0 ); 00500 _mortality = model_config->getAttributeAsDouble( "Mortality", 0.0 ); 00501 _significance = model_config->getAttributeAsDouble( "Significance", 0.0 ); 00502 _crossover_rate = model_config->getAttributeAsDouble( "FinalCrossoverRate", 0.0 ); 00503 _mutation_rate = model_config->getAttributeAsDouble( "FinalMutationRate", 0.0 ); 00504 _gapsize = model_config->getAttributeAsDouble( "FinalGapSize", 0.0 ); 00505 00506 // Need to read at least this parameter 00507 if ( ! getParameter("PopulationSize", &_popsize) ) { 00508 00509 Log::instance()->error("Could not read parameter PopulationSize from serialized model."); 00510 return; 00511 } 00512 00513 _offspring = new GarpRuleSet( 2 * _popsize ); 00514 _fittest = new GarpRuleSet( 2 * _popsize ); 00515 00516 /* 00517 * This code is commented out for now. Need to figure out how 00518 * to get the algorithm primed with its custom sampler after 00519 * it's deserialized. 00520 */ 00521 //_bioclimHistogram.initialize( _cachedOccs ); 00522 00523 ConstConfigurationPtr rules_config = model_config->getSubsection( "FittestRules" ); 00524 00525 //next line commented out since the var is not used after being declared 00526 //int nrules = rules_config->getAttributeAsInt( "Count", 0 ); 00527 00528 Configuration::subsection_list::const_iterator ss; 00529 for( ss = rules_config->getAllSubsections().begin(); 00530 ss != rules_config->getAllSubsections().end(); 00531 ++ss ) { 00532 00533 const ConstConfigurationPtr& c(*ss); 00534 GarpRule * rule = NULL; 00535 00536 string type = c->getAttribute( "Type" ); 00537 00538 Scalar pred = c->getAttributeAsDouble( "Prediction", 0.0 ); 00539 00540 Sample p_chrom1 = c->getAttributeAsSample( "Chromosome1" ); 00541 00542 Sample p_chrom2 = c->getAttributeAsSample( "Chromosome2" ); 00543 00544 Scalar *p_perf; 00545 int n_perf; 00546 c->getAttributeAsDoubleArray( "Performance", &p_perf, &n_perf ); 00547 00548 switch( type[0] ) { 00549 case 'd': 00550 rule = new RangeRule( pred, p_chrom1.size(), p_chrom1, p_chrom2, p_perf ); 00551 break; 00552 00553 case 'r': 00554 rule = new LogitRule( pred, p_chrom1.size(), p_chrom1, p_chrom2, p_perf ); 00555 break; 00556 00557 case '!': 00558 rule = new NegatedRangeRule( pred, p_chrom1.size(), p_chrom1, p_chrom2, p_perf ); 00559 break; 00560 00561 } 00562 00563 delete [] p_perf; 00564 00565 _fittest->add(rule); 00566 00567 } 00568 00569 } 00570 00571 /****************************************************************/ 00572 /***************** GARP Algorithm private methods ***************/ 00573 /****************************************************************/ 00574 00575 00576 /****************************************************************/ 00577 /***************** keepFittest **********************************/ 00578 00579 void Garp::keepFittest(GarpRuleSet * source, GarpRuleSet * target, 00580 PerfIndex perfIndex) 00581 { 00582 int i, n, converged, similarIndex; 00583 GarpRule * candidateRule, * similarRule; 00584 00585 converged = 0; 00586 00587 // step through source rule-set trying to insert rules into target 00588 n = source->numRules(); 00589 for (i = 0; i < n; i++) 00590 { 00591 candidateRule = source->get(i); 00592 similarIndex = target->findSimilar(candidateRule); 00593 //if ((similarIndex < -1) || (similarIndex >= target->numRules())) 00594 // Log::instance()->error("Index out of bounds (#8). Limits are (-1, %d), index is %d\n", target->numRules(), similarIndex); 00595 00596 if (similarIndex >= 0) 00597 { 00598 converged++; 00599 00600 // similar rule found replace it if better 00601 similarRule = target->get(similarIndex); 00602 if (candidateRule->getPerformance(perfIndex) > 00603 similarRule->getPerformance(perfIndex)) 00604 { 00605 // first create a clone of the rule, then replace the old one 00606 candidateRule = candidateRule->clone(); 00607 //target->replace(similarIndex, candidateRule); 00608 target->remove(similarIndex); 00609 target->insert(perfIndex, candidateRule); 00610 } 00611 } 00612 else 00613 { 00614 // no similar rule found: try to insert it into existing set 00615 // first create a clone of the rule, then insert it into 00616 // the target rs 00617 candidateRule = candidateRule->clone(); 00618 target->insert(perfIndex, candidateRule); 00619 } 00620 } 00621 00622 // update convergence value 00623 _improvements += converged; 00624 if (_improvements) 00625 { _convergence = ( _convergence + ( (double) converged ) / _improvements ) / 2.0; } 00626 else 00627 { _convergence = 1.0; } 00628 00629 /* 00630 printf("Convergence: %+7.4f at generation %5d (%3d; %6d; %+7.4f) similar=%d\n", _convergence, 00631 _gen, converged, _improvements, ( (double) converged ) / _improvements, similarIndex); 00632 */ 00633 00634 // TODO: update heuristic rates based on who entered the target rule-set 00635 } 00636 00637 /****************************************************************/ 00638 /***************** evaluate *************************************/ 00639 00640 void Garp::evaluate(GarpRuleSet * ruleset) 00641 { 00642 int i, n; 00643 00644 n = ruleset->numRules(); 00645 for (i = 0; i < n; i++) 00646 { 00647 ruleset->get(i)->evaluate(_cachedOccs); 00648 } 00649 00650 return; 00651 } 00652 00653 /****************************************************************/ 00654 /***************** colonize *************************************/ 00655 00656 void Garp::colonize(GarpRuleSet * ruleset, int numRules) 00657 { 00658 int i, p, dim; 00659 GarpRule * rule = 0; 00660 Random rnd; 00661 00662 dim = _samp->numIndependent(); 00663 00664 for (i = ruleset->numRules(); i < numRules; i++) 00665 { 00666 // pick the next rule to be generated 00667 p = rnd(3); 00668 00669 switch (p) 00670 { 00671 case 0: 00672 rule = new RangeRule(dim); 00673 rule->setPrediction(1.0); 00674 ((RangeRule *) rule)->initialize(_bioclimHistogram); 00675 break; 00676 00677 case 1: 00678 rule = new NegatedRangeRule(dim); 00679 rule->setPrediction(0.0); 00680 ((NegatedRangeRule *) rule)->initialize(_bioclimHistogram); 00681 break; 00682 00683 case 2: 00684 rule = new LogitRule(dim); 00685 Scalar pred = (rnd.get(0.0, 1.0) > 0.5) ? 1.0 : 0.0; 00686 rule->setPrediction(pred); 00687 ((LogitRule *) rule)->initialize(_regression); 00688 break; 00689 } 00690 00691 //Log::instance()->debug("[%c] ", rule->type()); 00692 ruleset->add(rule); 00693 } 00694 } 00695 00696 /****************************************************************/ 00697 /***************** select ***************************************/ 00698 00699 void Garp::select(GarpRuleSet * source, GarpRuleSet * target, 00700 double gapsize) 00701 { 00702 Random rnd; 00703 int * sample; 00704 int i, j, k, n, temp; 00705 double perfBest, perfWorst, perfAvg; 00706 double sum, ptr, factor, expected, rulePerf, size; 00707 perfBest = perfWorst = perfAvg = 0.0; 00708 GarpRule * pRuleBeingInserted; 00709 00710 source->performanceSummary(defaultPerfIndex, &perfBest, &perfWorst, &perfAvg); 00711 00712 //Log::instance()->debug( "Performances: %f %f %f.\n", perfBest, perfWorst, perfAvg ); 00713 00714 // normalizer for proportional selection probabilities 00715 if (perfAvg - perfWorst) 00716 factor = 1.0 / (perfAvg - perfWorst); 00717 else 00718 factor = 1.0; 00719 00720 // Stochastic universal sampling algorithm by James E. Baker 00721 k = 0; 00722 n = source->numRules(); 00723 sample = new int[_popsize + 1]; 00724 for (i = 0; i < _popsize; i++) 00725 sample[i] = i % n; 00726 00727 ptr = rnd.get(1.0); 00728 sum = 0.0; 00729 for (i = 0; i < n; i++) { 00730 rulePerf = source->get(i)->getPerformance(defaultPerfIndex); 00731 expected = (rulePerf - perfWorst) * factor; 00732 for (sum += expected; (sum > ptr) && (k <= _popsize); ptr++) { 00733 if ((k < 0) || (k > _popsize)) { 00734 Log::instance()->error("Index out of bounds (#6). Limits are (0, %d), index is %d\n", _popsize, k); 00735 throw AlgorithmException("Index out of bounds"); 00736 } 00737 sample[k++] = i; 00738 } 00739 } 00740 00741 /* 00742 FILE * f = stdout; 00743 fprintf(f, "Generation: %4d\n", _gen); 00744 for (i = 0; i < _popsize; i++) 00745 fprintf(f, "%+9.4f %3d\n", source->get(sample[i])->getPerformance(defaultPerfIndex), sample[i]); 00746 */ 00747 00748 // randomly shuffle pointers to new structures 00749 for (i = 0; i < _popsize; i++) 00750 { 00751 j = rnd.get (i , _popsize - 1); 00752 temp = sample[j]; 00753 sample[j] = sample[i]; 00754 sample[i] = temp; 00755 } 00756 00757 // finally, form the new population 00758 // Gapsize giving the proportion contribution 00759 // to the new population from the objBest archive set 00760 target->clear(); 00761 size = ((double) _popsize) * gapsize; 00762 00763 for (i = 0; i < size; i++) 00764 { 00765 pRuleBeingInserted = source->get(sample[i])->clone(); 00766 pRuleBeingInserted->forceEvaluation(); 00767 if (!target->add(pRuleBeingInserted)) { 00768 // target rs is full 00769 std::string error = "Garp::reproduce(): Target rule set is full"; 00770 Log::instance()->error(error.c_str()); 00771 throw AlgorithmException(error.c_str()); 00772 } 00773 } 00774 delete[] sample; 00775 } 00776 00777 /****************************************************************/ 00778 /***************** mutate ***************************************/ 00779 00780 void Garp::mutate(GarpRuleSet * ruleset) 00781 { 00782 int i, n; 00783 00784 double temperature = 2.0 / (double) _gen; 00785 n = ruleset->numRules(); 00786 for (i = 0; i < n; i++) 00787 ruleset->get(i)->mutate(temperature); 00788 } 00789 00790 /****************************************************************/ 00791 /***************** crossover ************************************/ 00792 00793 void Garp::crossover(GarpRuleSet * ruleset) 00794 { 00795 Random rnd; 00796 int nrules, genes, xcount, last, mom, dad, xpt1, xpt2; 00797 00798 genes = _samp->numIndependent(); 00799 nrules = ruleset->numRules(); 00800 last = (int) (_crossover_rate * (double) nrules); 00801 00802 for (xcount = 0; xcount < last; xcount += 2) 00803 { 00804 mom = rnd.get(nrules); 00805 dad = rnd.get(nrules); 00806 if (dad == mom) 00807 dad = (dad + 1) % nrules; 00808 00809 xpt1 = rnd.get(genes); 00810 xpt2 = rnd.get(genes); 00811 00812 ruleset->get(mom)->crossover(ruleset->get(dad), xpt1, xpt2); 00813 } 00814 } 00815 00816 00817 // ***************** 00818 void Garp::deleteTempDataMembers() 00819 { 00820 if (_offspring) 00821 delete _offspring; 00822 _offspring = NULL; 00823 } 00824 00825 /****************************************************************/ 00826 /**** This is a debug function that checks if a rule set is 00827 **** correctly sorted. If not it dumps the performance values 00828 **** for that rule set. 00829 **** It was used to debug Garp::keepFittest() (replace call bug) 00830 **** TODO: move this code to the test harness when we have one */ 00831 00832 void printPerfs(char * msg, int index, GarpRuleSet * ruleset) 00833 { 00834 for (int i = 1; i < ruleset->numRules(); i++) 00835 { 00836 if (ruleset->get(i - 1)->getPerformance((PerfIndex)8) < 00837 ruleset->get(i)->getPerformance((PerfIndex)8)) 00838 { 00839 printf("\nError: rule set out of sort order (Index: %d)\n", index); 00840 for (int i = 0; i < ruleset->numRules(); i++) 00841 { 00842 printf("[%2d]=%6.3f ", i, ruleset->get(i)->getPerformance((PerfIndex)8) ); 00843 if ((i + 1) % 5 == 0) 00844 printf("\n"); 00845 } 00846 } 00847 } 00848 } 00849