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>];