diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5c2ce39..90397cd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -24,7 +24,9 @@ find_package(Ceres CONFIG)
 find_package(Eigen3 REQUIRED)
 find_package(FastFloat REQUIRED)
 find_package(Threads REQUIRED)
+find_package(lbfgs REQUIRED)
 find_package(operon REQUIRED)
+find_package(outcome REQUIRED)
 find_package(pratt-parser REQUIRED)
 find_package(pybind11 REQUIRED)
 find_package(unordered_dense REQUIRED)
@@ -43,7 +45,6 @@ pybind11_add_module(
     pyoperon_pyoperon
     MODULE
     source/algorithm.cpp
-    source/autodiff.cpp
     source/benchmark.cpp
     source/creator.cpp
     source/crossover.cpp
@@ -129,8 +130,8 @@ target_link_libraries(pyoperon_pyoperon PRIVATE
 if (MSVC)
     target_compile_options(pyoperon_pyoperon PRIVATE "/std:c++latest")
 else ()
-    if (UNIX AND NOT MAC)
-        target_link_options(pyoperon_pyoperon PRIVATE "-Wl,--no-undefined")
+    if (UNIX AND NOT APPLE)
+        target_link_options(pyoperon_pyoperon PUBLIC "-Wl,--no-undefined")
     endif()
 endif()
 
diff --git a/example/operon-bindings.py b/example/operon-bindings.py
index af78e07..5f2dd1a 100644
--- a/example/operon-bindings.py
+++ b/example/operon-bindings.py
@@ -77,11 +77,12 @@
 crossover      = Operon.SubtreeCrossover(crossover_internal_probability, maxD, maxL)
 
 # define fitness evaluation
-interpreter    = Operon.Interpreter() # tree interpreter
+dtable         = Operon.DispatchTable()
 error_metric   = Operon.R2()          # use the coefficient of determination as fitness
-evaluator      = Operon.Evaluator(problem, interpreter, error_metric, True) # initialize evaluator, use linear scaling = True
+evaluator      = Operon.Evaluator(problem, dtable, error_metric, True) # initialize evaluator, use linear scaling = True
 evaluator.Budget = 1000 * 1000             # computational budget
-evaluator.LocalOptimizationIterations = 0  # number of local optimization iterations (coefficient tuning using gradient descent)
+
+optimizer      = Operon.Optimizer(dtable, problem, optimizer="lbfgs", likelihood="gaussian", iterations=10, batchsize=50)
 
 # define how new offspring are created
 generator      = Operon.BasicOffspringGenerator(evaluator, crossover, mutation, selector, selector)
diff --git a/example/operon-sklearn.py b/example/operon-sklearn.py
index 1a0f737..152d505 100644
--- a/example/operon-sklearn.py
+++ b/example/operon-sklearn.py
@@ -4,25 +4,30 @@
 import pandas as pd
 import numpy as np
 from sklearn.model_selection import train_test_split, cross_val_score
-from sklearn.metrics import r2_score, make_scorer
+from sklearn.metrics import r2_score, make_scorer, mean_squared_error
+from sklearn.ensemble import RandomForestRegressor
 
 from pyoperon.sklearn import SymbolicRegressor
 from pyoperon import R2, MSE, InfixFormatter, FitLeastSquares, Interpreter
 
-df_train = pd.read_csv('/home/bogdb/src/poetryenv/notebooks/sr-workshop/postprocessing/data/stage1/data/3946_extrapolation_easy_data_train.csv')
-df_test = pd.read_csv('/home/bogdb/src/poetryenv/notebooks/sr-workshop/postprocessing/data/stage1/data/3946_extrapolation_easy_data_train.csv')
+df = pd.read_csv('/home/bogdb/src/operon/data/Poly-10.csv')
+X = df.iloc[:,:-1]
+y = df.iloc[:, -1]
+X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.5)
 
-print(df_train.columns)
+y_pred = RandomForestRegressor(n_estimators=100).fit(X_train, y_train).predict(X_train)
+sErr = np.sqrt(mean_squared_error(y_train,  y_pred))
 
-D_train = np.asarray(df_train)
-D_test = np.asarray(df_test)
+# print(df_train.columns)
 
-X_train, y_train = D_train[:,:-1], D_train[:,-1]
-X_test, y_test = D_train[:,:-1], D_train[:,-1]
+# D_train = np.asarray(df_train)
+# D_test = np.asarray(df_test)
+
+# X_train, y_train = D_train[:,:-1], D_train[:,-1]
+# X_test, y_test = D_train[:,:-1], D_train[:,-1]
 
 from sympy import parse_expr
 import matplotlib.pyplot as plt
-import seaborn as sns
 
 reg = SymbolicRegressor(
         allowed_symbols= "add,sub,mul,div,constant,variable",
@@ -37,7 +42,8 @@
         initialization_max_length= 10,
         initialization_method= "btc",
         irregularity_bias= 0.0,
-        local_iterations= 0,
+        local_iterations= 5,
+        optimizer='lm',
         male_selector= "tournament",
         max_depth= 10,
         max_evaluations= 1000000,
@@ -54,26 +60,16 @@
         reinserter= "keep-best",
         time_limit= 900,
         tournament_size= 3,
+        uncertainty= [sErr]
         )
 
 print(X_train.shape, y_train.shape)
 
 reg.fit(X_train, y_train)
-values = [s['objective_values'] for s in reg.pareto_front_]
-for v in values:
-    print(v)
+res = [(s['objective_values'], s['tree'], s['minimum_description_length']) for s in reg.pareto_front_]
+for obj, expr, mdl in res:
+    print(obj, mdl, reg.get_model_string(expr, 16))
 
 m = reg.model_
-s = reg.get_model_string(m, 3, ['a'])
+s = reg.get_model_string(m, 3)
 print(s)
-
-
-fig, ax = plt.subplots(figsize=(18,8))
-ax.grid(True, linestyle='dotted')
-ax.set(xlabel='Obj 1', ylabel='Obj 2')
-sns.scatterplot(ax=ax, x=[x[1] for x in values], y=[x[0] for x in values])
-
-from pyoperon import RankSorter 
-rs = RankSorter()
-fronts = rs.Sort(values)
-print(fronts)
diff --git a/flake.nix b/flake.nix
index cb8ae5a..acec941 100644
--- a/flake.nix
+++ b/flake.nix
@@ -24,7 +24,7 @@
         stdenv_ = pkgs.overrideCC pkgs.llvmPackages_16.stdenv (
           pkgs.clang_16.override { gccForLibs = pkgs.gcc13.cc; }
         );
-        python = pkgs.python310;
+        python_ = pkgs.python310;
 
         operon = pkgs.callPackage ./nix/operon {
           enableShared = enableShared;
@@ -49,14 +49,14 @@
             cmake
             ninja
             pkg-config
-            python
-            python.pkgs.pybind11
+            python_
+            python_.pkgs.pybind11
           ];
 
           buildInputs = with pkgs; [
-            python.pkgs.setuptools
-            python.pkgs.wheel
-            python.pkgs.requests
+            python_.pkgs.setuptools
+            python_.pkgs.wheel
+            python_.pkgs.requests
             operon
           ] ++ operon.buildInputs;
         };
@@ -78,10 +78,10 @@
         devShells.default = stdenv_.mkDerivation {
           name = "pyoperon-dev";
           nativeBuildInputs = pyoperon.nativeBuildInputs;
-          buildInputs = pyoperon.buildInputs ++ (with pkgs; [ gdb valgrind ])
-                          ++ (with python.pkgs; [ scikit-build ] ) # cmake integration and release preparation
-                          ++ (with python.pkgs; [ numpy scikit-learn pandas ipdb sympy requests ])
-                          ++ (with pkgs; [ (pmlb.override { pythonPackages = python.pkgs; }) ]);
+          buildInputs = pyoperon.buildInputs ++ (with pkgs; [ gdb valgrind gcc13 ])
+                          ++ (with python_.pkgs; [ scikit-build ] ) # cmake integration and release preparation
+                          ++ (with python_.pkgs; [ numpy scikit-learn pandas ipdb sympy requests matplotlib ])
+                          ++ (with pkgs; [ (pmlb.override { pythonPackages = python_.pkgs; }) ]);
         };
 
         # backwards compatibility
diff --git a/nix/operon/default.nix b/nix/operon/default.nix
index d95b5e2..3a34d6c 100644
--- a/nix/operon/default.nix
+++ b/nix/operon/default.nix
@@ -36,11 +36,13 @@ stdenv.mkDerivation rec {
   pname = "operon";
   version = "0.3.1";
 
+  #src = /home/bogdb/src/operon-mdl-fix;
+
   src = fetchFromGitHub {
     owner = "heal-research";
     repo = "operon";
-    rev = "0e359494b0239f4427d9518097bb304641ea7990";
-    sha256 = "sha256-WXV295CF1fqG7KrOS+uSi7S7+yLAQ6r44rU+RTGivxA=";
+    rev = "382b68d83a6c693dca5852d251e7991250f09b3f";
+    hash = "sha256-zUlQ5bPjHWim7XA3JXYERiwM1YsA97dSPbRBBy2XD5Y=";
   };
 
   nativeBuildInputs = [ cmake git ];
@@ -52,10 +54,14 @@ stdenv.mkDerivation rec {
     eve
     fast_float
     git
+    lbfgs
+    ned14-outcome
+    ned14-quickcpplib
+    ned14-status-code
     pkg-config
     pratt-parser
-    unordered_dense
     taskflow
+    unordered_dense
     vstat
     xxHash
     (scnlib.override { enableShared = enableShared; })
diff --git a/pyoperon/sklearn.py b/pyoperon/sklearn.py
index 6382c8b..98bc209 100644
--- a/pyoperon/sklearn.py
+++ b/pyoperon/sklearn.py
@@ -38,6 +38,10 @@ def __init__(self,
         offspring_generator            = 'basic',
         reinserter                     = 'replace-worst',
         objectives                     = ['r2'],
+        optimizer                      = 'lbfgs',
+        optimizer_likelihood           = 'gaussian',
+        optimizer_likelihood_loginput  = False,
+        optimizer_batch_size           = 0,
         max_length                     = 50,
         max_depth                      = 10,
         initialization_method          = 'btc',
@@ -57,6 +61,7 @@ def __init__(self,
         irregularity_bias              = 0.0,
         epsilon                        = 1e-5,
         model_selection_criterion      = 'minimum_description_length',
+        uncertainty                    = [1],
         n_threads                      = 1,
         time_limit                     = None,
         random_state                   = None
@@ -72,6 +77,10 @@ def __init__(self,
         self.offspring_generator       = offspring_generator
         self.reinserter                = reinserter
         self.objectives                = objectives
+        self.optimizer                 = optimizer
+        self.optimizer_likelihood      = optimizer_likelihood
+        self.optimizer_likelihood_loginput = optimizer_likelihood_loginput
+        self.optimizer_batch_size      = optimizer_batch_size
         self.max_length                = max_length
         self.max_depth                 = max_depth
         self.initialization_method     = initialization_method
@@ -92,6 +101,7 @@ def __init__(self,
         self.epsilon                   = epsilon
         self.n_threads                 = n_threads
         self.model_selection_criterion = model_selection_criterion
+        self.uncertainty               = uncertainty
         self.time_limit                = time_limit
         self.random_state              = random_state
 
@@ -107,6 +117,10 @@ def __check_parameters(self):
         self.offspring_generator            = check(self.offspring_generator, 'basic')
         self.reinserter                     = check(self.reinserter, 'replace-worst')
         self.objectives                     = check(self.objectives, [ 'r2' ])
+        self.optimizer                      = check(self.optimizer, 'lbfgs')
+        self.optimizer_likelihood           = check(self.optimizer_likelihood, 'gaussian')
+        self.optimizer_likelihood_loginput  = check(self.optimizer_likelihood_loginput, False)
+        self.optimizer_batch_size           = check(self.optimizer_batch_size, 0)
         self.max_length                     = check(self.max_length, 50)
         self.max_depth                      = check(self.max_depth, 10)
         self.initialization_method          = check(self.initialization_method, 'btc')
@@ -126,6 +140,7 @@ def __check_parameters(self):
         self.irregularity_bias              = check(self.irregularity_bias, 0.0)
         self.epsilon                        = check(self.epsilon, 1e-5)
         self.model_selection_criterion      = check(self.model_selection_criterion, 'minimum_description_length')
+        self.uncertainty                    = check(self.uncertainty, [1])
         self.n_threads                      = check(self.n_threads, 1)
         self.time_limit                     = check(self.time_limit, sys.maxsize)
         self.random_state                   = check(self.random_state, random.getrandbits(64))
@@ -207,30 +222,30 @@ def __init_selector(self, selection_method, comp):
         raise ValueError('Unknown selection method {}'.format(selection_method))
 
 
-    def __init_evaluator(self, objective, problem, interpreter):
+    def __init_evaluator(self, objective, problem, dtable):
         if objective == 'r2':
             err = op.R2()
-            return op.Evaluator(problem, interpreter, err, True), err
+            return op.Evaluator(problem, dtable, err, True), err
 
         elif objective == 'c2':
             err = op.C2()
-            return op.Evaluator(problem, interpreter, err, False), err
+            return op.Evaluator(problem, dtable, err, False), err
 
         elif objective == 'nmse':
             err = op.NMSE()
-            return op.Evaluator(problem, interpreter, err, True), err
+            return op.Evaluator(problem, dtable, err, True), err
 
         elif objective == 'rmse':
             err = op.RMSE()
-            return op.Evaluator(problem, interpreter, err, True), err
+            return op.Evaluator(problem, dtable, err, True), err
 
         elif objective == 'mse':
             err = op.MSE()
-            return op.Evaluator(problem, interpreter, err, True), err
+            return op.Evaluator(problem, dtable, err, True), err
 
         elif objective == 'mae':
             err = op.MAE()
-            return op.Evaluator(problem, interpreter, err, True), err
+            return op.Evaluator(problem, dtable, err, True), err
 
         elif objective == 'length':
             return op.LengthEvaluator(problem), None
@@ -241,7 +256,7 @@ def __init_evaluator(self, objective, problem, interpreter):
         elif objective == 'diversity':
             return op.DiversityEvaluator(problem), None
 
-        raise ValueError('Unknown objective {}'.format(objectives))
+        raise ValueError('Unknown objective {}'.format(objective))
 
 
     def __init_generator(self, generator_name, evaluator, crossover, mutator, female_selector, male_selector):
@@ -375,31 +390,38 @@ def fit(self, X, y):
 
         single_objective      = True if len(self.objectives) == 1 else False
 
-        interpreter = op.Interpreter()
+        dtable = op.DispatchTable()
 
         # these lists are used as placeholders in order to extend the lifetimes of the objects
         error_metrics = [] # placeholder for the error metric
         evaluators = [] # placeholder for the evaluator(s)
 
+        optimizer = op.Optimizer(dtable=dtable, problem=problem, optimizer=self.optimizer, likelihood=self.optimizer_likelihood, iterations=self.local_iterations, batchsize=self.optimizer_batch_size, loginput=self.optimizer_likelihood_loginput)
+
         # evaluators for minimum description length and information criteria
-        mld_eval = op.MinimumDescriptionLengthEvaluator(problem, interpreter)
-        bic_eval = op.BayesianInformationCriterionEvaluator(problem, interpreter)
-        aik_eval = op.AkaikeInformationCriterionEvaluator(problem, interpreter)
+        mdl_eval = op.MinimumDescriptionLengthEvaluator(problem, dtable)
+        mdl_eval.Sigma = self.uncertainty
+
+        bic_eval = op.BayesianInformationCriterionEvaluator(problem, dtable)
+        aik_eval = op.AkaikeInformationCriterionEvaluator(problem, dtable)
+
+        for eval in [mdl_eval, bic_eval, aik_eval]:
+            eval.Optimizer = optimizer
 
         for obj in self.objectives:
-            eval_, err_  = self.__init_evaluator(obj, problem, interpreter)
+            eval_, err_  = self.__init_evaluator(obj, problem, dtable)
             eval_.Budget = self.max_evaluations
-            eval_.LocalOptimizationIterations = self.local_iterations
             evaluators.append(eval_)
             error_metrics.append(err_)
 
+        evaluators[0].Optimizer = optimizer
+
         if single_objective:
             evaluator = evaluators[0]
         else:
             evaluator = op.MultiEvaluator(problem)
             for eval_ in evaluators:
                 evaluator.Add(eval_)
-            evaluator.LocalOptimizationIterations = self.local_iterations
             evaluator.Budget = self.max_evaluations
 
         comparison            = op.SingleObjectiveComparison(0) if single_objective else op.CrowdedComparison()
@@ -453,9 +475,13 @@ def fit(self, X, y):
         def get_solution_stats(solution):
             """Takes a solution (operon individual) and computes a set of stats"""
             # perform linear scaling
-            y_pred = op.Evaluate(interpreter, solution.Genotype, ds, training_range)
+            y_pred = op.Evaluate(dtable, solution.Genotype, ds, training_range)
             scale, offset = op.FitLeastSquares(y_pred, y)
-            nodes = solution.Genotype.Nodes + [ op.Node.Constant(scale), op.Node.Mul(), op.Node.Constant(offset), op.Node.Add() ]
+            nodes = solution.Genotype.Nodes
+            if scale != 1:
+                nodes += [ op.Node.Constant(scale), op.Node.Mul() ]
+            if offset != 0:
+                nodes += [ op.Node.Constant(offset), op.Node.Add() ]
             solution.Genotype = op.Tree(nodes).UpdateNodes()
 
             # get solution variables
@@ -467,7 +493,7 @@ def get_solution_stats(solution):
                 'tree' : solution.Genotype,
                 'objective_values' : evaluator(rng, solution),
                 'mean_squared_error' : mean_squared_error(y, scale * y_pred + offset),
-                'minimum_description_length' : mld_eval(rng, solution)[0],
+                'minimum_description_length' : mdl_eval(rng, solution)[0],
                 'bayesian_information_criterion' : bic_eval(rng, solution)[0],
                 'akaike_information_criterion' : aik_eval(rng, solution)[0],
             }
diff --git a/source/autodiff.cpp b/source/autodiff.cpp
deleted file mode 100644
index 3b773db..0000000
--- a/source/autodiff.cpp
+++ /dev/null
@@ -1,56 +0,0 @@
-
-// SPDX-License-Identifier: MIT
-// SPDX-FileCopyrightText: Copyright 2019-2021 Heal Research
-
-#include <pybind11/numpy.h>
-#include <pybind11/pybind11.h>
-#include <pybind11/stl.h>
-
-#include <operon/interpreter/interpreter.hpp>
-#include <operon/autodiff/autodiff.hpp>
-#include "pyoperon/pyoperon.hpp"
-
-namespace py = pybind11;
-
-
-void InitAutodiff(py::module_ &m)
-{
-    using ReverseDerivativeCalculator = Operon::Autodiff::DerivativeCalculator<Operon::Interpreter, Operon::Autodiff::AutodiffMode::Reverse>;
-    using ForwardDerivativeCalculator = Operon::Autodiff::DerivativeCalculator<Operon::Interpreter, Operon::Autodiff::AutodiffMode::Forward>;
-
-    py::class_<ReverseDerivativeCalculator>(m, "ReverseDerivativeCalculator")
-        .def(py::init([](Operon::Interpreter const& interpreter) {
-                return ReverseDerivativeCalculator(interpreter, false);
-            }))
-        .def("__call__", [](ReverseDerivativeCalculator const& self, Operon::Tree const& tree, Operon::Dataset const& dataset, Operon::Range const range) {
-                auto coeff = tree.GetCoefficients();
-                auto result = py::array_t<Operon::Scalar>(static_cast<pybind11::ssize_t>(range.Size() * coeff.size()));
-                auto span = MakeSpan(result);
-                self.operator()<Eigen::RowMajor>(tree, dataset, range, { coeff }, span);
-                return result;
-            })
-        .def("__call__", [](ReverseDerivativeCalculator const& self, Operon::Tree const& tree, Operon::Dataset const& dataset, Operon::Range const range, py::array_t<Operon::Scalar> result) {
-                auto coeff = tree.GetCoefficients();
-                auto span = MakeSpan(result);
-                self.operator()<Eigen::RowMajor>(tree, dataset, range, { coeff }, span);
-            })
-        ;
-
-    py::class_<ForwardDerivativeCalculator>(m, "ForwardDerivativeCalculator")
-        .def(py::init([](Operon::Interpreter const& interpreter) {
-                return ForwardDerivativeCalculator(interpreter, false);
-            }))
-        .def("__call__", [](ForwardDerivativeCalculator const& self, Operon::Tree const& tree, Operon::Dataset const& dataset, Operon::Range const range) {
-                auto coeff = tree.GetCoefficients();
-                auto result = py::array_t<Operon::Scalar>(static_cast<pybind11::ssize_t>(range.Size() * coeff.size()));
-                auto span = MakeSpan(result);
-                self.operator()<Eigen::RowMajor>(tree, dataset, range, { coeff }, span);
-                return result;
-            })
-        .def("__call__", [](ForwardDerivativeCalculator const& self, Operon::Tree const& tree, Operon::Dataset const& dataset, Operon::Range const range, py::array_t<Operon::Scalar> result) {
-                auto coeff = tree.GetCoefficients();
-                auto span = MakeSpan(result);
-                self.operator()<Eigen::RowMajor>(tree, dataset, range, { coeff }, span);
-            })
-        ;
-}
diff --git a/source/benchmark.cpp b/source/benchmark.cpp
index 07615db..297603f 100644
--- a/source/benchmark.cpp
+++ b/source/benchmark.cpp
@@ -15,6 +15,9 @@
 
 #include "pyoperon/pyoperon.hpp"
 
+using TDispatch    = Operon::DispatchTable<Operon::Scalar>;
+using TInterpreter = Operon::Interpreter<Operon::Scalar, TDispatch>;
+
 void InitBenchmark(py::module_ &m)
 {
     // benchmark functionality
@@ -48,15 +51,16 @@ void InitBenchmark(py::module_ &m)
         auto nTotal = std::reduce(trees.begin(), trees.end(), 0UL, [](size_t partial, const auto& t) { return partial + t.Length(); });
 #else
         auto nTotal = std::transform_reduce(trees.begin(), trees.end(), 0UL, std::plus<> {}, [](auto& tree) { return tree.Length(); });
+
+        TDispatch dtable;
 #endif
-        Operon::Interpreter interpreter;
         tf::Executor executor(nThreads);
         std::vector<std::vector<Operon::Scalar>> values(executor.num_workers());
         for (auto& val : values) { val.resize(range.Size()); }
         tf::Taskflow taskflow;
         taskflow.for_each(trees.begin(), trees.end(), [&](auto const& tree) {
             auto& val = values[executor.this_worker_id()];
-            interpreter.operator()<Operon::Scalar>(tree, ds, range, val);
+            TInterpreter{dtable, ds, tree}.Evaluate({}, range, val);
         });
         executor.run(taskflow).wait();
         return nTotal * range.Size();
diff --git a/source/creator.cpp b/source/creator.cpp
index 582ce8d..3ee2881 100644
--- a/source/creator.cpp
+++ b/source/creator.cpp
@@ -12,21 +12,24 @@ void InitCreator(py::module_ &m)
     py::class_<Operon::CreatorBase> creatorBase(m, "CreatorBase");
 
     py::class_<Operon::BalancedTreeCreator, Operon::CreatorBase>(m, "BalancedTreeCreator")
-        .def(py::init([](Operon::PrimitiveSet const& grammar, std::vector<Operon::Hash> const& variables, double bias) {
-            return Operon::BalancedTreeCreator(grammar, Operon::Span<Operon::Hash const>(variables.data(), variables.size()), bias);
-        }), py::arg("grammar"), py::arg("variables"), py::arg("bias"))
+        .def(py::init<Operon::PrimitiveSet const&, std::vector<Operon::Hash>, double>(),
+                py::arg("grammar"),
+                py::arg("variables"),
+                py::arg("bias"))
         .def("__call__", &Operon::BalancedTreeCreator::operator())
         .def_property("IrregularityBias", &Operon::BalancedTreeCreator::GetBias, &Operon::BalancedTreeCreator::SetBias);
 
     py::class_<Operon::ProbabilisticTreeCreator, Operon::CreatorBase>(m, "ProbabilisticTreeCreator")
-        .def(py::init([](Operon::PrimitiveSet const& grammar, std::vector<Operon::Hash> const& variables, double bias) {
-            return Operon::ProbabilisticTreeCreator(grammar, Operon::Span<Operon::Hash const>(variables.data(), variables.size()), bias);
-        }), py::arg("grammar"), py::arg("variables"), py::arg("bias"))
-        .def("__call__", &Operon::ProbabilisticTreeCreator::operator());
+        .def(py::init<Operon::PrimitiveSet const&, std::vector<Operon::Hash>, double>(),
+                py::arg("grammar"),
+                py::arg("variables"),
+                py::arg("bias"))
+        .def("__call__", &Operon::ProbabilisticTreeCreator::operator())
+        .def_property("IrregularityBias", &Operon::ProbabilisticTreeCreator::GetBias, &Operon::ProbabilisticTreeCreator::SetBias);
 
     py::class_<Operon::GrowTreeCreator, Operon::CreatorBase>(m, "GrowTreeCreator")
-        .def(py::init([](Operon::PrimitiveSet const& grammar, std::vector<Operon::Hash> const& variables) {
-            return Operon::GrowTreeCreator(grammar, Operon::Span<Operon::Hash const>(variables.data(), variables.size()));
-        }), py::arg("grammar"), py::arg("variables"))
+        .def(py::init<Operon::PrimitiveSet const&, std::vector<Operon::Hash>>(),
+                py::arg("grammar"),
+                py::arg("variables"))
         .def("__call__", &Operon::GrowTreeCreator::operator());
 }
diff --git a/source/eval.cpp b/source/eval.cpp
index 5342566..70c82d6 100644
--- a/source/eval.cpp
+++ b/source/eval.cpp
@@ -10,28 +10,145 @@
 
 namespace py = pybind11;
 
+using TDispatch        = Operon::DefaultDispatch;
+using TInterpreter     = Operon::Interpreter<Operon::Scalar, TDispatch>;
+using TInterpreterBase = Operon::InterpreterBase<Operon::Scalar>;
+
+using TEvaluatorBase   = Operon::EvaluatorBase;
+using TEvaluator       = Operon::Evaluator<TDispatch>;
+using TMDLEvaluator    = Operon::MinimumDescriptionLengthEvaluator<TDispatch>;
+using TBICEvaluator    = Operon::BayesianInformationCriterionEvaluator<TDispatch>;
+using TAIKEvaluator    = Operon::AkaikeInformationCriterionEvaluator<TDispatch>;
+
+// likelihood
+using TGaussianLikelihood      = Operon::GaussianLikelihood<Operon::Scalar>;
+using TPoissonLikelihood       = Operon::PoissonLikelihood<Operon::Scalar, false>;
+using TPoissonLikelihoodLog    = Operon::PoissonLikelihood<Operon::Scalar, true>;
+
+// optimizer
+using TOptimizerBase           = Operon::OptimizerBase<TDispatch>;
+using TLMOptimizerTiny         = Operon::LevenbergMarquardtOptimizer<TDispatch, Operon::OptimizerType::Tiny>;
+using TLMOptimizerEigen        = Operon::LevenbergMarquardtOptimizer<TDispatch, Operon::OptimizerType::Eigen>;
+using TLMOptimizerCeres        = Operon::LevenbergMarquardtOptimizer<TDispatch, Operon::OptimizerType::Ceres>;
+
+using TLBGSOptimizerGauss      = Operon::LBFGSOptimizer<TDispatch, TGaussianLikelihood>;
+using TLBGSOptimizerPoisson    = Operon::LBFGSOptimizer<TDispatch, TPoissonLikelihood>;
+using TLBGSOptimizerPoissonLog = Operon::LBFGSOptimizer<TDispatch, TPoissonLikelihoodLog>;
+
 namespace detail {
-    template<typename T>
-    auto FitLeastSquares(py::array_t<T> lhs, py::array_t<T> rhs) -> std::pair<double, double>
-    {
-        auto s1 = MakeSpan(lhs);
-        auto s2 = MakeSpan(rhs);
-        return Operon::FitLeastSquares(s1, s2);
-    }
+template<typename T>
+auto FitLeastSquares(py::array_t<T> lhs, py::array_t<T> rhs) -> std::pair<double, double>
+{
+    auto s1 = MakeSpan(lhs);
+    auto s2 = MakeSpan(rhs);
+    return Operon::FitLeastSquares(s1, s2);
+}
+
+class Optimizer {
+    std::unique_ptr<TOptimizerBase> impl_;
+
+    public:
+        Optimizer(TDispatch const& dtable, Operon::Problem const& problem, std::string const& optName, std::string const& likName, std::size_t iter, std::size_t bsize, bool loginput = true) {
+            if (optName == "lm") {
+                impl_ = std::make_unique<TLMOptimizerEigen>(dtable, problem);
+            } else if (optName == "lbfgs") {
+                if (likName == "gaussian") {
+                    impl_ = std::make_unique<TLBGSOptimizerGauss>(dtable, problem);
+                } else if (likName == "poisson") {
+                    if (loginput) {
+                        impl_ = std::make_unique<TLBGSOptimizerPoissonLog>(dtable, problem);
+                    } else {
+                        impl_ = std::make_unique<TLBGSOptimizerPoisson>(dtable, problem);
+                    }
+                }
+            } else if (optName == "sgd") {
+                // TODO: add SGD optimizer
+            }
+            impl_->SetIterations(iter);
+            impl_->SetBatchSize(bsize);
+        }
+
+        auto SetBatchSize(std::size_t value) const { impl_->SetBatchSize(value); }
+        [[nodiscard]] auto BatchSize() const { return impl_->BatchSize(); }
+
+        auto SetIterations(std::size_t value) const { impl_->SetIterations(value); }
+        [[nodiscard]] auto Iterations() const { return impl_->Iterations(); }
+
+        [[nodiscard]] auto GetDispatchTable() const { return impl_->GetDispatchTable(); }
+        [[nodiscard]] auto GetProblem() const { return impl_->GetProblem(); }
+
+        [[nodiscard]] auto Optimize(Operon::RandomGenerator& rng, Operon::Tree const& tree) const {
+            return impl_->Optimize(rng, tree);
+        }
+
+        [[nodiscard]] auto ComputeLikelihood(Operon::Span<Operon::Scalar const> x, Operon::Span<Operon::Scalar const> y, Operon::Span<Operon::Scalar const> w) const {
+            return impl_->ComputeLikelihood(x, y, w);
+        }
+
+        [[nodiscard]] auto ComputeFisherMatrix(Operon::Span<Operon::Scalar const> pred, Operon::Span<Operon::Scalar const> jac, Operon::Span<Operon::Scalar const> sigma) const {
+            return impl_->ComputeFisherMatrix(pred, jac, sigma);
+        }
+
+        [[nodiscard]] auto OptimizerImpl() const { return impl_.get(); }
+};
 } // namespace detail
 
+void InitOptimizer(py::module_ &m)
+{
+    using detail::Optimizer;
+    using Operon::Problem;
+    using std::string;
+    using std::size_t;
+
+    py::class_<Operon::OptimizerSummary>(m, "OptimizerSummary")
+        .def_readwrite("InitialCost", &Operon::OptimizerSummary::InitialCost)
+        .def_readwrite("FinalCost", &Operon::OptimizerSummary::FinalCost)
+        .def_readwrite("Iterations", &Operon::OptimizerSummary::Iterations)
+        .def_readwrite("FunctionEvaluations", &Operon::OptimizerSummary::FunctionEvaluations)
+        .def_readwrite("JacobianEvaluations", &Operon::OptimizerSummary::JacobianEvaluations)
+        .def_readwrite("Success", &Operon::OptimizerSummary::Success);
+
+    py::class_<detail::Optimizer>(m, "Optimizer")
+        .def(py::init<TDispatch const&, Problem const&, string const&, string const, size_t, size_t, bool>()
+                , py::arg("dtable")
+                , py::arg("problem")
+                , py::arg("optimizer") = "lbfgs"
+                , py::arg("likelihood") = "gaussian"
+                , py::arg("iterations") = 10
+                , py::arg("batchsize") = TDispatch::BatchSize<Operon::Scalar>
+                , py::arg("loginput") = false)
+        .def_property("BatchSize", &Optimizer::SetBatchSize, &Optimizer::BatchSize)
+        .def_property("Iterations", &Optimizer::SetIterations, &Optimizer::Iterations)
+        .def_property_readonly("DispatchTable", &Optimizer::GetDispatchTable)
+        .def_property_readonly("Problem", &Optimizer::GetProblem)
+        .def_property_readonly("OptimizerImpl", &Optimizer::OptimizerImpl)
+        .def("Optimize", &Optimizer::Optimize)
+        .def("ComputeLikelihood", &Optimizer::ComputeLikelihood)
+        .def("ComputeFisherMatrix", &Optimizer::ComputeFisherMatrix);
+}
+
 void InitEval(py::module_ &m)
 {
     // free functions
     // we use a lambda to avoid defining a fourth arg for the defaulted C++ function arg
-    m.def("Evaluate", [](Operon::Interpreter const& i, Operon::Tree const& t, Operon::Dataset const& d, Operon::Range r) {
+    m.def("Evaluate", [](Operon::Tree const& t, Operon::Dataset const& d, Operon::Range r) {
+        auto result = py::array_t<Operon::Scalar>(static_cast<pybind11::ssize_t>(r.Size()));
+        auto span = MakeSpan(result);
+        py::gil_scoped_release release;
+        TDispatch dtable;
+        TInterpreter{dtable, d, t}.Evaluate({}, r, span);
+        py::gil_scoped_acquire acquire;
+        return result;
+        }, py::arg("tree"), py::arg("dataset"), py::arg("range"));
+
+    m.def("Evaluate", [](TDispatch const& dtable, Operon::Tree const& t, Operon::Dataset const& d, Operon::Range r) {
         auto result = py::array_t<Operon::Scalar>(static_cast<pybind11::ssize_t>(r.Size()));
         auto span = MakeSpan(result);
         py::gil_scoped_release release;
-        i.operator()<Operon::Scalar>(t, d, r, span);
+        TInterpreter{dtable, d, t}.Evaluate({}, r, span);
         py::gil_scoped_acquire acquire;
         return result;
-        }, py::arg("interpreter"), py::arg("tree"), py::arg("dataset"), py::arg("range"));
+        }, py::arg("dtable"), py::arg("tree"), py::arg("dataset"), py::arg("range"));
 
     m.def("EvaluateTrees", [](std::vector<Operon::Tree> const& trees, Operon::Dataset const& ds, Operon::Range range, py::array_t<Operon::Scalar> result, size_t nthread) {
             auto span = MakeSpan(result);
@@ -40,10 +157,12 @@ void InitEval(py::module_ &m)
             py::gil_scoped_acquire acquire;
             }, py::arg("trees"), py::arg("dataset"), py::arg("range"), py::arg("result").noconvert(), py::arg("nthread") = 1);
 
-    m.def("CalculateFitness", [](Operon::Interpreter const& i, Operon::Tree const& t, Operon::Dataset const& d, Operon::Range r, std::string const& target, std::string const& metric) {
-        auto estimated = i.operator()<Operon::Scalar>(t, d, r);
+    m.def("CalculateFitness", [](Operon::Tree const& t, Operon::Dataset const& d, Operon::Range r, std::string const& target, std::string const& metric) {
         auto values = d.GetValues(target).subspan(r.Start(), r.Size());
 
+        TDispatch dtable;
+        auto estimated = TInterpreter{dtable, d, t}.Evaluate({}, r);
+
         if (metric == "c2") { return Operon::C2{}(estimated, values); }
         if (metric == "r2") { return Operon::R2{}(estimated, values); }
         if (metric == "mse") { return Operon::MSE{}(estimated, values); }
@@ -52,9 +171,9 @@ void InitEval(py::module_ &m)
         if (metric == "mae") { return Operon::MAE{}(estimated, values); }
         throw std::runtime_error("Invalid fitness metric"); 
 
-    }, py::arg("interpreter"), py::arg("tree"), py::arg("dataset"), py::arg("range"), py::arg("target"), py::arg("metric") = "rsquared");
+    }, py::arg("tree"), py::arg("dataset"), py::arg("range"), py::arg("target"), py::arg("metric") = "rsquared");
 
-    m.def("CalculateFitness", [](Operon::Interpreter const& i, std::vector<Operon::Tree> const& trees, Operon::Dataset const& d, Operon::Range r, std::string const& target, std::string const& metric) {
+    m.def("CalculateFitness", [](std::vector<Operon::Tree> const& trees, Operon::Dataset const& d, Operon::Range r, std::string const& target, std::string const& metric) {
         std::unique_ptr<Operon::ErrorMetric> error;
         if (metric == "c2") { error = std::make_unique<Operon::C2>(); }
         else if (metric == "r2") { error = std::make_unique<Operon::R2>(); }
@@ -64,18 +183,20 @@ void InitEval(py::module_ &m)
         else if (metric == "mae") { error = std::make_unique<Operon::MAE>(); }
         else { throw std::runtime_error("Unsupported error metric"); }
 
+        TDispatch dtable;
+
         auto result = py::array_t<double>(static_cast<pybind11::ssize_t>(trees.size()));
         auto buf = result.request();
         auto values = d.GetValues(target).subspan(r.Start(), r.Size());
 
         // TODO: make this run in parallel with taskflow
         std::transform(trees.begin(), trees.end(), static_cast<double*>(buf.ptr), [&](auto const& t) -> double {
-            auto estimated = i.operator()<Operon::Scalar>(t, d, r);
+            auto estimated = TInterpreter{dtable, d, t}.Evaluate({}, r);
             return (*error)(estimated, values);
         });
 
         return result;
-    }, py::arg("interpreter"), py::arg("trees"), py::arg("dataset"), py::arg("range"), py::arg("target"), py::arg("metric") = "rsquared");
+    }, py::arg("trees"), py::arg("dataset"), py::arg("range"), py::arg("target"), py::arg("metric") = "rsquared");
 
 
     m.def("FitLeastSquares", [](py::array_t<float> lhs, py::array_t<float> rhs) -> std::pair<double, double> {
@@ -86,16 +207,18 @@ void InitEval(py::module_ &m)
         return detail::FitLeastSquares<double>(lhs, rhs);
     });
 
-    using DispatchTable = Operon::DispatchTable<Operon::Scalar, Operon::Dual>;
-
     // dispatch table
-    py::class_<DispatchTable>(m, "DispatchTable")
+    py::class_<TDispatch>(m, "DispatchTable")
         .def(py::init<>());
 
+    py::class_<TInterpreterBase>(m, "InterpreterBase");
+
     // interpreter
-    py::class_<Operon::Interpreter>(m, "Interpreter")
-        .def(py::init<>())
-        .def(py::init<DispatchTable>());
+    py::class_<TInterpreter, TInterpreterBase>(m, "Interpreter")
+        .def(py::init<TDispatch const&, Operon::Dataset const&, Operon::Tree const&>())
+        .def("Evaluate", [](TInterpreter const& self, Operon::Span<Operon::Scalar const> coeff, Operon::Range range, Operon::Span<Operon::Scalar> result){
+            return self.Evaluate(coeff, range, result);
+        });
 
     // error metric
     py::class_<Operon::ErrorMetric>(m, "ErrorMetric")
@@ -111,18 +234,21 @@ void InitEval(py::module_ &m)
     py::class_<Operon::C2, Operon::ErrorMetric>(m, "C2").def(py::init<>());
 
     // evaluator
-    py::class_<Operon::EvaluatorBase>(m, "EvaluatorBase")
-        .def_property("LocalOptimizationIterations", &Operon::EvaluatorBase::LocalOptimizationIterations, &Operon::EvaluatorBase::SetLocalOptimizationIterations)
-        .def_property("Budget", &Operon::EvaluatorBase::Budget, &Operon::EvaluatorBase::SetBudget)
-        .def_property_readonly("TotalEvaluations", &Operon::EvaluatorBase::TotalEvaluations)
-        // .def("__call__", &Operon::EvaluatorBase::operator())
+    py::class_<TEvaluatorBase>(m, "EvaluatorBase")
+        .def_property("Budget", &TEvaluatorBase::Budget, &Operon::EvaluatorBase::SetBudget)
+        .def_property_readonly("TotalEvaluations", &TEvaluatorBase::TotalEvaluations)
+        //.def("__call__", &TEvaluatorBase::operator())
         .def("__call__", [](Operon::EvaluatorBase const& self, Operon::RandomGenerator& rng, Operon::Individual& ind) { return self(rng, ind, {}); })
-        .def_property_readonly("CallCount", [](Operon::EvaluatorBase& self) { return self.CallCount.load(); })
-        .def_property_readonly("ResidualEvaluations", [](Operon::EvaluatorBase& self) { return self.ResidualEvaluations.load(); })
-        .def_property_readonly("JacobianEvaluations", [](Operon::EvaluatorBase& self) { return self.JacobianEvaluations.load(); });
-
-    py::class_<Operon::Evaluator, Operon::EvaluatorBase>(m, "Evaluator")
-        .def(py::init<Operon::Problem&, Operon::Interpreter&, Operon::ErrorMetric const&, bool>());
+        .def_property_readonly("CallCount", [](TEvaluatorBase& self) { return self.CallCount.load(); })
+        .def_property_readonly("ResidualEvaluations", [](TEvaluatorBase& self) { return self.ResidualEvaluations.load(); })
+        .def_property_readonly("JacobianEvaluations", [](TEvaluatorBase& self) { return self.JacobianEvaluations.load(); });
+
+    py::class_<TEvaluator, Operon::EvaluatorBase>(m, "Evaluator")
+        .def(py::init<Operon::Problem&, TDispatch const&, Operon::ErrorMetric const&, bool>())
+        //.def_property("Optimizer", &TEvaluator::GetOptimizer, &TEvaluator::SetOptimizer);
+        .def_property("Optimizer", nullptr, [](TEvaluator& self, detail::Optimizer const& opt) {
+           self.SetOptimizer(opt.OptimizerImpl()); 
+        });
 
     py::class_<Operon::UserDefinedEvaluator, Operon::EvaluatorBase>(m, "UserDefinedEvaluator")
         .def(py::init<Operon::Problem&, std::function<typename Operon::EvaluatorBase::ReturnType(Operon::RandomGenerator*, Operon::Individual&)> const&>());
@@ -152,24 +278,14 @@ void InitEval(py::module_ &m)
         .def(py::init<Operon::EvaluatorBase&>())
         .def_property("AggregateType", &Operon::AggregateEvaluator::GetAggregateType, &Operon::AggregateEvaluator::SetAggregateType);
 
-    py::class_<Operon::MinimumDescriptionLengthEvaluator, Operon::Evaluator>(m, "MinimumDescriptionLengthEvaluator")
-        .def(py::init<Operon::Problem&, Operon::Interpreter&>())
-        .def_property("LocalOptimizationIterations",
-                &Operon::MinimumDescriptionLengthEvaluator::LocalOptimizationIterations,   // getter
-                &Operon::MinimumDescriptionLengthEvaluator::SetLocalOptimizationIterations<std::size_t> // setter
-        );
-
-    py::class_<Operon::BayesianInformationCriterionEvaluator, Operon::Evaluator>(m, "BayesianInformationCriterionEvaluator")
-        .def(py::init<Operon::Problem&, Operon::Interpreter&>())
-        .def_property("LocalOptimizationIterations",
-                &Operon::BayesianInformationCriterionEvaluator::LocalOptimizationIterations,   // getter
-                &Operon::BayesianInformationCriterionEvaluator::SetLocalOptimizationIterations<std::size_t> // setter
-        );
-
-    py::class_<Operon::AkaikeInformationCriterionEvaluator, Operon::Evaluator>(m, "AkaikeInformationCriterionEvaluator")
-        .def(py::init<Operon::Problem&, Operon::Interpreter&>())
-        .def_property("LocalOptimizationIterations",
-                &Operon::AkaikeInformationCriterionEvaluator::LocalOptimizationIterations,   // getter
-                &Operon::AkaikeInformationCriterionEvaluator::SetLocalOptimizationIterations<std::size_t> // setter
-        );
+    py::class_<TMDLEvaluator, TEvaluator>(m, "MinimumDescriptionLengthEvaluator")
+        .def(py::init<Operon::Problem&, TDispatch const&>())
+        .def_property("Sigma", nullptr /*get*/ , &TMDLEvaluator::SetSigma /*set*/);
+        // .def("__call__", &TMDLEvaluator::operator());
+
+    py::class_<TBICEvaluator, TEvaluator>(m, "BayesianInformationCriterionEvaluator")
+        .def(py::init<Operon::Problem&, TDispatch const&>());
+
+    py::class_<TAIKEvaluator, TEvaluator>(m, "AkaikeInformationCriterionEvaluator")
+        .def(py::init<Operon::Problem&, TDispatch const&>());
 }
diff --git a/source/optimizer.cpp b/source/optimizer.cpp
index f366f9d..7da5e96 100644
--- a/source/optimizer.cpp
+++ b/source/optimizer.cpp
@@ -1,43 +1,10 @@
 // SPDX-License-Identifier: MIT
 // SPDX-FileCopyrightText: Copyright 2019-2021 Heal Research
 
-#include <operon/autodiff/autodiff.hpp>
+#include <memory>
+#include <operon/optimizer/likelihood/gaussian_likelihood.hpp>
+#include <operon/optimizer/likelihood/poisson_likelihood.hpp>
 #include <operon/optimizer/optimizer.hpp>
 #include "pyoperon/pyoperon.hpp"
 
 namespace py = pybind11;
-
-void InitOptimizer(py::module_ &m)
-{
-    using DerivativeCalculator = Operon::Autodiff::DerivativeCalculator<Operon::Interpreter>;
-    using EigenOptimizer = Operon::NonlinearLeastSquaresOptimizer<DerivativeCalculator, Operon::OptimizerType::Eigen>;
-    using TinyOptimizer = Operon::NonlinearLeastSquaresOptimizer<DerivativeCalculator, Operon::OptimizerType::Tiny>;
-    using CeresOptimizer = Operon::NonlinearLeastSquaresOptimizer<DerivativeCalculator, Operon::OptimizerType::Ceres>;
-
-    py::class_<Operon::OptimizerSummary>(m, "OptimizerSummary")
-        .def_readwrite("InitialCost", &Operon::OptimizerSummary::InitialCost)
-        .def_readwrite("FinalCost", &Operon::OptimizerSummary::FinalCost)
-        .def_readwrite("Iterations", &Operon::OptimizerSummary::Iterations)
-        .def_readwrite("FunctionEvaluations", &Operon::OptimizerSummary::FunctionEvaluations)
-        .def_readwrite("JacobianEvaluations", &Operon::OptimizerSummary::JacobianEvaluations)
-        .def_readwrite("Success", &Operon::OptimizerSummary::Success);
-
-    py::class_<EigenOptimizer>(m, "Optimizer")
-        .def(py::init<DerivativeCalculator const&, Operon::Tree const&, Operon::Dataset const&>())
-        .def("Optimize", &EigenOptimizer::Optimize); 
-
-    m.def("Optimize", [](DerivativeCalculator const& dc, Operon::Tree const& tree, Operon::Dataset const& ds, py::array_t<Operon::Scalar const> target, Operon::Range range, std::size_t iterations) {
-        EigenOptimizer optimizer(dc, tree, ds);
-        Operon::OptimizerSummary summary{};
-        auto coeff = optimizer.Optimize(MakeSpan(target), range, iterations, summary);
-        return std::make_tuple(coeff, summary);
-    });
-
-    m.def("Optimize", [](Operon::Interpreter const& interpreter, Operon::Tree const& tree, Operon::Dataset const& ds, py::array_t<Operon::Scalar const> target, Operon::Range range, size_t iterations) {
-        DerivativeCalculator dc{interpreter};
-        EigenOptimizer optimizer(dc, tree, ds);
-        Operon::OptimizerSummary summary{};
-        auto coeff = optimizer.Optimize(MakeSpan(target), range, iterations, summary);
-        return std::make_tuple(coeff, summary);
-    });
-}
diff --git a/source/pyoperon.cpp b/source/pyoperon.cpp
index f4af7ea..2dd894d 100644
--- a/source/pyoperon.cpp
+++ b/source/pyoperon.cpp
@@ -20,7 +20,6 @@ PYBIND11_MODULE(pyoperon, m)
     py::bind_vector<std::vector<Operon::Individual>>(m, "IndividualCollection");
 
     InitAlgorithm(m);
-    InitAutodiff(m);
     InitBenchmark(m);
     InitCreator(m);
     InitCrossover(m);
diff --git a/source/tree.cpp b/source/tree.cpp
index d70a279..7bded56 100644
--- a/source/tree.cpp
+++ b/source/tree.cpp
@@ -25,7 +25,7 @@ void InitTree(py::module_ &m)
         .def_property_readonly("Nodes", static_cast<Operon::Vector<Operon::Node>& (Operon::Tree::*)()&>(&Operon::Tree::Nodes))
         .def_property_readonly("Nodes", static_cast<Operon::Vector<Operon::Node> const& (Operon::Tree::*)() const&>(&Operon::Tree::Nodes))
         //.def_property_readonly("Nodes", static_cast<Operon::Vector<Operon::Node>&& (Operon::Tree::*)() &&>(&Operon::Tree::Nodes))
-        .def_property_readonly("Indices", &Operon::Tree::Indices)
+        .def_property_readonly("Indices", [](Operon::Tree const& tree, std::size_t i) { return tree.Indices(i); })
         .def_property_readonly("Children", py::overload_cast<std::size_t>(&Operon::Tree::Children))
         .def_property_readonly("Children", py::overload_cast<std::size_t>(&Operon::Tree::Children, py::const_))
         .def_property_readonly("Length", &Operon::Tree::Length)