Merge pull request #36 from heal-research/local-search-refactor

Refactor local search

CMakeLists.txt

@@ -66,6 +66,7 @@ add_library(
  source/operators/generator/brood.cpp
  source/operators/generator/os.cpp
  source/operators/generator/poly.cpp
+ source/operators/local_search.cpp
  source/operators/mutation.cpp
  source/operators/non_dominated_sorter/best_order_sort.cpp
  source/operators/non_dominated_sorter/deductive_sort.cpp

cli/source/operator_factory.cpp

@@ -4,14 +4,14 @@
 #include "operator_factory.hpp"
 #include <stdexcept> // for runtime_error
 #include <fmt/format.h> // for format
-#include <tuple>
 #include <scn/scan.h>
 #include "operon/interpreter/dispatch_table.hpp"
 #include "operon/operators/creator.hpp" // for CreatorBase, BalancedTreeC...
 #include "operon/operators/evaluator.hpp" // for Evaluator, EvaluatorBase
 #include "operon/operators/generator.hpp" // for OffspringGeneratorBase
 #include "operon/operators/reinserter.hpp" // for OffspringGeneratorBase
 #include "operon/operators/selector.hpp"
+#include "operon/operators/local_search.hpp"
 #include "operon/optimizer/optimizer.hpp"
 
 #include <cxxopts.hpp>
@@ -123,8 +123,10 @@ auto ParseEvaluator(std::string const& str, Problem& problem, DefaultDispatch& d
  evaluator = std::make_unique<Operon::Evaluator<T>>(problem, dtable, Operon::RMSE{}, scale);
  } else if (str == "mae") {
  evaluator = std::make_unique<Operon::Evaluator<T>>(problem, dtable, Operon::MAE{}, scale);
- } else if (str == "mdl") {
- evaluator = std::make_unique<Operon::MinimumDescriptionLengthEvaluator<T>>(problem, dtable);
+ } else if (str == "mdl_gauss") {
+ evaluator = std::make_unique<Operon::MinimumDescriptionLengthEvaluator<T, GaussianLikelihood<Operon::Scalar>>>(problem, dtable);
+ } else if (str == "mdl_poisson") {
+ evaluator = std::make_unique<Operon::MinimumDescriptionLengthEvaluator<T, PoissonLikelihood<Operon::Scalar>>>(problem, dtable);
  } else if (str == "gauss") {
  evaluator = std::make_unique<Operon::GaussianLikelihoodEvaluator<T>>(problem, dtable);
  } else {
@@ -133,13 +135,13 @@ auto ParseEvaluator(std::string const& str, Problem& problem, DefaultDispatch& d
  return evaluator;
 }
 
-auto ParseGenerator(std::string const& str, EvaluatorBase& eval, CrossoverBase& cx, MutatorBase& mut, SelectorBase& femSel, SelectorBase& maleSel) -> std::unique_ptr<OffspringGeneratorBase>
+auto ParseGenerator(std::string const& str, EvaluatorBase& eval, CrossoverBase& cx, MutatorBase& mut, SelectorBase& femSel, SelectorBase& maleSel, CoefficientOptimizer const* coeffOptimizer = nullptr) -> std::unique_ptr<OffspringGeneratorBase>
 {
  std::unique_ptr<OffspringGeneratorBase> generator;
  auto tok = Split(str, ':');
  auto name = tok[0];
  if (name == "basic") {
- generator = std::make_unique<BasicOffspringGenerator>(eval, cx, mut, femSel, maleSel);
+ generator = std::make_unique<BasicOffspringGenerator>(eval, cx, mut, femSel, maleSel, coeffOptimizer);
  } else if (name == "os") {
  size_t maxSelectionPressure{100};
  double comparisonFactor{0};
@@ -149,16 +151,16 @@ auto ParseGenerator(std::string const& str, EvaluatorBase& eval, CrossoverBase&
  if (tok.size() > 2) {
  comparisonFactor = scn::scan<double>(tok[2], "{}")->value();
  }
- generator = std::make_unique<OffspringSelectionGenerator>(eval, cx, mut, femSel, maleSel);
+ generator = std::make_unique<OffspringSelectionGenerator>(eval, cx, mut, femSel, maleSel, coeffOptimizer);
  dynamic_cast<OffspringSelectionGenerator*>(generator.get())->MaxSelectionPressure(maxSelectionPressure);
  dynamic_cast<OffspringSelectionGenerator*>(generator.get())->ComparisonFactor(comparisonFactor);
  } else if (name == "brood") {
- generator = std::make_unique<BroodOffspringGenerator>(eval, cx, mut, femSel, maleSel);
+ generator = std::make_unique<BroodOffspringGenerator>(eval, cx, mut, femSel, maleSel, coeffOptimizer);
  size_t broodSize{BroodOffspringGenerator::DefaultBroodSize};
  if (tok.size() > 1) { broodSize = scn::scan<size_t>(tok[1], "{}")->value(); }
  dynamic_cast<BroodOffspringGenerator*>(generator.get())->BroodSize(broodSize);
  } else if (name == "poly") {
- generator = std::make_unique<PolygenicOffspringGenerator>(eval, cx, mut, femSel, maleSel);
+ generator = std::make_unique<PolygenicOffspringGenerator>(eval, cx, mut, femSel, maleSel, coeffOptimizer);
  size_t polygenicSize{PolygenicOffspringGenerator::DefaultBroodSize};
  if (tok.size() > 1) { polygenicSize = scn::scan<size_t>(tok[1], "{}")->value(); }
  dynamic_cast<PolygenicOffspringGenerator*>(generator.get())->PolygenicSize(polygenicSize);
@@ -168,7 +170,7 @@ auto ParseGenerator(std::string const& str, EvaluatorBase& eval, CrossoverBase&
  return generator;
 }
 
-auto ParseOptimizer(std::string const& /*str*/, Problem const& /*problem*/, DefaultDispatch const& /*dtable*/) -> std::unique_ptr<OptimizerBase<DefaultDispatch>> {
+auto ParseOptimizer(std::string const& /*str*/, Problem const& /*problem*/, DefaultDispatch const& /*dtable*/) -> std::unique_ptr<OptimizerBase> {
  throw std::runtime_error("not implemented");
 }
 

cli/source/operator_factory.hpp

@@ -26,6 +26,7 @@ namespace Operon { class SelectorBase; }
 namespace Operon { struct CreatorBase; }
 namespace Operon { struct CrossoverBase; }
 namespace Operon { struct ErrorMetric; }
+namespace Operon { class CoefficientOptimizer; }
 namespace Operon { struct MutatorBase; }
 namespace Operon { struct Variable; }
 
@@ -41,9 +42,9 @@ auto ParseEvaluator(std::string const& str, Problem& problem, DefaultDispatch& d
 
 auto ParseErrorMetric(std::string const& str) -> std::tuple<std::unique_ptr<Operon::ErrorMetric>, bool>;
 
-auto ParseGenerator(std::string const& str, EvaluatorBase& eval, CrossoverBase& cx, MutatorBase& mut, SelectorBase& femSel, SelectorBase& maleSel) -> std::unique_ptr<OffspringGeneratorBase>;
+auto ParseGenerator(std::string const& str, EvaluatorBase& eval, CrossoverBase& cx, MutatorBase& mut, SelectorBase& femSel, SelectorBase& maleSel, CoefficientOptimizer const* cOpt) -> std::unique_ptr<OffspringGeneratorBase>;
 
-auto ParseOptimizer(std::string const& str, Problem const& problem, DefaultDispatch const& dtable) -> std::unique_ptr<OptimizerBase<DefaultDispatch>>;
+auto ParseOptimizer(std::string const& str, Problem const& problem, DefaultDispatch const& dtable) -> std::unique_ptr<OptimizerBase>;
 
 } // namespace Operon
 

cli/source/operon_gp.cpp

@@ -222,7 +222,8 @@ auto main(int argc, char** argv) -> int
 
  auto optimizer = std::make_unique<Operon::LevenbergMarquardtOptimizer<decltype(dtable), Operon::OptimizerType::Eigen>>(dtable, problem);
  optimizer->SetIterations(config.Iterations);
- dynamic_cast<Operon::Evaluator<Operon::DefaultDispatch>*>(evaluator.get())->SetOptimizer(optimizer.get());
+
+ Operon::CoefficientOptimizer cOpt{*optimizer, config.LamarckianProbability};
 
  EXPECT(problem.TrainingRange().Size() > 0);
 
@@ -231,7 +232,7 @@ auto main(int argc, char** argv) -> int
  auto femaleSelector = Operon::ParseSelector(result["female-selector"].as<std::string>(), comp);
  auto maleSelector = Operon::ParseSelector(result["male-selector"].as<std::string>(), comp);
 
- auto generator = Operon::ParseGenerator(result["offspring-generator"].as<std::string>(), *evaluator, crossover, mutator, *femaleSelector, *maleSelector);
+ auto generator = Operon::ParseGenerator(result["offspring-generator"].as<std::string>(), *evaluator, crossover, mutator, *femaleSelector, *maleSelector, &cOpt);
  auto reinserter = Operon::ParseReinserter(result["reinserter"].as<std::string>(), comp);
 
  Operon::RandomGenerator random(config.Seed);

cli/source/operon_nsgp.cpp

@@ -51,6 +51,8 @@ auto main(int argc, char** argv) -> int
  config.Iterations = result["iterations"].as<size_t>();
  config.CrossoverProbability = result["crossover-probability"].as<Operon::Scalar>();
  config.MutationProbability = result["mutation-probability"].as<Operon::Scalar>();
+ config.LocalSearchProbability = result["local-search-probability"].as<Operon::Scalar>();
+ config.LamarckianProbability = result["lamarckian-probability"].as<Operon::Scalar>();
  config.TimeLimit = result["timelimit"].as<size_t>();
  config.Seed = std::random_device {}();
 
@@ -229,7 +231,6 @@ auto main(int argc, char** argv) -> int
 
  auto optimizer = std::make_unique<Operon::LevenbergMarquardtOptimizer<decltype(dtable), Operon::OptimizerType::Eigen>>(dtable, problem);
  optimizer->SetIterations(config.Iterations);
- dynamic_cast<Operon::Evaluator<Operon::DefaultDispatch>*>(errorEvaluator.get())->SetOptimizer(optimizer.get());
  Operon::LengthEvaluator lengthEvaluator(problem, maxLength);
 
  Operon::MultiEvaluator evaluator(problem);
@@ -243,8 +244,9 @@ auto main(int argc, char** argv) -> int
 
  auto femaleSelector = Operon::ParseSelector(result["female-selector"].as<std::string>(), comp);
  auto maleSelector = Operon::ParseSelector(result["male-selector"].as<std::string>(), comp);
+ Operon::CoefficientOptimizer cOpt{*optimizer, config.LamarckianProbability};
 
- auto generator = Operon::ParseGenerator(result["offspring-generator"].as<std::string>(), evaluator, crossover, mutator, *femaleSelector, *maleSelector);
+ auto generator = Operon::ParseGenerator(result["offspring-generator"].as<std::string>(), evaluator, crossover, mutator, *femaleSelector, *maleSelector, &cOpt);
  auto reinserter = Operon::ParseReinserter(result["reinserter"].as<std::string>(), comp);
 
  Operon::RandomGenerator random(config.Seed);
@@ -373,6 +375,8 @@ auto main(int argc, char** argv) -> int
 
  using T = std::tuple<std::string, double, std::string>;
  auto const* format = ":>#8.3g"; // see https://fmt.dev/latest/syntax.html
+
+ auto [resEval, jacEval, callCount, cfTime ] = evaluator.Stats();
  std::array stats {
  T{ "iteration", gp.Generation(), ":>" },
  T{ "r2_tr", r2Train, format },
@@ -383,10 +387,10 @@ auto main(int argc, char** argv) -> int
  T{ "nmse_te", nmseTest, format },
  T{ "avg_fit", avgQuality, format },
  T{ "avg_len", avgLength, format },
- T{ "eval_cnt", evaluator.CallCount , ":>" },
- T{ "res_eval", evaluator.ResidualEvaluations, ":>" },
- T{ "jac_eval", evaluator.JacobianEvaluations, ":>" },
- T{ "opt_time", evaluator.CostFunctionTime,  ":>" },
+ T{ "eval_cnt", callCount, ":>" },
+ T{ "res_eval", resEval, ":>" },
+ T{ "jac_eval", jacEval, ":>" },
+ T{ "opt_time", cfTime, ":>" },
  T{ "seed", config.Seed, ":>" },
  T{ "elapsed", elapsed, ":>"},
  };

cli/source/util.cpp

@@ -192,6 +192,8 @@ auto InitOptions(std::string const& name, std::string const& desc, int width) ->
  ("offspring-generator", "OffspringGenerator operator, with optional parameters separated by : (eg --offspring-generator brood:10:10)", cxxopts::value<std::string>()->default_value("basic"))
  ("reinserter", "Reinsertion operator merging offspring in the recombination pool back into the population", cxxopts::value<std::string>()->default_value("keep-best"))
  ("enable-symbols", "Comma-separated list of enabled symbols ("+symbols+")", cxxopts::value<std::string>())
+ ("local-search-probability", "Probability for local search", cxxopts::value<Operon::Scalar>()->default_value("1.0"))
+ ("lamarckian-probability", "Probability that the local search improvements are saved back into the chromosome", cxxopts::value<Operon::Scalar>()->default_value("1.0"))
  ("disable-symbols", "Comma-separated list of disabled symbols ("+symbols+")", cxxopts::value<std::string>())
  ("symbolic", "Operate in symbolic mode - no coefficient tuning or coefficient mutation", cxxopts::value<bool>()->default_value("false"))
  ("show-primitives", "Display the primitive set used by the algorithm")

compile_commands.json

@@ -0,0 +1 @@
+build/compile_commands.json

include/operon/algorithms/config.hpp

@@ -14,10 +14,12 @@ struct GeneticAlgorithmConfig {
  size_t PopulationSize;
  size_t PoolSize;
  size_t Seed; // random seed
- size_t TimeLimit; // time limit
- double CrossoverProbability;
- double MutationProbability;
- double Epsilon; // used when comparing fitness values
+ size_t TimeLimit{~std::size_t{0}}; // time limit
+ double CrossoverProbability{1.0};
+ double MutationProbability{0.25};
+ double LocalSearchProbability{1.0};
+ double LamarckianProbability{1.0};
+ double Epsilon{0}; // used when comparing fitness values
 };
 } // namespace Operon
 

include/operon/interpreter/dispatch_table.hpp

@@ -185,6 +185,8 @@ struct DispatchTable {
  }(Ext{});
 
 public:
+ using SupportedTypes = Tup;
+
  template<typename T>
  static constexpr typename Ext::value_type BatchSize = Sizes[TypeIndex<T>];