diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md index 8c41ce03..14c1c4d6 100644 --- a/.github/copilot-instructions.md +++ b/.github/copilot-instructions.md @@ -17,7 +17,6 @@ Bingo is an open-source Python package developed by NASA for symbolic regression - `sympy`: Symbolic mathematics - `scikit-learn>=1.1`: Machine learning utilities - `dill>=0.2.9`: Serialization -- `smcpy>=0.1.4`: Statistical methods - `pybind11`: Python-C++ bindings ### Build System @@ -34,7 +33,6 @@ bingo/ │ ├── evaluation/ # Fitness evaluation │ ├── evolutionary_algorithms/ # EA implementations │ ├── evolutionary_optimizers/ # Optimization strategies (Island, Archipelago) -│ ├── local_optimizers/ # Local optimization methods │ ├── selection/ # Selection operators (Tournament, etc.) │ ├── stats/ # Statistics and tracking │ ├── expressions/ # Expression representations and variation @@ -74,7 +72,7 @@ bingo/ - **Variation strategies**: - `VarOr`: Mutation OR crossover OR replication - `VarAnd`: Mutation AND crossover -- **Fitness evaluation**: Vector-based with aggregation metrics (MAE, RMSE, MSE, etc.) +- **Fitness evaluation**: Fitness functions return scalar lower-is-better values ### Code Organization Principles - Each module should have a clear, single responsibility @@ -216,8 +214,8 @@ pytest tests 6. Add usage examples if applicable ### Adding a New Fitness Function -1. Create class inheriting from `FitnessFunction` or `VectorBasedFunction` -2. Implement `evaluate_fitness_vector(self, individual)` method +1. Create a class inheriting from `FitnessFunction` +2. Implement `__call__(self, individual)` to return a scalar lower-is-better value 3. Add docstring explaining the metric 4. Add tests in `tests/unit/evaluation/` diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index d9bcd728..a4d51add 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -104,6 +104,26 @@ jobs: - name: Clean-install and import source distribution run: python -m pip install --upgrade pip && python -m pip install dist/*.tar.gz && python -c "import bingo" + smoke-test-evidence-extra: + name: Smoke test Evidence extra + needs: build + runs-on: ubuntu-latest + steps: + - uses: actions/download-artifact@v4 + with: + name: distributions + path: dist + - uses: actions/setup-python@v5 + with: + python-version: "3.13" + - name: Clean-install and import Evidence extra + shell: bash + run: | + python -m pip install --upgrade pip + wheel=$(printf '%s\n' dist/*cp313-*manylinux*.whl) + python -m pip install --only-binary=:all: "${wheel}[evidence]" + python -c "import bingo.symbolic_regression; import smcpy" + installed-tests: name: Installed artifact tests (Python ${{ matrix.python-version }}) needs: build @@ -114,8 +134,10 @@ jobs: include: - python-version: "3.13" wheel-tag: cp313 + dependencies: ".[MPI,ONNX,TESTS,evidence]" - python-version: "3.14" wheel-tag: cp314 + dependencies: ".[MPI,ONNX,TESTS,evidence]" steps: - uses: actions/checkout@v4 with: @@ -123,7 +145,7 @@ jobs: - uses: ./.github/actions/setup-bingo with: python-version: ${{ matrix.python-version }} - dependencies: ".[MPI,ONNX,TESTS]" + dependencies: ${{ matrix.dependencies }} - uses: actions/download-artifact@v4 with: name: distributions @@ -138,7 +160,7 @@ jobs: release-gate: name: Verify release tag and artifact version - needs: [metadata, smoke-test-wheel, smoke-test-sdist, installed-tests] + needs: [metadata, smoke-test-wheel, smoke-test-sdist, smoke-test-evidence-extra, installed-tests] if: startsWith(github.ref, 'refs/tags/') runs-on: ubuntu-latest outputs: diff --git a/README.md b/README.md index bea3465e..6fd32e1d 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ Bingo is an open source package for performing symbolic regression, though it can be used as a general purpose evolutionary optimization package. ## Key Features -* Integrated local optimization strategies +* Expression-specific symbolic-regression fitting * Parallel island evolution strategy implemented with mpi4py * Coevolution of fitness predictors diff --git a/bingo/chromosomes/chromosome.py b/bingo/chromosomes/chromosome.py index 4eb2981d..3b05d7e4 100644 --- a/bingo/chromosomes/chromosome.py +++ b/bingo/chromosomes/chromosome.py @@ -103,41 +103,3 @@ def distance(self, other): Distance from self to other """ raise NotImplementedError - - def needs_local_optimization(self): - """Does the `Chromosome` need local optimization - - Returns - ------- - bool - Whether `Chromosome` needs optimization - """ - raise NotImplementedError( - "This Chromosome cannot be used in local " - "optimization until its local optimization " - "interface has been implemented" - ) - - def get_number_local_optimization_params(self): - """Get number of parameters in local optimization - - Returns - ------- - int - Number of parameters to be optimized - """ - return 0 - - def set_local_optimization_params(self, params): - """Set local optimization parameters - - Parameters - ---------- - params : list-like of numeric - Values to set the parameters to - """ - raise NotImplementedError( - "This Chromosome cannot be used in local " - "optimization until its local optimization " - "interface has been implemented" - ) diff --git a/bingo/chromosomes/multiple_floats.py b/bingo/chromosomes/multiple_floats.py index 0695e446..356f1b41 100644 --- a/bingo/chromosomes/multiple_floats.py +++ b/bingo/chromosomes/multiple_floats.py @@ -14,60 +14,7 @@ class MultipleFloatChromosome(MultipleValueChromosome): ---------- values : list of floats The genetic information stored in an individual chromosome. - needs_opt_list : list of ints - (optional) The indices of the `individual_list` in a - `chromosomes` object that are subject local optimization. - This list may be empty """ - def __init__(self, values, needs_opt_list=None): - super().__init__(values) - if needs_opt_list is None: - needs_opt_list = [] - self._needs_opt_list = needs_opt_list - - def needs_local_optimization(self): - """Does the individual need local optimization - - Returns - ------- - bool - Individual needs optimization - """ - if not self._needs_opt_list: - return False - return True - - def get_number_local_optimization_params(self): - """Get number of parameters in local optimization - - Returns - ------- - int - number of parameters to be optimized - """ - return len(self._needs_opt_list) - - def set_local_optimization_params(self, params): - """Set local optimization parameters - - Parameters - ---------- - params : list-like of numeric - Values to set the parameters - """ - for param, index in zip(params, self._needs_opt_list): - self.values[index] = param - - def get_local_optimization_params(self): - """Gets local optimization_params - - Returns - ------- - list - """ - return [self.values[i] for i in self._needs_opt_list] - - class MultipleFloatChromosomeGenerator(MultipleValueChromosomeGenerator): """Generation of a population of Multi-Value chromosomes @@ -77,18 +24,7 @@ class MultipleFloatChromosomeGenerator(MultipleValueChromosomeGenerator): A function that returns a randomly generated float value. values_per_chromosome : int The number of values that each chromosome will hold - needs_opt_list : list of ints - The indices of the `individual_list` in a `chromosomes` object - that are subject local optimization. This list may be empty """ - def __init__(self, random_value_function, values_per_chromosome, - needs_opt_list=None): - super().__init__(random_value_function, values_per_chromosome) - if needs_opt_list is None: - needs_opt_list = [] - self._check_opt_list_contains_feasible_values(needs_opt_list) - self._needs_opt_list = self._remove_duplicates(needs_opt_list) - def __call__(self): """Generation of a population of size `population_size` of Multi-Value chromosomes with lists that contain @@ -99,21 +35,7 @@ def __call__(self): list of chromosomes : The chromosomes which their values are generated by `random_value_function` with the optimization list - `needs_opt_list`. + generated by `random_value_function`. """ random_list = self._generate_list(self._values_per_chromosome) - return MultipleFloatChromosome(random_list, self._needs_opt_list) - - def _check_opt_list_contains_feasible_values(self, list_of_indices): - if not all(isinstance(x, int) for x in list_of_indices): - raise ValueError("The list of optimization indices must be \ - unsigned integers.") - if list_of_indices and (min(list_of_indices) < 0 or - max(list_of_indices) >= self._values_per_chromosome): - raise ValueError("The list of optimization indices must be within \ - the length of the list of values.") - - @staticmethod - def _remove_duplicates(list_of_ints): - set_of_ints = set(list_of_ints) - return sorted(list(set_of_ints)) + return MultipleFloatChromosome(random_list) diff --git a/bingo/chromosomes/multiple_values.py b/bingo/chromosomes/multiple_values.py index fced2a5e..047ce4cd 100644 --- a/bingo/chromosomes/multiple_values.py +++ b/bingo/chromosomes/multiple_values.py @@ -55,16 +55,6 @@ def distance(self, other): dist = sum(v1 != v2 for v1, v2 in zip(self.values, other.values)) return dist - def get_number_local_optimization_params(self): - raise NotImplementedError - - def needs_local_optimization(self): - raise NotImplementedError - - def set_local_optimization_params(self, params): - raise NotImplementedError - - class MultipleValueChromosomeGenerator(Generator): """Generation of a population of Multi-Value chromosomes diff --git a/bingo/evaluation/fitness_function.py b/bingo/evaluation/fitness_function.py index fe621c0a..32bdbce6 100644 --- a/bingo/evaluation/fitness_function.py +++ b/bingo/evaluation/fitness_function.py @@ -1,78 +1,7 @@ -"""The definition of fitness evaluations for individuals. - -This module defines the basis of fitness evaluation in bingo evolutionary -analyses. It defines a fitness function class and a version of the fitness -function which is built around the idea that there is a vector of fitness -values that can be aggregated in numerous different ways. -""" +"""The definition of fitness evaluations for individuals.""" from abc import ABCMeta, abstractmethod -import numpy as np - - -# Fitness metric functions, outside of FitnessFunction for use in GradientMixin -def mean_absolute_error(vector, individual=None): # pylint: disable=unused-argument - """Calculate the mean absolute error of an error vector""" - return np.mean(np.abs(vector)) - - -def root_mean_squared_error(vector, individual=None): # pylint: disable=unused-argument - """Calculate the root mean squared error of an error vector""" - return np.sqrt(np.mean(np.square(vector))) - - -def mean_squared_error(vector, individual=None): # pylint: disable=unused-argument - """Calculate the mean squared error of an error vector""" - return np.mean(np.square(vector)) - - -def negative_nmll_laplace(vector, individual): - """Calculate the negative normalized marginal log likelihood (NMLL) of an error vector - - The normalized marginal log likelihood is a Bayesian model selection criterion - that balances model fit with complexity. It approximates the marginal likelihood - by integrating over the posterior distribution of model parameters using a - Laplace approximation. - - The NMLL is calculated as: - NMLL = - k * ln(b)/2 - (1 - b) * ln(L̂) - - where: - - b = 1/sqrt(n) is a normalization factor - - L̂ = maximized value of the likelihood function - - k = number of model parameters estimated by the model - - n = number of data points - - This metric penalizes model complexity more heavily than BIC for small sample - sizes and provides a probabilistic framework for model comparison. Lower values - indicate better models. - """ - n = len(vector) - k = individual.get_number_local_optimization_params() + 1 - b = 1 / np.sqrt(n) - mse = np.mean(np.square(vector)) - log_like = -n / 2 * np.log(mse) - n / 2 - n / 2 * np.log(2 * np.pi) - nmll_laplace = (1 - b) * log_like + np.log(b) / 2 * k - return -nmll_laplace - - -def bic(vector, individual): - """Calculate the Bayesian Information Criterion (BIC) of an error vector - - BIC = k * ln(n) - 2 * ln(L̂) - - where: - - k = number of parameters estimated by the model - - n = number of data points - - L̂ = maximized value of the likelihood function - """ - n = len(vector) - k = individual.get_number_local_optimization_params() + 1 - mse = np.mean(np.square(vector)) - log_likelihood = - n / 2 * np.log(mse) - n / 2 - n / 2 * np.log(2 * np.pi) - return k * np.log(n) - 2 * log_likelihood - class FitnessFunction(metaclass=ABCMeta): """Fitness evaluation metric for individuals. @@ -117,72 +46,3 @@ def __call__(self, individual): fitness of the individual """ raise NotImplementedError - - -class VectorBasedFunction(FitnessFunction, metaclass=ABCMeta): - """Fitness evaluation based on vectorized fitness - - An aggregation metric is needed to quantify Fitness relative to a vector - of fitness (error) measures. - - Parameters - ---------- - training_data : TrainingData - data that is used in fitness evaluation. - metric : str - String defining the measure of error to use. Available options are: - 'mean absolute error'/'mae', 'mean squared error'/'mse', - 'root mean squared error'/'rmse', "negative nmll laplace", and "bic" - """ - - def __init__(self, training_data=None, metric="mae"): - super().__init__(training_data) - - if metric in ["mean absolute error", "mae"]: - self._metric = mean_absolute_error - elif metric in ["mean squared error", "mse"]: - self._metric = mean_squared_error - elif metric in ["root mean squared error", "rmse"]: - self._metric = root_mean_squared_error - elif metric in ["negative nmll laplace"]: - self._metric = negative_nmll_laplace - elif metric in ["bic"]: - self._metric = bic - else: - raise ValueError("Invalid metric for Fitness Function") - - def __call__(self, individual): - """Vector based fitness evaluation - - Evaluate the fitness of an individual as based on a vector of fitness - (error) values. The metric defined in the constructor is used to - aggregate the vector fitness into a single fitness value - - Parameters - ---------- - individual : Chromosome - individual for which fitness will be calculated - - Returns - ------- - fitness : numeric - fitness of the individual - """ - fitness_vector = self.evaluate_fitness_vector(individual) - return self._metric(fitness_vector, individual) - - @abstractmethod - def evaluate_fitness_vector(self, individual): - """Calculate a vector of fitness values for the passed in individual - - Parameters - ---------- - individual : Chromosome - individual for which fitness will be calculated - - Returns - ------- - vector_fitness : array of numeric - a vector of fitness values for the passed in individual - """ - raise NotImplementedError diff --git a/bingo/evaluation/gradient_mixin.py b/bingo/evaluation/gradient_mixin.py deleted file mode 100644 index 1e544566..00000000 --- a/bingo/evaluation/gradient_mixin.py +++ /dev/null @@ -1,174 +0,0 @@ -"""Mixin classes used to extend fitness functions to be able to use -gradient- and jacobian-based continuous local optimization methods. - -This module defines the basis of gradient and jacobian partial derivatives -of fitness functions used in bingo evolutionary analyses. -""" - -from abc import ABCMeta, abstractmethod -import numpy as np - -from .fitness_function import ( - mean_absolute_error, - mean_squared_error, - root_mean_squared_error, - negative_nmll_laplace, - bic, -) - - -class GradientMixin(metaclass=ABCMeta): - """Mixin for using gradients for fitness functions - - An abstract base class/mixin used to implement the gradients - of fitness functions. - """ - - @abstractmethod - def get_fitness_and_gradient(self, individual): - """Fitness function evaluation and gradient - - Get the fitness of the individual and the gradient - of this function with respect to the individual's constants. - - Parameters - ---------- - individual : Chromosome - individual for which the fitness and gradient will be calculated for - - Returns - ------- - fitness, gradient : - fitness of the individual and the gradient of this function - with respect to the individual's constants - """ - raise NotImplementedError - - -class VectorGradientMixin(GradientMixin): - """Mixin for using gradients and jacobians in vector based fitness functions - - An abstract base class/mixin used to implement the gradients and jacobians - of vector based fitness functions. - - Parameters - ---------- - training_data : TrainingData - data that is used in fitness evaluation (passed to parent). - metric : str - String defining the measure of error to use. Available options are: - 'mean absolute error', 'mean squared error', - 'root mean squared error', "negative nmll laplace", and "bic" - """ - - def __init__(self, training_data=None, metric="mae"): - super().__init__(training_data, metric) - - if metric in ["mean absolute error", "mae"]: - self._metric = mean_absolute_error - self._metric_derivative = ( - VectorGradientMixin._mean_absolute_error_derivative - ) - elif metric in ["mean squared error", "mse"]: - self._metric = mean_squared_error - self._metric_derivative = VectorGradientMixin._mean_squared_error_derivative - elif metric in ["root mean squared error", "rmse"]: - self._metric = root_mean_squared_error - self._metric_derivative = ( - VectorGradientMixin._root_mean_squared_error_derivative - ) - elif metric in ["negative nmll laplace"]: - self._metric = negative_nmll_laplace - self._metric_derivative = ( - VectorGradientMixin._negative_nmll_laplace_derivative - ) - elif metric in ["bic"]: - self._metric = bic - self._metric_derivative = VectorGradientMixin._bic_derivative - else: - raise ValueError("Invalid metric for vector gradient mixin") - - def get_fitness_and_gradient(self, individual): - """Fitness evaluation and gradient of vector based fitness - function using metric (i.e. the fitness function returns - a vector that is converted into a scalar using its metric function) - - Get the fitness of the individual and the gradient - of this function with respect to the individual's constants. - - Parameters - ---------- - individual : chromosomes - individual for which the fitness and gradient will be calculated for - - Returns - ------- - fitness, gradient : - fitness of the individual and the gradient of this function - with respect to the individual's constants - """ - fitness_vector, jacobian = self.get_fitness_vector_and_jacobian(individual) - return self._metric(fitness_vector, individual), self._metric_derivative( - fitness_vector, jacobian.transpose() - ) - - @abstractmethod - def get_fitness_vector_and_jacobian(self, individual): - r"""Returns the vectorized fitness of this individual and - the jacobian of this vector fitness function with - respect to the individual's constants - - jacobian = [[:math:`df_1/dc_1`, :math:`df_1/dc_2`, ...], - [:math:`df_2/dc_1`, :math:`df_2/dc_2`, ...], - ...] - where :math:`f_\#` is the fitness function corresponding with the - #th fitness vector entry and :math:`c_\#` is the corresponding - constant of the individual - - Parameters - ---------- - individual : chromosomes - individual used for vectorized fitness evaluation and jacobian - calculation - - Returns - ------- - fitness_vector, jacobian : - the vectorized fitness of the individual and - the partial derivatives of each fitness function with respect - to the individual's constants - """ - raise NotImplementedError - - @staticmethod - def _mean_absolute_error_derivative(fitness_vector, fitness_partials): - return np.mean(np.sign(fitness_vector) * fitness_partials, axis=1) - - @staticmethod - def _mean_squared_error_derivative(fitness_vector, fitness_partials): - return 2 * np.mean(fitness_vector * fitness_partials, axis=1) - - @staticmethod - def _root_mean_squared_error_derivative(fitness_vector, fitness_partials): - return ( - 1 - / np.sqrt(np.mean(np.square(fitness_vector))) - * np.mean(fitness_vector * fitness_partials, axis=1) - ) - - @staticmethod - def _negative_nmll_laplace_derivative(fitness_vector, fitness_partials): - n = len(fitness_vector) - b = 1 / np.sqrt(n) - dmse = 2 * np.mean(fitness_vector * fitness_partials, axis=1) - mse = np.mean(np.square(fitness_vector)) - dll = -0.5 * n / mse * dmse - dnmll = (1 - b) * dll - return -dnmll - - @staticmethod - def _bic_derivative(fitness_vector, fitness_partials): - n = len(fitness_vector) - mse = np.mean(np.square(fitness_vector)) - dmse = 2 * np.mean(fitness_vector * fitness_partials, axis=1) - return n / mse * dmse diff --git a/bingo/expressions/agraph/cppagraph/bindings/bind_expression.cpp b/bingo/expressions/agraph/cppagraph/bindings/bind_expression.cpp index 38ded362..edc7c3d6 100644 --- a/bingo/expressions/agraph/cppagraph/bindings/bind_expression.cpp +++ b/bingo/expressions/agraph/cppagraph/bindings/bind_expression.cpp @@ -447,7 +447,22 @@ void bind_expression(py::module_& m) { py::arg("required_params") = py::none()) .def_property_readonly("is_fitted", - &AGraphExpression::is_fitted) + &AGraphExpression::is_fitted) + + .def("commit_fit", + [](AGraphExpression& self, py::object constants) -> AGraphExpression& { + self.commit_fit(iterable_to_double_vec(constants)); + return self; + }, + py::arg("constants"), + py::return_value_policy::reference_internal) + + .def("clear_fit", + [](AGraphExpression& self) -> AGraphExpression& { + self.clear_fit(); + return self; + }, + py::return_value_policy::reference_internal) .def("__sklearn_is_fitted__", &AGraphExpression::is_fitted) diff --git a/bingo/expressions/agraph/cppagraph/include/cppagraph/expression.h b/bingo/expressions/agraph/cppagraph/include/cppagraph/expression.h index 5332d342..e3067562 100644 --- a/bingo/expressions/agraph/cppagraph/include/cppagraph/expression.h +++ b/bingo/expressions/agraph/cppagraph/include/cppagraph/expression.h @@ -227,6 +227,12 @@ class AGraphExpression { */ bool is_fitted(); + /** Validate and atomically install fitted constants. */ + void commit_fit(std::vector constants); + + /** Clear fittedness without changing constants or raw structure. */ + void clear_fit(); + /** Lifecycle state access for serialization and raw-state reconstruction. */ bool fit_attempted() const; void set_fit_attempted(bool v); diff --git a/bingo/expressions/agraph/cppagraph/src/expression.cpp b/bingo/expressions/agraph/cppagraph/src/expression.cpp index bf45b599..a3e1c822 100644 --- a/bingo/expressions/agraph/cppagraph/src/expression.cpp +++ b/bingo/expressions/agraph/cppagraph/src/expression.cpp @@ -384,10 +384,6 @@ void AGraphExpression::fit(const RowMatrixXd& X, double tolerance, int max_iter) { validate_explicit_data(X, y); if (modified_) update(); - // A fitting attempt establishes the fitted state for the current raw - // structure, even when the solver does not numerically converge. - fit_attempted_ = true; - if (constants_.empty()) return; const Eigen::Index m = X.rows(); @@ -459,10 +455,11 @@ void AGraphExpression::fit(const RowMatrixXd& X, } } } catch (...) { - // Don't crash on bad fits — keep current params. + // Preserve the prior lifecycle when the solver cannot produce a result. + return; } - set_constants(std::vector( + commit_fit(std::vector( params.data(), params.data() + params.size())); } @@ -471,11 +468,9 @@ void AGraphExpression::fit_implicit(const RowMatrixXd& X, double tolerance, int max_iter) { validate_implicit_data(X, dx_dt); if (modified_) update(); - // A fitting attempt establishes the fitted state regardless of convergence. - fit_attempted_ = true; - if (constants_.empty()) return; + const auto original_constants = constants_; const Eigen::Index n = static_cast(constants_.size()); @@ -555,10 +550,12 @@ void AGraphExpression::fit_implicit(const RowMatrixXd& X, } } } catch (...) { - // Don't crash on bad fits — keep current params. + set_constants(original_constants); + return; } - set_constants(std::vector( + set_constants(original_constants); + commit_fit(std::vector( params.data(), params.data() + params.size())); } @@ -686,6 +683,25 @@ bool AGraphExpression::is_fitted() { return fit_attempted_ || constants_.empty(); } +void AGraphExpression::commit_fit(std::vector constants) { + if (modified_) update(); + if (constants.size() != constants_.size()) { + throw std::invalid_argument( + "constants must have one entry per simplified expression constant"); + } + if (!std::all_of(constants.begin(), constants.end(), + [](double value) { return std::isfinite(value); })) { + throw std::invalid_argument("fitted constants must be finite"); + } + set_constants(std::move(constants)); + fit_attempted_ = true; +} + +void AGraphExpression::clear_fit() { + if (modified_) update(); + if (!constants_.empty()) fit_attempted_ = false; +} + bool AGraphExpression::fit_attempted() const { return fit_attempted_; } diff --git a/bingo/expressions/agraph/pyagraph/expression.py b/bingo/expressions/agraph/pyagraph/expression.py index a5a2a339..08d7c9ad 100644 --- a/bingo/expressions/agraph/pyagraph/expression.py +++ b/bingo/expressions/agraph/pyagraph/expression.py @@ -611,8 +611,6 @@ def fit(self, X, y, *, tolerance=1e-5): y = np.asarray(y, dtype=float).ravel() if X.shape[0] != y.size: raise ValueError("X and y must have the same number of samples") - self._fit_attempted = True - if len(self.constants) == 0: return self @@ -620,12 +618,10 @@ def fit(self, X, y, *, tolerance=1e-5): cached = CachedEvaluator(self._command_array, X, self._integers) def residuals(params): - self.constants = params - return cached.forward_eval(self._constants).ravel() - y + return cached.forward_eval(tuple(params)).ravel() - y def jacobian(params): - self.constants = params - _, jac = cached.forward_eval_with_const_derivative(self._constants) + _, jac = cached.forward_eval_with_const_derivative(tuple(params)) return jac try: @@ -636,7 +632,7 @@ def jacobian(params): method="lm", tol=tolerance, ) - self.constants = result.x + self.commit_fit(result.x) except Exception: # noqa: broad-except — don't crash on bad fits pass @@ -671,12 +667,11 @@ def fit_implicit(self, X, dx_dt, *, tolerance=1e-5): dx_dt = np.atleast_2d(np.asarray(dx_dt, dtype=float)) if X.shape != dx_dt.shape: raise ValueError("X and dx_dt must have the same shape") - self._fit_attempted = True - if len(self.constants) == 0: return self x0 = np.array(self.constants, dtype=float) + original_constants = self.constants def residuals(params): self.constants = params @@ -689,9 +684,10 @@ def residuals(params): ftol=tolerance, xtol=tolerance, ) - self.constants = result.x + self.constants = original_constants + self.commit_fit(result.x) except Exception: # noqa: broad-except — don't crash on bad fits - pass + self.constants = original_constants return self @@ -883,6 +879,29 @@ def _set_fit_attempted(self, value): """Restore the fitted lifecycle after constructing from raw state.""" self._fit_attempted = bool(value) + def commit_fit(self, constants): + """Atomically install validated fitted constants for this structure.""" + if self._modified: + self._update() + constants = tuple(float(value) for value in constants) + if len(constants) != len(self._constants): + raise ValueError( + "constants must have one entry per simplified expression constant" + ) + if not np.all(np.isfinite(constants)): + raise ValueError("fitted constants must be finite") + self.constants = constants + self._fit_attempted = True + return self + + def clear_fit(self): + """Clear fittedness without changing constants or raw structure.""" + if self._modified: + self._update() + if self._constants: + self._fit_attempted = False + return self + # ------------------------------------------------------------------ # # Simplification / utility # # ------------------------------------------------------------------ # diff --git a/bingo/local_optimizers/__init__.py b/bingo/local_optimizers/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/bingo/local_optimizers/local_opt_fitness.py b/bingo/local_optimizers/local_opt_fitness.py deleted file mode 100644 index 0e7b2963..00000000 --- a/bingo/local_optimizers/local_opt_fitness.py +++ /dev/null @@ -1,74 +0,0 @@ -"""Fitness evaluation with local optimization - -This module contains the implementation of a fitness function wrapper -that will perform local optimization of a `Chromosome` as necessary -using a `LocalOptimizer` before evaluating it. -""" - -from ..evaluation.fitness_function import FitnessFunction - - -class LocalOptFitnessFunction(FitnessFunction): - """Fitness function wrapper for individuals that want local optimization - - A class for fitness evaluation of individuals that may or may - not need local optimization before evaluation. - - Parameters - ---------- - fitness_function : `FitnessFunction` - A `FitnessFunction` for evaluating the fitness of a `Chromosome`. - optimizer : `LocalOptimizer` - An optimizer that will perform local optimization on a - `Chromosome` before evaluation as needed. - - Attributes - ---------- - eval_count : int - the number of evaluations that have been performed by the wrapped - fitness function - training_data : `TrainingData` - data that can be used in the wrapped fitness function - """ - - def __init__(self, fitness_function, optimizer): - # pylint: disable=super-init-not-called - self._fitness_function = fitness_function - self.optimizer = optimizer - - @property - def training_data(self): - """TrainingData : data that can be used in fitness evaluations""" - return self._fitness_function.training_data - - @training_data.setter - def training_data(self, value): - self._fitness_function.training_data = value - - @property - def eval_count(self): - """int : the number of evaluations that have been performed""" - return self._fitness_function.eval_count - - @eval_count.setter - def eval_count(self, value): - self._fitness_function.eval_count = value - - def __call__(self, individual): - """Evaluates the fitness of the individual. Provides local optimization - on the individual if necessary. - - Parameters - ---------- - individual : `Chromosome` - Individual to calculate the fitness of. Local optimization - is performed if necessary before evaluation. - - Returns - ------- - float - The fitness of the individual - """ - if individual.needs_local_optimization(): - self.optimizer(individual) - return self._fitness_function(individual) diff --git a/bingo/local_optimizers/local_optimizer.py b/bingo/local_optimizers/local_optimizer.py deleted file mode 100644 index b3675fe0..00000000 --- a/bingo/local_optimizers/local_optimizer.py +++ /dev/null @@ -1,44 +0,0 @@ -"""This module contains the abstract definition of an optimizer -that can be used for local optimization of a `Chromosome`. -""" - -from abc import ABCMeta, abstractmethod - - -class LocalOptimizer(metaclass=ABCMeta): - """An abstract base class for optimizing a `Chromosome`. - """ - @property - @abstractmethod - def objective_fn(self): - """function to minimize, must take a `Chromosome` as input - and return a number""" - raise NotImplementedError - - @objective_fn.setter - @abstractmethod - def objective_fn(self, value): - raise NotImplementedError - - @property - @abstractmethod - def options(self): - """dict : optimizer's options""" - raise NotImplementedError - - @options.setter - @abstractmethod - def options(self, value): - raise NotImplementedError - - @abstractmethod - def __call__(self, individual): - """Performs local optimization of the individual - based on minimizing this object's objective_fn. - - Parameters - ---------- - individual : `Chromosome` - The individual who will be optimized. - """ - raise NotImplementedError diff --git a/bingo/local_optimizers/normalized_marginal_likelihood.py b/bingo/local_optimizers/normalized_marginal_likelihood.py deleted file mode 100644 index 8c5a1a7e..00000000 --- a/bingo/local_optimizers/normalized_marginal_likelihood.py +++ /dev/null @@ -1,82 +0,0 @@ -"""Normalized marginal likelihood calculation using SMCPy - -This module contains the implementation of a fitness function wrapper -that will perform probabilistic local optimization of a `Chromosome` using -SMCPy. The normaized marginal likelihood from the SMC optimization is returned. -""" - -import numpy as np -from bingo.local_optimizers.smcpy_optimizer import SmcpyOptimizer -from ..evaluation.fitness_function import FitnessFunction - - -class NormalizedMarginalLikelihood(FitnessFunction): - """Normalized marginal likelihood calculation using SMCPy - - A class for fitness evaluation of individuals that have local optimization - parameters - - Parameters - ---------- - fitness_function : `FitnessFunction` - A `FitnessFunction` for evaluating the fitness of a `Chromosome`. - deterministic_optimizer : `LocalOptimizer` - An optimizer that will perform deterministic local optimization on a - `Chromosome`. Used in proposals of `SmcpyOptimizer` - **kwargs: - other keyword arguments are passed to the SmcpyOptimizer initialization - - Attributes - ---------- - eval_count : int - the number of evaluations that have been performed by the wrapped - fitness function - training_data : `TrainingData` - data that can be used in the wrapped fitness function - """ - - def __init__( - self, fitness_function, deterministic_optimizer, log_scale=True, **kwargs - ): - # pylint: disable=super-init-not-called - self._log_scale = log_scale - self.optimizer = SmcpyOptimizer( - fitness_function, deterministic_optimizer, **kwargs - ) - - @property - def training_data(self): - """TrainingData : data that can be used in fitness evaluations""" - return self.optimizer.training_data - - @training_data.setter - def training_data(self, value): - self.optimizer.training_data = value - - @property - def eval_count(self): - """int : the number of evaluations that have been performed""" - return self.optimizer.eval_count - - @eval_count.setter - def eval_count(self, value): - self.optimizer.eval_count = value - - def __call__(self, individual): - """Evaluates the normalized marginal likelihood of the individual. - - Parameters - ---------- - individual : `Chromosome` - Individual to calculate the normalized marginal likelihood of. - Probabilistic local optimization is performed during evaluation. - - Returns - ------- - float - The *negative* normalized marginal likelihood - """ - log_nml = self.optimizer(individual)[0] - if self._log_scale: - return -log_nml - return -np.exp(log_nml) diff --git a/bingo/local_optimizers/scipy_optimizer.py b/bingo/local_optimizers/scipy_optimizer.py deleted file mode 100644 index aab473e6..00000000 --- a/bingo/local_optimizers/scipy_optimizer.py +++ /dev/null @@ -1,239 +0,0 @@ -"""Local optimization using scipy - -Specifies `ScipyOptimizer` which is a class for local optimization of -`Chromosome`s using scipy's minimize or root methods. Also -specifies ROOT_SET, a set of methods that will use scipy's root method; -MINIMIZE_SET, a set of methods that will use scipy's minimize method; -and JACOBIAN_SET, a set of methods that will use jacobian information. -""" - -import numpy as np -from scipy import optimize - -from .local_optimizer import LocalOptimizer -from ..evaluation.gradient_mixin import GradientMixin, VectorGradientMixin - -ROOT_SET = { - # "hybr", - "lm" - # "broyden1", - # "broyden2", - # "anderson", - # "linearmixing", - # "diagbroyden", - # "excitingmixing", - # "krylov", - # "df-sane" -} - -MINIMIZE_SET = { - "Nelder-Mead", - "Powell", - "CG", - "BFGS", - # "Newton-CG", - "L-BFGS-B", - "TNC", - # "COBYLA", - "SLSQP" - # "trust-constr" - # "dogleg", - # "trust-ncg", - # "trust-exact", - # "trust-krylov" -} - -JACOBIAN_SET = { - "CG", - "BFGS", - # "Newton-CG", - "L-BFGS-B", - "TNC", - "SLSQP", - # "trust-constr" - # "dogleg", - # "trust-ncg", - # "trust-exact", - # "trust-krylov", - # "hybr", - "lm" -} - - -class ScipyOptimizer(LocalOptimizer): - """An optimizer that uses scipy.minimize or scipy.root - for local optimization - - A class for optimizing the parameters of a `Chromosome` using - either scipy.minimize or scipy.root depending on the method - specified. - - Parameters - ---------- - objective_fn - A function to minimize which can be evaluated by passing in a - `Chromosome`. - options - Additional arguments for optimization. - e.g. (..., tol=1e-8, options={"maxiter": 1000}) - - e.g. param_init_bounds: iterable - [low, high) bounds that are used to initialize params, - formatted as an iterable - defaults to [-10000, 10000) - - e.g. method: string - method to use for optimization (e.g. BFGS, lm, etc.) - defaults to BFGS - - e.g. tol: float - tolerance used for method - defaults to 1e-6 - - e.g. options : dict - method-specific options (e.g. maxiter, ftol, etc.) - - Attributes - ---------- - objective_fn - A function to minimize which can be evaluated by passing in a - `Chromosome` - options : dict - Additional arguments for clo options - - Raises - ------ - KeyError - `method` must be a method supported by scipy - TypeError - `objective_function` must suit the specified method - """ - def __init__(self, objective_fn, **options): - self.options = options - self.objective_fn = objective_fn - - @property - def objective_fn(self): - """function to minimize, must take a `Chromosome` as input - and return a number""" - return self._objective_fn - - @objective_fn.setter - def objective_fn(self, obj_fn): - self._jacobian_capable = isinstance(obj_fn, VectorGradientMixin) - self._gradient_capable = isinstance(obj_fn, GradientMixin) - self._objective_fn = obj_fn - self._verify_objective_fn(obj_fn, self.options["method"]) - - @staticmethod - def _verify_objective_fn(objective_fn, method): - if method in ROOT_SET and not hasattr(objective_fn, - "evaluate_fitness_vector"): - raise TypeError(f"{method} requires VectorBasedFunction \ - as a fitness function") - - @property - def options(self): - """dict : optimizer options (e.g. param_init_bounds, method, - tol, options, etc.)""" - return self._options - - @options.setter - def options(self, kwargs): - self._options = kwargs - - # set default param init bounds to [-10000, 10000) if not included - if "param_init_bounds" not in self._options.keys(): - self._options["param_init_bounds"] = [-10000, 10000] - - # set default method to BFGS if not included - if "method" not in self._options.keys(): - self._options["method"] = "BFGS" - self._verify_method(self._options["method"]) - - # set default tol to 1e-6 if not included - if "tol" not in self._options.keys(): - self._options["tol"] = 1e-6 - - # scipy_options = normal options w/o param_init_bounds - self._scipy_options = {k: v for k, v in self._options.items() if - k != "param_init_bounds"} - - @staticmethod - def _verify_method(method): - if method not in ROOT_SET and method not in MINIMIZE_SET: - raise KeyError(f"{method} is not a listed method") - - def __call__(self, individual): - """Performs local optimization of the individual - based on minimizing this object's objective_fn. - - Parameters - ---------- - individual : `Chromosome` - The individual who will be optimized. - """ - num_params = individual.get_number_local_optimization_params() - c_0 = np.random.uniform(*self.options["param_init_bounds"], num_params) - params = self._run_method_for_optimization( - self._sub_routine_for_obj_fn, individual, c_0) - individual.set_local_optimization_params(params) - - def _sub_routine_for_obj_fn(self, params, individual): - individual.set_local_optimization_params(params) - - if self.options["method"] in ROOT_SET: - return self.objective_fn.evaluate_fitness_vector(individual) - return self.objective_fn(individual) - - def _run_method_for_optimization(self, sub_routine, individual, params): - backend, jacobian = self._get_scipy_backend_and_jacobian_fn() - try: - optimize_result = backend( - sub_routine, - params, - args=individual, - jac=jacobian, - **self._scipy_options - ) - return optimize_result.x - except TypeError: # issue with too many constants using root method - old_method = self.options["method"] - - self.options["method"] = "BFGS" # use minimize method instead - self._scipy_options["method"] = "BFGS" - - backend, jacobian = self._get_scipy_backend_and_jacobian_fn() - optimize_result = backend( - sub_routine, - params, - args=individual, - jac=jacobian, - **self._scipy_options - ) - self.options["method"] = old_method # reset to old method - self._scipy_options["method"] = old_method - return optimize_result.x - - def _get_scipy_backend_and_jacobian_fn(self): - def get_just_jacobian(_, indv): - return self.objective_fn.get_fitness_vector_and_jacobian(indv)[1] - - def get_just_gradient(_, indv): - return self.objective_fn.get_fitness_and_gradient(indv)[1] - - backend = optimize.minimize - jacobian = False - - jacobian_method = self.options["method"] in JACOBIAN_SET - - if self.options["method"] in ROOT_SET: - backend = optimize.root - if jacobian_method and self._jacobian_capable: - jacobian = get_just_jacobian - - else: # MINIMIZE_SET - if jacobian_method and self._gradient_capable: - jacobian = get_just_gradient - - return backend, jacobian diff --git a/bingo/local_optimizers/smcpy_optimizer.py b/bingo/local_optimizers/smcpy_optimizer.py deleted file mode 100644 index 10e21db5..00000000 --- a/bingo/local_optimizers/smcpy_optimizer.py +++ /dev/null @@ -1,513 +0,0 @@ -"""A module for probabilistic calibration of parameters. - -Probabilistic calibration of model parameters can be useful in cases where data -is sparse and/or noisy. Using a calibration of this type can allow for a better -estimate of true fitness while being a bit more robust to overfitting. -""" - -import numpy as np -from scipy.stats import multivariate_normal as mvn -from scipy.stats import invgamma - -from smcpy import VectorMCMC, VectorMCMCKernel, AdaptiveSampler, ImproperUniform -from smcpy.paths import GeometricPath -from smcpy.proposals import MultivarIndependent - -from .local_optimizer import LocalOptimizer - - -class SmcpyOptimizer(LocalOptimizer): - """An optimizer that uses SMCPy for probabilistic parameter calibration - - A class for probabilistic parameter calibration for the parameters of a - `Chromosome` using SMCPy - - Parameters - ---------- - objective_fn : VectorBasedFunction, VectorGradientMixin - A `VectorBasedFunction` with `VectorGradientMixin` (e.g., - ExplicitRegression). It should produce a vector where the target value - is 0. - deterministic_optimizer : LocalOptimizer - A deterministic local optimizer e.g., `ScipyOptimizer` - num_particles : int - The number of particles to use in the SMC approximation - mcmc_steps : int - The number of MCMC steps to perform with each SMC update - ess_threshold : float (0-1) - The effective sample size (ratio) below which SMC particles will be - resampled - std : float - (Optional) The fixed noise level, if it is known - num_multistarts : int - (Optional) The number of deterministic optimizations performed when - developing the SMC proposal - - Attributes - ---------- - objective_fn - A function to minimize which can be evaluated by passing in a - `Chromosome` - options : dict - Additional arguments for clo options - - """ - - def __init__( - self, - objective_fn, - deterministic_optimizer, - num_particles=150, - mcmc_steps=12, - ess_threshold=0.75, - std=None, - num_multistarts=1, - reuse_starting_point=True, - ): - - self._num_particles = num_particles - self._mcmc_steps = mcmc_steps - self._ess_threshold = ess_threshold - self._std = std - self._num_multistarts = num_multistarts - self._objective_fn = objective_fn - self._deterministic_optimizer = deterministic_optimizer - self._reuse_starting_point = reuse_starting_point - - self._norm_phi = self._calculate_norm_phi() - - def _calculate_norm_phi(self): - num_observations = len(self.training_data) - return 1 / np.sqrt(num_observations) - - @property - def objective_fn(self): - """A `VectorBasedFunction` with `VectorGradientMixin` (e.g., - ExplicitRegression). It should produce a vector where the target value - is 0.""" - return self._objective_fn - - @objective_fn.setter - def objective_fn(self, value): - self._objective_fn = value - - @property - def training_data(self): - """Training data used in objective function""" - return self._objective_fn.training_data - - @training_data.setter - def training_data(self, value): - self._objective_fn.training_data = value - self._norm_phi = self._calculate_norm_phi() - - @property - def eval_count(self): - """int : the number of evaluations that have been performed""" - return self._objective_fn.eval_count - - @eval_count.setter - def eval_count(self, value): - self._objective_fn.eval_count = value - - @property - def options(self): - """dict : optimizer's options""" - return { - "num_particles": self._num_particles, - "mcmc_steps": self._mcmc_steps, - "ess_threshold": self._ess_threshold, - "std": self._std, - "num_multistarts": self._num_multistarts, - } - - @options.setter - def options(self, value): - if "num_particles" in value: - self._num_particles = value["num_particles"] - if "mcmc_steps" in value: - self._mcmc_steps = value["mcmc_steps"] - if "ess_threshold" in value: - self._ess_threshold = value["ess_threshold"] - if "std" in value: - self._std = value["std"] - if "num_multistarts" in value: - self._num_multistarts = value["num_multistarts"] - - def __call__(self, individual): - try: - proposal = self._generate_proposals(individual) - except (ValueError, np.linalg.LinAlgError, RuntimeError) as e: - return np.nan, "proposal error", e - - param_names = self._get_parameter_names(individual) - priors = [ImproperUniform() for _ in range(len(param_names))] - if self._std is None: - priors.append(ImproperUniform(0, None)) - param_names.append("std_dev") - - path = GeometricPath(proposal=proposal, required_phi=self._norm_phi) - vmcmc = VectorMCMC( - lambda x: self.evaluate_model(x, individual), - np.zeros(len(self.training_data)), - priors, - log_like_args=self._std, - ) - kernel = VectorMCMCKernel(vmcmc, param_order=param_names, path=path) - smc = AdaptiveSampler(kernel, show_progress_bar=False) - - try: - step_list, marginal_log_likes = smc.sample( - self._num_particles, - self._mcmc_steps, - self._ess_threshold, - ) - except (ValueError, np.linalg.LinAlgError, ZeroDivisionError) as e: - # print(e) - return np.nan, "sample error", e - - max_idx = np.argmax(step_list[-1].log_likes) - maps = step_list[-1].params[max_idx] - individual.set_local_optimization_params(maps[:-1]) - - norm_phi = 1 / np.sqrt(len(self.training_data)) - norm_phi_index = np.argmin(np.abs(np.array(smc._phi_sequence) - norm_phi)) - log_nml = marginal_log_likes[-1] - marginal_log_likes[norm_phi_index] - - return log_nml, step_list, vmcmc - - def _generate_proposals(self, individual): - param_names = self._get_parameter_names(individual) - num_multistarts = self._num_multistarts - - param_dists = [] - cov_estimates = [] - mix_dists = [] - if not param_names: - cov_estimates.append(self._estimate_covariance(individual)) - else: - for i in range(3 * num_multistarts): - try: - do_det_opt = not self._reuse_starting_point or i != 0 - mean, cov, var_ols, ssqe = self._estimate_covariance( - individual, do_det_opt - ) - cov = 0.5 * (cov + cov.T) # ensuring symmetry - evals, evecs = np.linalg.eig(cov) - # the below approximation attempts to correct for cov matrices - # that are not positive semidefinite - if np.min(evals) < 0: - diag = np.diag(evals) - diag[diag < 0] = 0 - cov = evecs.dot(diag).dot(evecs.T) - dists = mvn(mean, cov, allow_singular=True) - except (ValueError, np.linalg.LinAlgError) as _: - continue - cov_estimates.append((mean, cov, var_ols, ssqe)) - param_dists.append(dists) - if len(param_dists) == num_multistarts: - break - if not param_dists: - raise RuntimeError( - "Could not generate any valid proposal distributions" - ) - - mix_dists.append(MixtureDist(*param_dists)) - - len_data = len(self.training_data) - scale_data = np.sqrt( - np.mean(np.square(self.training_data.y)) - ) # TODO can we do this differently without knowing what the training data is? - noise_dists = [] - for _, _, var_ols, ssqe in cov_estimates: - shape = (0.01 + len_data) / 2 - scale = max((0.01 * var_ols + ssqe) / 2, 1e-12 * scale_data) - noise_dists.append(SqrtInvGamma(shape, scale=scale)) - param_names.append("std_dev") - - mix_dists.append(MixtureDist(*noise_dists)) - - return MultivarIndependent(*mix_dists) - - @staticmethod - def _get_parameter_names(individual): - num_params = individual.get_number_local_optimization_params() - return [f"p{i}" for i in range(num_params)] - - def _estimate_covariance(self, individual, do_det_opt=True): - if do_det_opt: - self._deterministic_optimizer(individual) - - # # RALPH data approx method - f, f_deriv = self._objective_fn.get_fitness_vector_and_jacobian(individual) - ssqe = np.sum((f) ** 2) - var_ols = ssqe / len(f) - cov = var_ols * np.linalg.inv(f_deriv.T.dot(f_deriv)) - - # LAPLACE approx - # f, g = self._objective_fn.get_fitness_vector_and_jacobian( - # individual - # ) - # h = np.squeeze( - # individual.evaluate_with_local_opt_hessian_at( - # self.objective_fn.training_data.x - # )[1].detach().numpy(), - # 1 - # ) - # A = 2*np.sum(np.einsum('...i,...j->...ij', g, g) - # + np.expand_dims(f, axis=(1,2))*h, axis=0) - # ssqe = np.sum((f) ** 2) - # var_ols = ssqe / len(f) - # # try: - # cov = np.linalg.inv(A) - # # except np.linalg.LinAlgError: - # # # print(A) - # # # A = A+np.empty_like(A)*1e-8 # adding nugget for invertability - # # # print(A) - # # # cov = np.linalg.inv(A) - # # # print(cov) - # # cov = np.linalg.pinv(A) - # # # print(cov) - - return individual.constants, cov, var_ols, ssqe - - def evaluate_model(self, params, individual): - """ - Evaluate a model with given parameters and return fitness vector. - - This method sets the local optimization parameters for an individual, - evaluates its fitness using the objective function, and reshapes the - result to ensure consistent dimensionality for further processing. - - Parameters - ---------- - params : ndarray - Model parameters to evaluate. Expected shape is (n_params,) or - (n_params, n_models). Will be transposed before setting on individual. - individual : object - Individual model object that implements `set_local_optimization_params` - method. Represents the model structure or configuration to evaluate. - - Returns - ------- - ndarray - Fitness evaluation results with shape (n_models, n_training_samples) - where n_training_samples is the length of the training data. If the - original result is 1D, it will be reshaped to ensure 2D output. - - Examples - -------- - >>> # Assuming self is an instance with _objective_fn and training data - >>> params = np.array([[1.0, 2.0], [3.0, 4.0]]) # 2 parameters, 2 models - >>> result = self.evaluate_model(params, individual) - >>> result.shape - (2, 100) # 2 models evaluated on 100 training samples - """ - individual.set_local_optimization_params(params.T) - result = self._objective_fn.evaluate_fitness_vector(individual).T - if len(result.shape) < 2: - # TODO, would it be better to remove the flatten in explicit - # regression and add a flatten to the scipy wrapper? - result = result.reshape(-1, len(self._objective_fn.training_data)) - return result - - -class MixtureDist: - """ - A mixture distribution class that combines multiple probability distributions. - - This class represents a mixture model where samples are drawn uniformly at random - from one of the component distributions. Each component distribution has equal - weight (1/n where n is the number of components). - - Parameters - ---------- - *args : tuple of distribution objects - Variable number of probability distribution objects. Each distribution - should have `rvs()` and `pdf()` methods compatible with scipy.stats - distributions. - - Attributes - ---------- - _dists : tuple - Tuple of component probability distributions. - - Examples - -------- - >>> from scipy.stats import norm, uniform - >>> mixture = MixtureDist(norm(0, 1), uniform(-2, 4)) - >>> samples = mixture.rvs(1000, random_state=42) - >>> log_probs = mixture.logpdf(samples) - """ - - def __init__(self, *args): - self._dists = args - - def rvs(self, num_samples, random_state=None): - """ - Generate random samples from the mixture distribution. - - For each sample, randomly selects one of the component distributions with - equal probability and draws a sample from it. - - Parameters - ---------- - num_samples : int - Number of samples to generate. - random_state : int, optional - Random seed for reproducible results. If None, uses a random seed. - - Returns - ------- - ndarray - Array of shape (num_samples, dim) where dim is the dimensionality - of the component distributions. For 1D distributions, returns - shape (num_samples, 1). - - Examples - -------- - >>> mixture = MixtureDist(norm(0, 1), norm(5, 2)) - >>> samples = mixture.rvs(100, random_state=42) - >>> samples.shape - (100, 1) - """ - candidate_samples = np.zeros( - ( - len(self._dists), - num_samples, - self._dists[0].dim if hasattr(self._dists[0], "dim") else 1, - ) - ) - - for i, d in enumerate(self._dists): - candidate_samples[i, :, :] = d.rvs( - num_samples, random_state=random_state - ).reshape(num_samples, -1) - - if random_state is None: - rng = np.random.default_rng(seed=np.random.randint(np.iinfo(np.int16).max)) - else: - rng = np.random.default_rng(random_state) - dist_indices = rng.integers(0, len(self._dists), num_samples) - sample_indices = np.arange(0, num_samples) - - return candidate_samples[dist_indices, sample_indices, :] - - def logpdf(self, x): - """ - Compute the log probability density function of the mixture distribution. - - The PDF of a mixture distribution is the average of the PDFs of the - component distributions: pdf(x) = (1/n) * sum(pdf_i(x)) where n is the - number of components. - - Parameters - ---------- - x : ndarray - Input samples of shape (num_samples, dim) where dim matches the - dimensionality of component distributions. - - Returns - ------- - ndarray - Log probability densities of shape (num_samples, 1). - - Examples - -------- - >>> import numpy as np - >>> mixture = MixtureDist(norm(0, 1), norm(5, 2)) - >>> x = np.array([[0.5], [2.0], [4.5]]) - >>> log_probs = mixture.logpdf(x) - """ - num_samples = x.shape[0] - pdfs = np.zeros((len(self._dists), num_samples, 1)) - - for i, d in enumerate(self._dists): - pdfs[i, :, :] = d.pdf(x).reshape(num_samples, 1) - - return np.log(pdfs.sum(axis=0) / len(self._dists)) - - -class SqrtInvGamma: - """ - Square root of inverse gamma distribution. - - This class represents a distribution where if X ~ InvGamma(shape, scale), - then Y = sqrt(X) follows this distribution. This is useful when you need - the square root transformation of an inverse gamma random variable. - - Parameters - ---------- - shape : float - Shape parameter of the underlying inverse gamma distribution. - Must be positive. - scale : float - Scale parameter of the underlying inverse gamma distribution. - Must be positive. - - Examples - -------- - >>> sqrt_inv_gamma = SqrtInvGamma(shape=2.0, scale=1.0) - >>> samples = sqrt_inv_gamma.rvs(1000, random_state=42) - >>> probabilities = sqrt_inv_gamma.pdf(samples) - """ - - def __init__(self, shape, scale): - self._dist = invgamma(shape, scale=scale) - - def rvs(self, *args, **kwargs): - """ - Generate random samples from the square root inverse gamma distribution. - - Parameters - ---------- - *args : tuple - Positional arguments passed to the underlying inverse gamma - distribution's rvs method (e.g., size, random_state). - **kwargs : dict - Keyword arguments passed to the underlying inverse gamma - distribution's rvs method. - - Returns - ------- - ndarray or float - Square root of inverse gamma random samples. Shape depends on - the size parameter passed. - - Examples - -------- - >>> dist = SqrtInvGamma(shape=2.0, scale=1.0) - >>> samples = dist.rvs(size=100, random_state=42) - """ - return np.sqrt(self._dist.rvs(*args, **kwargs)) - - def pdf(self, x, *args, **kwargs): - """ - Compute the probability density function of the square root inverse gamma distribution. - - Uses the transformation formula: if Y = sqrt(X) where X ~ InvGamma(shape, scale), - then pdf_Y(y) = pdf_X(y^2) where pdf_X is the inverse gamma PDF. - - Parameters - ---------- - x : ndarray or float - Points at which to evaluate the PDF. Must be non-negative. - *args : tuple - Additional positional arguments passed to the underlying - inverse gamma PDF method. - **kwargs : dict - Additional keyword arguments passed to the underlying - inverse gamma PDF method. - - Returns - ------- - ndarray or float - Probability density values at the input points. - - Examples - -------- - >>> dist = SqrtInvGamma(shape=2.0, scale=1.0) - >>> x = np.linspace(0.1, 3.0, 100) - >>> pdf_values = dist.pdf(x) - """ - return self._dist.pdf(np.square(x), *args, **kwargs) diff --git a/bingo/selection/probabilistic_tournament.py b/bingo/selection/probabilistic_tournament.py index dc9df16b..07741591 100644 --- a/bingo/selection/probabilistic_tournament.py +++ b/bingo/selection/probabilistic_tournament.py @@ -14,10 +14,8 @@ class ProbabilisticTournament(Selection): """Tournament selection using probabilistic model selection - Individuals are chosen with a probability equal to the relative vale of - their fitness. When used in conjunction with `NormalizedMarginalLikelihood` - this results in selection with Bayesian Model Selection (Based on the - Fractional Bayes Factor) + Individuals are chosen with a probability equal to the relative value of + their fitness. Parameters ---------- diff --git a/bingo/symbolic_regression/__init__.py b/bingo/symbolic_regression/__init__.py index 2a1e7ff5..898c5d14 100644 --- a/bingo/symbolic_regression/__init__.py +++ b/bingo/symbolic_regression/__init__.py @@ -1,7 +1,36 @@ """Public symbolic-regression objectives and scikit-learn estimator.""" +from .custom_regression import CustomRegression +from .evidence import EvidenceResult, SmcEvidenceEstimator from .explicit_regression import ExplicitRegression +from .fitting import ( + FitResult, + ResidualMeasure, + ScalarMeasure, + ScipyFitter, + explicit_residuals, + expression_loss, + implicit_loss, + implicit_residuals, +) from .implicit_regression import ImplicitRegression +from .objective_data import ObjectiveData from .symbolic_regressor import SymbolicRegressor -__all__ = ["SymbolicRegressor", "ExplicitRegression", "ImplicitRegression"] +__all__ = [ + "SymbolicRegressor", + "CustomRegression", + "SmcEvidenceEstimator", + "EvidenceResult", + "ExplicitRegression", + "ImplicitRegression", + "ObjectiveData", + "ScipyFitter", + "FitResult", + "ResidualMeasure", + "ScalarMeasure", + "explicit_residuals", + "implicit_residuals", + "expression_loss", + "implicit_loss", +] diff --git a/bingo/symbolic_regression/custom_regression.py b/bingo/symbolic_regression/custom_regression.py new file mode 100644 index 00000000..033ededd --- /dev/null +++ b/bingo/symbolic_regression/custom_regression.py @@ -0,0 +1,77 @@ +"""Advanced independently configurable expression-regression objective.""" + +import numpy as np + +from ._expression_regression_objective import _ExpressionRegressionObjective +from .fitting import FitResult +from .objective_data import ObjectiveData + + +class CustomRegression(_ExpressionRegressionObjective): + """Fit Expressions with a selected fitter and rank them with a scalar loss. + + Fitter, measure, and loss callables may be any compatible Python callable + during serial evaluation. Multiprocessing evaluation and checkpointing + require those callables to be pickleable. + + Parameters + ---------- + data : ObjectiveData + Aligned arrays provided to the fitter, fitting measure, and loss. + fitter : callable + Called as ``fitter(expression, data, fitting_measure)`` and returns a + :class:`FitResult`. + fitting_measure : callable + Called as ``fitting_measure(expression, data, constants)`` by fitter. + loss : callable + Lower-is-better scalar called as ``loss(expression, data)``. + """ + + def __init__(self, data, fitter, fitting_measure, loss): + if not isinstance(data, ObjectiveData): + raise TypeError("CustomRegression data must be an ObjectiveData instance") + if not callable(fitter): + raise TypeError("CustomRegression fitter must be callable") + if not callable(fitting_measure): + raise TypeError("CustomRegression fitting_measure must be callable") + if not callable(loss): + raise TypeError("CustomRegression loss must be callable") + super().__init__(data) + self._fitter = fitter + self._fitting_measure = fitting_measure + self._loss = loss + + def __call__(self, individual): + """Fit and rank an individual, committing only after successful ranking.""" + expression = individual.expression + data = self._objective_data + if expression.is_fitted or not expression.constants: + loss = self._expression_loss(expression, data) + else: + constants = self._fit_constants(expression, data) + trial = expression.copy() + trial.commit_fit(constants) + loss = self._expression_loss(trial, data) + expression.commit_fit(constants) + self.eval_count += 1 + return loss + + def _fit_expression(self, expression, data): + expression.commit_fit(self._fit_constants(expression, data)) + + def _fit_constants(self, expression, data): + if not expression.constants: + return np.empty(0) + result = self._fitter(expression, data, self._fitting_measure) + if not isinstance(result, FitResult): + raise TypeError("fitters must return a FitResult") + constants = np.asarray(result.constants, dtype=float) + if constants.ndim != 1: + raise ValueError("FitResult constants must be a one-dimensional array") + return constants + + def _expression_loss(self, expression, data): + value = np.asarray(self._loss(expression, data), dtype=float) + if value.ndim != 0: + raise ValueError("CustomRegression loss must return a scalar") + return float(value) diff --git a/bingo/symbolic_regression/evidence.py b/bingo/symbolic_regression/evidence.py new file mode 100644 index 00000000..72382e4b --- /dev/null +++ b/bingo/symbolic_regression/evidence.py @@ -0,0 +1,433 @@ +"""Sequential-Monte-Carlo evidence estimation for Expressions.""" + +# The public estimator intentionally has several sampling controls. SMCPy must +# be imported lazily because it is an optional evidence-estimation dependency. +# pylint: disable=broad-exception-caught,import-outside-toplevel,protected-access,too-many-arguments,too-many-instance-attributes,too-many-locals,too-many-positional-arguments + +from dataclasses import dataclass, field +import hashlib + +import numpy as np +from scipy import optimize +from scipy.stats import invgamma, multivariate_normal + +from .fitting import FitResult, ScipyFitter + + +class _UserCodeError(Exception): + """Mark an exception raised by a user-provided fitter or measure.""" + + def __init__(self, error): + super().__init__(str(error)) + self.error = error + + +@dataclass(frozen=True) +class EvidenceResult: + """The outcome of an Evidence estimation attempt. + + ``smc_nmll`` is higher-is-better. Failed estimates have negative infinite + NMLL and do not contain a MAP estimate. + + Parameters + ---------- + smc_nmll : float + Higher-is-better normalized marginal log likelihood. + map_constants : tuple of float or None + Posterior MAP Expression constants when estimation succeeds. + success : bool + Whether sampling and Evidence calculation succeeded. + message : str, optional + Failure or status message. + diagnostics : dict, optional + Estimator diagnostics. + posterior : object, optional + Posterior samples when explicitly requested. + """ + + smc_nmll: float + map_constants: tuple[float, ...] | None + success: bool + message: str | None = None + diagnostics: dict = field(default_factory=dict) + posterior: object | None = None + + +def smc_nmll_loss(result): + """Adapt an Evidence result to Bingo's lower-is-better Loss convention. + + Parameters + ---------- + result : EvidenceResult + Evidence estimation outcome. + + Returns + ------- + float + The negated NMLL, or infinity when estimation failed. + """ + if not result.success or not np.isfinite(result.smc_nmll): + return np.inf + return -result.smc_nmll + + +class SmcEvidenceEstimator: + """Estimate Expression evidence with an SMC sampler and Laplace proposal. + + Parameters + ---------- + num_particles, mcmc_steps, ess_threshold : optional + SMCPy sampling controls. + std : float, optional + Known observation noise. When omitted it is inferred by SMC. + num_multistarts : int, optional + Number of Laplace proposal components. + seed : int, optional + Root seed. The expression's raw structure derives the sampler stream. + return_posterior : bool, optional + Include posterior steps in :class:`EvidenceResult` when true. + fitter : callable, optional + Fitter used only to construct the Laplace proposal. + """ + + def __init__( + self, + *, + num_particles=150, + mcmc_steps=12, + ess_threshold=0.75, + std=None, + num_multistarts=1, + seed=0, + return_posterior=False, + fitter=None, + ): + self.num_particles = num_particles + self.mcmc_steps = mcmc_steps + self.ess_threshold = ess_threshold + self.std = std + self.num_multistarts = num_multistarts + self.seed = seed + self.return_posterior = return_posterior + self.fitter = fitter if fitter is not None else ScipyFitter("least_squares") + + def estimate(self, expression, data, residual_measure): + """Estimate evidence, installing posterior MAP constants on success. + + Parameters + ---------- + expression : Expression + Expression whose constants are estimated. + data : ObjectiveData + Aligned data consumed by ``residual_measure``. + residual_measure : callable + Returns residuals for an Expression, data, and candidate constants. + + Returns + ------- + EvidenceResult + The Evidence estimate and optional posterior samples. + + Raises + ------ + Exception + Any exception raised by a user-supplied fitter or residual measure. + """ + try: + from smcpy import ( + AdaptiveSampler, + ImproperUniform, + VectorMCMC, + VectorMCMCKernel, + ) + from smcpy.paths import GeometricPath + from smcpy.proposals import MultivarIndependent + except ImportError as error: + return self._failure("SMCPy is required for evidence estimation", error) + try: + proposal, diagnostics = self._generate_proposal( + expression, data, residual_measure, MultivarIndependent + ) + except _UserCodeError as error: + raise error.error from error + except (TypeError, ValueError, np.linalg.LinAlgError, RuntimeError) as error: + return self._failure("proposal error", error) + + parameter_count = len(expression.constants) + priors = [ImproperUniform() for _ in range(parameter_count)] + parameter_names = [f"c{i}" for i in range(parameter_count)] + if self.std is None: + priors.append(ImproperUniform(0, None)) + parameter_names.append("std_dev") + + def evaluate_model(params): + params = np.asarray(params, dtype=float) + if params.ndim == 1: + params = params.reshape(1, -1) + constants = params[:, :parameter_count].T + residuals = np.asarray( + [residual_measure(expression, data, candidate) for candidate in constants.T], + dtype=float, + ) + return residuals + + try: + with self._sampling_seed(expression): + generator = self._generator(expression) + path = GeometricPath( + proposal=proposal, required_phi=1 / np.sqrt(len(data)) + ) + mcmc = VectorMCMC( + evaluate_model, + np.zeros(len(data)), + priors, + log_like_args=self.std, + ) + mcmc.rng = generator + sampler = AdaptiveSampler( + VectorMCMCKernel( + mcmc, + param_order=parameter_names, + path=path, + rng=generator, + ), + show_progress_bar=False, + ) + steps, marginal_log_likes = sampler.sample( + self.num_particles, + self.mcmc_steps, + self.ess_threshold, + resample_rng=lambda _, size: generator.uniform(0, 1, size), + ) + final_step = steps[-1] + map_parameters = np.asarray( + final_step.params[np.argmax(final_step.log_likes)], dtype=float + ) + map_constants = tuple(map_parameters[:parameter_count]) + phi_index = np.argmin( + np.abs(np.asarray(sampler._phi_sequence) - 1 / np.sqrt(len(data))) + ) + smc_nmll = float(marginal_log_likes[-1] - marginal_log_likes[phi_index]) + if not np.isfinite(smc_nmll) or not np.all(np.isfinite(map_constants)): + raise ValueError("SMC returned non-finite evidence or MAP constants") + except _UserCodeError as error: + raise error.error from error + except Exception as error: # SMCPy exposes several backend-specific errors. + return self._failure("sample error", error, diagnostics) + + try: + expression.commit_fit(map_constants) + except (TypeError, ValueError) as error: + return self._failure("MAP commit error", error, diagnostics) + return EvidenceResult( + smc_nmll, + map_constants, + True, + diagnostics=diagnostics, + posterior=steps if self.return_posterior else None, + ) + + def _generate_proposal(self, expression, data, residual_measure, proposal_type): + components = [] + noise_components = [] + generator = self._generator(expression) + for index in range(3 * self.num_multistarts): + trial = expression.copy() + if index: + trial.constants = np.asarray(expression.constants) + generator.normal( + scale=0.01, size=len(expression.constants) + ) + try: + result = self._call_user_code(self.fitter, trial, data, residual_measure) + if not isinstance(result, FitResult): + raise TypeError("evidence proposal fitter must return a FitResult") + constants = np.asarray(result.constants, dtype=float) + residuals = self._residuals(expression, data, residual_measure, constants) + covariance = self._laplace_covariance( + expression, data, residual_measure, constants, residuals + ) + components.append( + multivariate_normal(constants, covariance, allow_singular=True) + ) + ssqe = float(np.dot(residuals, residuals)) + variance = ssqe / len(residuals) + residual_scale = np.sqrt(np.mean(np.square(residuals))) + noise_components.append( + SqrtInvGamma( + (0.01 + len(residuals)) / 2, + max( + (0.01 * variance + ssqe) / 2, + 1e-12 * max(residual_scale, 1.0), + ), + ) + ) + except (TypeError, ValueError, np.linalg.LinAlgError): + continue + if len(components) == self.num_multistarts: + break + if not components: + raise RuntimeError("could not generate a proposal distribution") + distributions = [MixtureDistribution(*components)] + if self.std is None: + distributions.append(MixtureDistribution(*noise_components)) + return proposal_type(*distributions), {"proposal_components": len(components)} + + @staticmethod + def _residuals(expression, data, measure, constants): + residuals = np.asarray( + SmcEvidenceEstimator._call_user_code(measure, expression, data, constants), + dtype=float, + ) + if residuals.ndim != 1 or not np.all(np.isfinite(residuals)): + raise ValueError("residual measure must return finite one-dimensional values") + return residuals + + def _laplace_covariance(self, expression, data, measure, constants, residuals): + jacobian = getattr(measure, "jacobian", None) + if jacobian is None: + jacobian = self._numerical_jacobian(expression, data, measure, constants) + else: + try: + jacobian = np.asarray( + self._call_user_code(jacobian, expression, data, constants), dtype=float + ) + except AttributeError: + jacobian = self._numerical_jacobian( + expression, data, measure, constants + ) + if jacobian.shape != (len(residuals), len(constants)): + raise ValueError("residual Jacobian has an invalid shape") + curvature = 2 * jacobian.T @ jacobian + residual_hessian = getattr(measure, "residual_hessian", None) + if residual_hessian is not None: + try: + hessians = np.asarray( + self._call_user_code( + residual_hessian, expression, data, constants + ), + dtype=float, + ) + except AttributeError: + hessians = None + if hessians is None: + return self._regularized_covariance(curvature) + if hessians.shape != (len(residuals), len(constants), len(constants)): + raise ValueError("residual Hessian has an invalid shape") + curvature += 2 * np.einsum("i,ijk->jk", residuals, hessians) + return self._regularized_covariance(curvature) + + @staticmethod + def _regularized_covariance(curvature): + curvature = 0.5 * (curvature + curvature.T) + eigenvalues, eigenvectors = np.linalg.eigh(curvature) + scale = max(float(np.max(np.abs(eigenvalues))), 1.0) + eigenvalues = np.maximum(eigenvalues, scale * 1e-12) + covariance = (eigenvectors / eigenvalues) @ eigenvectors.T + return 0.5 * (covariance + covariance.T) + + def _numerical_jacobian(self, expression, data, measure, constants): + return optimize._numdiff.approx_derivative( + lambda values: self._residuals(expression, data, measure, values), constants + ) + + @staticmethod + def _call_user_code(callable_, *args): + try: + return callable_(*args) + except Exception as error: + raise _UserCodeError(error) from error + + def _generator(self, expression): + if self.seed is None: + return np.random.default_rng() + digest = hashlib.blake2b( + str(self.seed).encode() + self._structure_bytes(expression), digest_size=16 + ).digest() + return np.random.default_rng(int.from_bytes(digest, "big")) + + def _sampling_seed(self, expression): + return _SamplingSeed(self._generator(expression)) + + @staticmethod + def _structure_bytes(expression): + commands = np.asarray(expression.raw_command_array, dtype=np.uint8) + integers = np.asarray(getattr(expression, "raw_integers", ()), dtype=np.int64) + return ( + repr(commands.shape).encode() + + commands.tobytes() + + integers.tobytes() + ) + + @staticmethod + def _failure(message, error, diagnostics=None): + return EvidenceResult( + -np.inf, + None, + False, + f"{message}: {error}", + diagnostics or {}, + ) + + +class _SamplingSeed: + """Temporarily seed libraries that rely on NumPy's legacy global state.""" + + def __init__(self, generator): + self._seed = int(generator.integers(0, 2**32, dtype=np.uint32)) + self._state = None + + def __enter__(self): + self._state = np.random.get_state() + np.random.seed(self._seed) + + def __exit__(self, *_): + np.random.set_state(self._state) + + +class MixtureDistribution: + """Equal-weight mixture compatible with SMCPy's independent proposal.""" + + def __init__(self, *distributions): + self._distributions = distributions + + def rvs(self, num_samples, random_state=None): + """Draw equally from each component distribution.""" + generator = ( + np.random.default_rng(random_state) + if random_state is not None + else np.random + ) + if random_state is None: + indices = generator.randint(len(self._distributions), size=num_samples) + else: + indices = generator.integers(len(self._distributions), size=num_samples) + samples = [ + np.asarray(distribution.rvs(num_samples, random_state=random_state)).reshape( + num_samples, -1 + ) + for distribution in self._distributions + ] + return np.asarray(samples)[indices, np.arange(num_samples)] + + def logpdf(self, values): + """Return the equal-weight mixture log density.""" + densities = np.array( + [distribution.pdf(values) for distribution in self._distributions] + ) + return np.log(np.mean(densities, axis=0)).reshape(-1, 1) + + +class SqrtInvGamma: + """Square-root transformed inverse-gamma distribution for noise proposals.""" + + def __init__(self, shape, scale): + self._distribution = invgamma(shape, scale=scale) + + def rvs(self, *args, **kwargs): + """Draw standard deviations from the transformed distribution.""" + return np.sqrt(self._distribution.rvs(*args, **kwargs)) + + def pdf(self, values): + """Evaluate the transformed density.""" + values = np.asarray(values) + return 2 * values * self._distribution.pdf(np.square(values)) diff --git a/bingo/symbolic_regression/explicit_regression.py b/bingo/symbolic_regression/explicit_regression.py index 3277546d..9a58f4df 100644 --- a/bingo/symbolic_regression/explicit_regression.py +++ b/bingo/symbolic_regression/explicit_regression.py @@ -3,9 +3,10 @@ import numpy as np from ._expression_regression_objective import _ExpressionRegressionObjective +from .objective_data import ObjectiveData -class _ExplicitObjectiveData: +class _ExplicitObjectiveData(ObjectiveData): """Aligned explicit-regression arrays kept private by the objective.""" def __init__(self, X, y): @@ -15,21 +16,34 @@ def __init__(self, X, y): if self.X.ndim != 2: raise TypeError("Explicit regression X must be a 2D array") self.y = np.asarray(y, dtype=float).ravel() - if len(self.X) != len(self.y): - raise ValueError("Explicit regression X and y must have equal length") - - def __getitem__(self, items): - return _ExplicitObjectiveData(self.X[items], self.y[items]) - - def __len__(self): - return len(self.X) + super().__init__(self.X, self.y) class ExplicitRegression(_ExpressionRegressionObjective): - """Lower-is-better explicit-regression loss for evolvable Expressions.""" + """Lower-is-better explicit-regression loss for evolvable Expressions. + + Parameters + ---------- + X : array-like + Predictor values with samples along the first axis. + y : array-like + Target values aligned with ``X``. + loss : str, optional + Named Expression loss used for ranking. + fit_tolerance : float, optional + Convergence tolerance for Levenberg-Marquardt fitting. + + Raises + ------ + TypeError + If ``X`` cannot be represented as a two-dimensional array. + ValueError + If ``X`` and ``y`` have unequal sample counts. + """ def __init__(self, X, y, loss="mse", fit_tolerance=1e-5): - super().__init__(_ExplicitObjectiveData(X, y)) + data = _ExplicitObjectiveData(X, y) + super().__init__(data) self._loss = loss self._fit_tolerance = fit_tolerance diff --git a/bingo/symbolic_regression/fitting.py b/bingo/symbolic_regression/fitting.py new file mode 100644 index 00000000..993479fb --- /dev/null +++ b/bingo/symbolic_regression/fitting.py @@ -0,0 +1,323 @@ +"""Pure fitting measures and configurable SciPy-based expression fitting.""" + +# Built-in analytic expression derivatives are intentionally internal hooks. +# They are shared by the Python and C++ Expression backends. +# pylint: disable=protected-access + +from dataclasses import dataclass + +import numpy as np +from scipy import optimize + + +@dataclass(frozen=True) +class FitResult: + """Constants selected by a fitter without exposing solver-specific results. + + Parameters + ---------- + constants : array-like + Candidate simplified Expression constants. + success : bool + Whether the fitting algorithm converged numerically. + message : str, optional + Human-readable solver status. + """ + + constants: object + success: bool + message: str | None = None + + +class ResidualMeasure: + """A vector-valued fitting measure with optional residual derivatives. + + ``value(expression, data, constants)`` must return one residual per sample. + Optional ``jacobian`` and ``residual_hessian`` callables use the same + arguments and return derivatives with respect to ``constants``. + + Parameters + ---------- + value : callable + Returns a one-dimensional residual vector. + jacobian : callable, optional + Returns the residual Jacobian with respect to constants. + residual_hessian : callable, optional + Returns one constant Hessian per residual. + """ + + def __init__(self, value, *, jacobian=None, residual_hessian=None): + self.value = value + self.jacobian = jacobian + self.residual_hessian = residual_hessian + + def __call__(self, expression, data, constants): + return self.value(expression, data, constants) + + +class ScalarMeasure: + """A scalar-valued fitting measure with optional constant derivatives. + + Parameters + ---------- + value : callable + Returns a scalar fitting value. + gradient : callable, optional + Returns its gradient with respect to constants. + hessian : callable, optional + Returns its Hessian with respect to constants. + """ + + def __init__(self, value, *, gradient=None, hessian=None): + self.value = value + self.gradient = gradient + self.hessian = hessian + + def __call__(self, expression, data, constants): + return self.value(expression, data, constants) + + +def explicit_residuals(): + """Return ordinary explicit-regression residuals and their Jacobian. + + Returns + ------- + ResidualMeasure + A pure measure over two-array ``ObjectiveData(X, y)``. + """ + + def value(expression, data, constants): + X, y = data.arrays + return expression.predict(X, constants=constants) - y + + def jacobian(expression, data, constants): + X, _ = data.arrays + trial = expression.copy() + trial.constants = constants + _, derivative = trial._evaluate_with_const_gradient(np.atleast_2d(X)) + return derivative + + def residual_hessian(expression, data, constants): + X, _ = data.arrays + trial = expression.copy() + trial.constants = constants + _, _, hessian = trial._evaluate_with_const_hessian(np.atleast_2d(X)) + return hessian + + return ResidualMeasure( + value, jacobian=jacobian, residual_hessian=residual_hessian + ) + + +def implicit_residuals(*, required_params=None): + """Return normalized implicit-regression residuals. + + Parameters + ---------- + required_params : int, optional + Implicit-regression anti-triviality guard. + + Returns + ------- + ResidualMeasure + A pure measure over two-array ``ObjectiveData(X, dx_dt)``. + """ + + def value(expression, data, constants): + X, dx_dt = data.arrays + trial = expression.copy() + trial.constants = constants + _, df_dx = trial.gradient(X) + dot_product = df_dx * dx_dt + if required_params is not None: + n_params_used = (np.abs(dot_product) > 1e-16).sum(axis=1) + if not np.any(n_params_used >= required_params): + return np.full(X.shape[0], np.inf) + denominator = np.sum(np.abs(dot_product), axis=1) + with np.errstate(divide="ignore", invalid="ignore"): + residual = np.sum(dot_product, axis=1) / denominator + residual[~np.isfinite(denominator)] = np.inf + return residual + + return ResidualMeasure(value) + + +def expression_loss(kind="mse"): + """Adapt a named explicit Expression loss into a scalar fitting measure. + + Parameters + ---------- + kind : str, optional + A loss kind accepted by :meth:`Expression.loss`. + + Returns + ------- + ScalarMeasure + A pure measure over two-array ``ObjectiveData(X, y)``. + """ + + def value(expression, data, constants): + X, y = data.arrays + trial = expression.copy() + trial.constants = constants + return trial.loss(X, y, kind=kind) + + return ScalarMeasure(value) + + +def implicit_loss(*, required_params=None): + """Adapt Expression implicit loss into a scalar fitting measure. + + Parameters + ---------- + required_params : int, optional + Implicit-regression anti-triviality guard. + + Returns + ------- + ScalarMeasure + A pure measure over two-array ``ObjectiveData(X, dx_dt)``. + """ + + def value(expression, data, constants): + X, dx_dt = data.arrays + trial = expression.copy() + trial.constants = constants + return trial.implicit_loss(X, dx_dt, required_params=required_params) + + return ScalarMeasure(value) + + +class ScipyFitter: + """Fit an Expression with a selected SciPy root or minimize method. + + Root methods consume :class:`ResidualMeasure` instances. Minimize methods + consume :class:`ScalarMeasure` instances. Plain callables are classified by + their value at the fitter-owned initial constants. + + Parameters + ---------- + method : str, optional + A :func:`scipy.optimize.root` or :func:`scipy.optimize.minimize` + method. Root methods require residual measures; all other methods are + passed to ``minimize`` and require scalar measures. + tolerance : float, optional + Solver convergence tolerance. + """ + + _ROOT_METHODS = frozenset( + { + "hybr", + "lm", + "df-sane", + "broyden1", + "broyden2", + "anderson", + "linearmixing", + "diagbroyden", + "excitingmixing", + "krylov", + } + ) + + def __init__(self, method="lm", *, tolerance=1e-5): + self.method = method + self.tolerance = tolerance + + def __call__(self, expression, data, measure): + initial = np.asarray(expression.constants, dtype=float) + if self.method == "least_squares": + return self._least_squares(expression, data, measure, initial) + if self.method in self._ROOT_METHODS: + return self._root(expression, data, measure, initial) + return self._minimize(expression, data, measure, initial) + + def _least_squares(self, expression, data, measure, initial): + if isinstance(measure, ScalarMeasure): + raise TypeError("SciPy least_squares requires a residual fitting measure") + + def value(constants): + residuals = np.asarray(measure(expression, data, constants), dtype=float) + if residuals.ndim != 1: + raise ValueError("residual fitting measures must return a one-dimensional vector") + return residuals + + jacobian = getattr(measure, "jacobian", None) + if not np.all(np.isfinite(value(initial))): + return FitResult(initial, False, "non-finite residuals at initial constants") + result = optimize.least_squares( + value, + initial, + jac=( + lambda constants: np.asarray( + jacobian(expression, data, constants), dtype=float + ) + ) + if jacobian + else "2-point", + ftol=self.tolerance, + xtol=self.tolerance, + ) + return FitResult(result.x, bool(result.success), str(result.message)) + + def _root(self, expression, data, measure, initial): + if isinstance(measure, ScalarMeasure): + raise TypeError("SciPy root methods require a residual fitting measure") + + def value(constants): + residuals = np.asarray(measure(expression, data, constants), dtype=float) + if residuals.ndim != 1: + raise ValueError("residual fitting measures must return a one-dimensional vector") + return residuals + + jacobian = getattr(measure, "jacobian", None) + if not np.all(np.isfinite(value(initial))): + return FitResult(initial, False, "non-finite residuals at initial constants") + result = optimize.root( + value, + initial, + jac=( + lambda constants: np.asarray( + jacobian(expression, data, constants), dtype=float + ) + ) + if jacobian + else None, + method=self.method, + tol=self.tolerance, + ) + return FitResult(result.x, bool(result.success), str(result.message)) + + def _minimize(self, expression, data, measure, initial): + if isinstance(measure, ResidualMeasure): + raise TypeError("SciPy minimize methods require a scalar fitting measure") + + def value(constants): + result = np.asarray(measure(expression, data, constants), dtype=float) + if result.ndim != 0: + raise ValueError("scalar fitting measures must return a scalar") + return float(result) + + gradient = getattr(measure, "gradient", None) + hessian = getattr(measure, "hessian", None) + result = optimize.minimize( + value, + initial, + method=self.method, + jac=( + lambda constants: np.asarray( + gradient(expression, data, constants), dtype=float + ) + ) + if gradient + else None, + hess=( + lambda constants: np.asarray( + hessian(expression, data, constants), dtype=float + ) + ) + if hessian + else None, + tol=self.tolerance, + ) + return FitResult(result.x, bool(result.success), str(result.message)) diff --git a/bingo/symbolic_regression/implicit_regression.py b/bingo/symbolic_regression/implicit_regression.py index 3a623a4a..a21bd8d0 100644 --- a/bingo/symbolic_regression/implicit_regression.py +++ b/bingo/symbolic_regression/implicit_regression.py @@ -3,9 +3,10 @@ import numpy as np from ._expression_regression_objective import _ExpressionRegressionObjective +from .objective_data import ObjectiveData -class _ImplicitObjectiveData: +class _ImplicitObjectiveData(ObjectiveData): """Aligned implicit-regression arrays kept private by the objective.""" def __init__(self, X, dx_dt): @@ -21,23 +22,37 @@ def __init__(self, X, dx_dt): raise TypeError("Implicit regression dx_dt must be a 2D array") if self.X.shape != self.dx_dt.shape: raise ValueError("Implicit regression X and dx_dt must have equal shape") - - def __getitem__(self, items): - return _ImplicitObjectiveData(self.X[items], self.dx_dt[items]) - - def __len__(self): - return len(self.X) + super().__init__(self.X, self.dx_dt) class ImplicitRegression(_ExpressionRegressionObjective): - """Lower-is-better implicit-regression loss for evolvable Expressions.""" + """Lower-is-better implicit-regression loss for evolvable Expressions. + + Parameters + ---------- + X : array-like + State values with samples along the first axis. + dx_dt : array-like + State derivatives aligned with and shaped like ``X``. + required_params : int, optional + Minimum number of active state derivatives required to avoid a trivial + implicit solution. + + Raises + ------ + TypeError + If ``X`` or ``dx_dt`` cannot be represented as two-dimensional arrays. + ValueError + If ``X`` and ``dx_dt`` have unequal shapes. + """ def __init__(self, X, dx_dt, required_params=None): - super().__init__(_ImplicitObjectiveData(X, dx_dt)) + data = _ImplicitObjectiveData(X, dx_dt) + super().__init__(data) self._required_params = required_params def _fit_expression(self, expression, data): - expression.fit_implicit(data.X, data.dx_dt, tolerance=1e-5) + expression.fit_implicit(data.X, data.dx_dt) def _expression_loss(self, expression, data): return expression.implicit_loss( diff --git a/bingo/symbolic_regression/objective_data.py b/bingo/symbolic_regression/objective_data.py new file mode 100644 index 00000000..ccaafa74 --- /dev/null +++ b/bingo/symbolic_regression/objective_data.py @@ -0,0 +1,63 @@ +"""Aligned arrays used privately by expression regression objectives.""" + +import numpy as np + + +class ObjectiveData: + """A collection of arrays that share a sample axis. + + Parameters + ---------- + *arrays : array-like + Arrays with equal lengths along their first axes. + + Raises + ------ + TypeError + If an array has no first axis. + ValueError + If the arrays have unequal first-axis lengths. + """ + + def __init__(self, *arrays): + self._arrays = tuple(np.asarray(array) for array in arrays) + if not self._arrays: + return + try: + lengths = {array.shape[0] for array in self._arrays} + except IndexError as error: + raise TypeError("Objective data arrays must have a first axis") from error + if len(lengths) != 1: + raise ValueError("Objective data arrays must have equal first-axis length") + + @property + def arrays(self): + """The aligned arrays in their construction order.""" + return self._arrays + + def __getitem__(self, items): + """Return a collection containing each aligned array indexed by ``items``. + + Parameters + ---------- + items : int, slice, or array-like + Index applied to every aligned array. + + Returns + ------- + ObjectiveData + The indexed aligned arrays. Integer indices retain a sample axis. + """ + if isinstance(items, (int, np.integer)): + items = slice(items, items + 1) + return type(self)(*(array[items] for array in self._arrays)) + + def __len__(self): + """Return the number of aligned samples. + + Returns + ------- + int + The common first-axis length, or zero when no arrays were supplied. + """ + return len(self._arrays[0]) if self._arrays else 0 diff --git a/conda_environment.yml b/conda_environment.yml index 4aed92c8..7c5d7620 100644 --- a/conda_environment.yml +++ b/conda_environment.yml @@ -13,5 +13,3 @@ dependencies: - pip - h5py - pandas - - pip: - - smcpy diff --git a/docs/adr/0002-expression-specific-fitting.md b/docs/adr/0002-expression-specific-fitting.md index 4fb2d102..86810ac7 100644 --- a/docs/adr/0002-expression-specific-fitting.md +++ b/docs/adr/0002-expression-specific-fitting.md @@ -8,10 +8,11 @@ and fittedness; preserving a second generic protocol would duplicate those semantics without a demonstrated non-expression use case. `ExplicitRegression` and `ImplicitRegression` remain opinionated public -objectives. Advanced workflows use `CustomRegression`, which independently -selects a fitting policy, fitting measure, and ranking loss over aligned -Objective data. Evidence estimation is separate from fitting, although a -successful estimator may atomically install posterior MAP constants. +objectives and delegate their default fitting to `Expression.fit()` and +`Expression.fit_implicit()`. Advanced workflows use `CustomRegression`, which +independently selects a fitting policy, fitting measure, and ranking loss over +aligned Objective data. Evidence estimation is separate from fitting, although +a successful estimator may atomically install posterior MAP constants. ## Consequences @@ -19,4 +20,4 @@ This is a hard API switch without compatibility aliases. Generic Chromosome local fitting has no replacement, while Expression workflows gain custom SciPy fitting and SMC NMLL evidence estimation. Laplace NMLL and SMC NMLL denote the same higher-is-better normalized marginal log-likelihood concept estimated by -different methods. \ No newline at end of file +different methods. diff --git a/docs/issues/expression-fitting-and-evidence/issues/02-custom-regression-fitting.md b/docs/issues/expression-fitting-and-evidence/issues/02-custom-regression-fitting.md index bdccd89b..123ae29e 100644 --- a/docs/issues/expression-fitting-and-evidence/issues/02-custom-regression-fitting.md +++ b/docs/issues/expression-fitting-and-evidence/issues/02-custom-regression-fitting.md @@ -51,7 +51,7 @@ their current observable behavior. exposing SciPy result objects as the public contract. - [ ] Existing LM explicit fitting, implicit least-squares fitting, `required_params`, score/Loss direction, and `SymbolicRegressor` behavior are - unchanged after the preset Objectives adopt the shared machinery. + unchanged; preset Objectives continue to delegate fitting to their Expression. - [ ] Custom callable portability is documented: serial workflows accept any compatible callable, while multiprocessing and checkpoints require pickleable callables. @@ -60,4 +60,4 @@ their current observable behavior. ## Blocked by -- [01 - Build Expression fitting foundations](01-expression-fitting-foundations.md) \ No newline at end of file +- [01 - Build Expression fitting foundations](01-expression-fitting-foundations.md) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 9488a954..6d364d1e 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -7,6 +7,16 @@ To install Bingo, simply use pip: pip install bingo-nasa +Evidence Estimation +------------------- + +Install the optional Evidence-estimation dependency when using +``SmcEvidenceEstimator``: + +.. code-block:: console + + pip install "bingo-nasa[evidence]" + To use parallel island evolution, install the MPI extra: .. code-block:: console diff --git a/docs/source/migration.rst b/docs/source/migration.rst index 6ba08570..b04946c0 100644 --- a/docs/source/migration.rst +++ b/docs/source/migration.rst @@ -1,10 +1,12 @@ Expressions Migration ===================== -Bingo's symbolic-regression API is Expression-based. The only public names in -``bingo.symbolic_regression`` are ``SymbolicRegressor``, -``ExplicitRegression``, and ``ImplicitRegression``. Import expression -generation and variation from ``bingo.expressions``: +Bingo's symbolic-regression API is Expression-based. Public objectives and +fitting tools include ``SymbolicRegressor``, ``ExplicitRegression``, +``ImplicitRegression``, ``CustomRegression``, ``ObjectiveData``, +``ScipyFitter``, and ``FitResult``. Evidence estimation is provided by +``SmcEvidenceEstimator`` and ``EvidenceResult``. Import expression generation +and variation from ``bingo.expressions``: .. code-block:: python @@ -32,6 +34,32 @@ arrays directly: ``ExplicitRegression(X, y)`` and ``metric`` with ``loss``. The removed ``clo_alg`` and ``clo_threshold`` settings are replaced by expression-owned fitting configured with ``fit_tolerance``. +Generic local optimization +-------------------------- + +``bingo.local_optimizers`` has been removed with no compatibility module. +``LocalOptimizer``, ``ScipyOptimizer``, ``SmcpyOptimizer``, +``LocalOptFitnessFunction``, and ``NormalizedMarginalLikelihood`` have no +generic replacement. Generic evolutionary optimization continues to evaluate +Chromosomes directly through ``Evaluation``. + +The Chromosome local-optimization protocol +(``needs_local_optimization``, ``get_number_local_optimization_params``, +``get_local_optimization_params``, and ``set_local_optimization_params``) and +the ``needs_opt_list`` arguments to ``MultipleFloatChromosome`` and +``MultipleFloatChromosomeGenerator`` have been removed with no replacement. + +Expression users should use ``ExplicitRegression``, ``ImplicitRegression``, or +``CustomRegression`` to select an Expression fitting policy. Use +``ScipyFitter`` for configurable SciPy fitting. For Evidence estimation, use +``SmcEvidenceEstimator`` and install the optional ``evidence`` dependency. + +``VectorBasedFunction``, ``GradientMixin``, and ``VectorGradientMixin`` have +been removed. Implement custom objectives by subclassing ``FitnessFunction`` +and returning a scalar lower-is-better fitness value. Expression +``loss(..., "laplace_nmll")`` remains available for Laplace normalized marginal +log-likelihood. + Loss and serialization ---------------------- diff --git a/examples/Tutorial_1_One_Max.ipynb b/examples/Tutorial_1_One_Max.ipynb index 588c7244..43177778 100644 --- a/examples/Tutorial_1_One_Max.ipynb +++ b/examples/Tutorial_1_One_Max.ipynb @@ -274,11 +274,18 @@ " island.evolve(1)\n", " print(\"Best individual in generation\", i, \": \", island.get_best_individual())" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "pysips_dev", "language": "python", "name": "python3" }, @@ -292,7 +299,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.2" + "version": "3.13.5" } }, "nbformat": 4, diff --git a/examples/Tutorial_2_Zero_Min.ipynb b/examples/Tutorial_2_Zero_Min.ipynb index 90dc2e9f..83d39d9d 100644 --- a/examples/Tutorial_2_Zero_Min.ipynb +++ b/examples/Tutorial_2_Zero_Min.ipynb @@ -14,7 +14,7 @@ "metadata": {}, "source": [ "### Chromosome\n", - "The basic unit of bingo evolutionary analyses are Chromosomes. The chromosome used in this example is a `MultipleFloatChromosome`. The `MultipleFloatChromosome` contains a list of floating point values. It also has optional use of local optimization for some of those values." + "The basic unit of bingo evolutionary analyses are Chromosomes. The chromosome used in this example is a `MultipleFloatChromosome`, which contains a list of floating point values." ] }, { @@ -36,13 +36,13 @@ "### Chromosome Generator\n", "Chromosomes are created with a Generator. Generation of `MultipleValueChromosome` requires a function that returns floats to populate the list of values. In this example, that function is `get_random_float`.\n", "\n", - "The Generator is initialized with the random value function, along with the desired size of the float list, and an optional list of indices on which to perform local optimization. \n", + "The Generator is initialized with the random value function and desired size of the float list.\n", "The Generator is used to generate populations of Chromosomes on Islands." ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -55,7 +55,7 @@ "def get_random_float():\n", " return np.random.random_sample()\n", "\n", - "generator = MultipleFloatChromosomeGenerator(get_random_float, VALUE_LIST_SIZE, [1, 3, 4])" + "generator = MultipleFloatChromosomeGenerator(get_random_float, VALUE_LIST_SIZE)" ] }, { @@ -66,8 +66,7 @@ "source": [ "# Example of Generator\n", "chromosome = generator()\n", - "print(chromosome)\n", - "print(chromosome.get_number_local_optimization_params())" + "print(chromosome)" ] }, { @@ -80,7 +79,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -127,18 +126,16 @@ "metadata": {}, "source": [ "### Fitness and Evaluation\n", - "In order to Evaluate Chromosomes and assign them a fitness value, first we must define a `FitnessFunction`. For the Zero Min Problem, this Fitness Function calculates fitness by finding the norm of all the values in a Chromosome's list of values. Once a `FitnessFunction` has been defined, it can be passed to an Evaluation to be applied to a population. In this example, we also wrap the `FitnessFunction` with LocalOptFitnessFunction to perform local optimization on indicies specified in the Generator class." + "In order to evaluate Chromosomes and assign them a fitness value, first we must define a `FitnessFunction`. For the Zero Min Problem, this Fitness Function calculates fitness by finding the norm of all the values in a Chromosome's list of values. Once a `FitnessFunction` has been defined, it can be passed to an Evaluation to be applied to a population." ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from bingo.evaluation.fitness_function import FitnessFunction\n", - "from bingo.local_optimizers.scipy_optimizer import ScipyOptimizer\n", - "from bingo.local_optimizers.local_opt_fitness import LocalOptFitnessFunction\n", "from bingo.evaluation.evaluation import Evaluation\n", "\n", "class ZeroMinFitnessFunction(FitnessFunction):\n", @@ -147,9 +144,7 @@ "\n", " \n", "fitness = ZeroMinFitnessFunction()\n", - "optimizer = ScipyOptimizer(fitness)\n", - "local_opt_fitness = LocalOptFitnessFunction(fitness, optimizer)\n", - "evaluator = Evaluation(local_opt_fitness) # evaluates a population (list of chromosomes)" + "evaluator = Evaluation(fitness) # evaluates a population (list of chromosomes)" ] }, { @@ -159,11 +154,8 @@ "outputs": [], "source": [ "# Example of fitness\n", - "chromosome = MultipleFloatChromosome([1., 1., 1., 1., 1., 1.], \n", - " needs_opt_list=[0, 3]) # perform local optimization on these indices\n", + "chromosome = MultipleFloatChromosome([1., 1., 1., 1., 1., 1.])\n", "print(fitness(chromosome))\n", - "print(chromosome)\n", - "print(local_opt_fitness(chromosome))\n", "print(chromosome)" ] }, @@ -171,7 +163,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Notice that the values in the chromosome at indices 0 and 3 become very near zero. This occurs as part of the local optimization." + "Evolution selects chromosomes with values closer to zero over successive generations." ] }, { @@ -179,20 +171,18 @@ "metadata": {}, "source": [ "### Selection\n", - "For this example, we use Tournament Selection to select `GOAL_POPULATION_SIZE` individuals to advance to the next generation." + "For this example, we use Tournament Selection to select individuals to advance to the next generation." ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from bingo.selection.tournament import Tournament\n", "\n", - "GOAL_POPULATION_SIZE = 25\n", - "\n", - "selection = Tournament(GOAL_POPULATION_SIZE)" + "selection = Tournament(tournament_size=2)" ] }, { @@ -205,7 +195,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -213,7 +203,7 @@ "\n", "MUTATION_PROBABILITY = 0.4\n", "CROSSOVER_PROBABILITY = 0.4\n", - "NUM_OFFSPRING = GOAL_POPULATION_SIZE\n", + "NUM_OFFSPRING = 25\n", "\n", "evo_alg = MuPlusLambda(evaluator,\n", " selection,\n", @@ -234,7 +224,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -255,12 +245,12 @@ "metadata": {}, "source": [ "### Island\n", - "An `Island` is where evolution takes place in bingo analyses. The `Island` class takes as arguments an Evolutionary Algorithm, a Generator with which to generate an initial population, and thesize of the population on the island. The `Island` will create a population and then execute generational steps of the Evolutionary Algorithm to evolve the population." + "An `Island` is where evolution takes place in bingo analyses. The `Island` class takes as arguments an Evolutionary Algorithm, a Generator with which to generate an initial population, and the size of the population on the island. The `Island` will create a population and then execute generational steps of the Evolutionary Algorithm to evolve the population." ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -357,7 +347,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -372,7 +362,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -441,7 +431,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "pysips_dev", "language": "python", "name": "python3" }, @@ -455,7 +445,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.2" + "version": "3.13.5" } }, "nbformat": 4, diff --git a/examples/Tutorial_3_Archipelagos_and_Logging.ipynb b/examples/Tutorial_3_Archipelagos_and_Logging.ipynb index 3c08c598..a9747166 100644 --- a/examples/Tutorial_3_Archipelagos_and_Logging.ipynb +++ b/examples/Tutorial_3_Archipelagos_and_Logging.ipynb @@ -54,7 +54,7 @@ "def get_random_float():\n", " return np.random.random_sample()\n", "\n", - "generator = MultipleFloatChromosomeGenerator(get_random_float, VALUE_LIST_SIZE, [1, 3, 4])" + "generator = MultipleFloatChromosomeGenerator(get_random_float, VALUE_LIST_SIZE)" ] }, { @@ -95,8 +95,6 @@ "outputs": [], "source": [ "from bingo.evaluation.fitness_function import FitnessFunction\n", - "from bingo.local_optimizers.scipy_optimizer import ScipyOptimizer\n", - "from bingo.local_optimizers.local_opt_fitness import LocalOptFitnessFunction\n", "from bingo.evaluation.evaluation import Evaluation\n", "\n", "class ZeroMinFitnessFunction(FitnessFunction):\n", @@ -105,9 +103,7 @@ "\n", " \n", "fitness = ZeroMinFitnessFunction()\n", - "optimizer = ScipyOptimizer(fitness)\n", - "local_opt_fitness = LocalOptFitnessFunction(fitness, optimizer)\n", - "evaluator = Evaluation(local_opt_fitness) # evaluates a population (list of chromosomes)" + "evaluator = Evaluation(fitness) # evaluates a population (list of chromosomes)" ] }, { @@ -127,9 +123,7 @@ "source": [ "from bingo.selection.tournament import Tournament\n", "\n", - "GOAL_POPULATION_SIZE = 25\n", - "\n", - "selection = Tournament(GOAL_POPULATION_SIZE)" + "selection = Tournament(tournament_size=2)" ] }, { @@ -151,7 +145,7 @@ "\n", "MUTATION_PROBABILITY = 0.4\n", "CROSSOVER_PROBABILITY = 0.4\n", - "NUM_OFFSPRING = GOAL_POPULATION_SIZE\n", + "NUM_OFFSPRING = 25\n", "\n", "evo_alg = MuPlusLambda(evaluator,\n", " selection,\n", @@ -491,11 +485,19 @@ "from IPython.display import HTML\n", "HTML(animate_data(best_indv_values).to_jshtml())" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7b2ca8b", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pysips_dev", "language": "python", "name": "python3" }, @@ -509,7 +511,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.13" + "version": "3.13.5" }, "pycharm": { "stem_cell": { diff --git a/examples/ZeroMinExample.py b/examples/ZeroMinExample.py index bd5123aa..765e65c8 100644 --- a/examples/ZeroMinExample.py +++ b/examples/ZeroMinExample.py @@ -8,9 +8,6 @@ from bingo.selection.tournament import Tournament from bingo.evaluation.evaluation import Evaluation from bingo.evolutionary_optimizers.island import Island -from bingo.local_optimizers.scipy_optimizer import ScipyOptimizer -from bingo.local_optimizers.local_opt_fitness \ - import LocalOptFitnessFunction from bingo.chromosomes.multiple_values \ import SinglePointCrossover, SinglePointMutation from bingo.chromosomes.multiple_floats import MultipleFloatChromosomeGenerator @@ -30,9 +27,7 @@ def main(): mutation = SinglePointMutation(get_random_float) selection = Tournament(10) fitness = ZeroMinFitnessFunction() - optimizer = ScipyOptimizer(fitness) - local_opt_fitness = LocalOptFitnessFunction(fitness, optimizer) - evaluator = Evaluation(local_opt_fitness) + evaluator = Evaluation(fitness) ea = MuPlusLambda(evaluator, selection, crossover, mutation, 0.4, 0.4, 20) generator = MultipleFloatChromosomeGenerator(get_random_float, 8) island = Island(ea, generator, 25) diff --git a/examples/get_started.ipynb b/examples/get_started.ipynb index 394ca903..9b7c481a 100644 --- a/examples/get_started.ipynb +++ b/examples/get_started.ipynb @@ -39,7 +39,7 @@ "source": [ "import numpy as np\n", "X = np.linspace(-10, 10).reshape((-1, 1))\n", - "y = (3.0 * X ** 2 - 2.0 * X).ravel()" + "y = (3.0 * X ** 2 - 2.0 * X).ravel()" ] }, { @@ -115,7 +115,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pysips_dev", "language": "python", "name": "python3" }, @@ -129,7 +129,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.5" + "version": "3.13.5" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index 6f6fe58b..cb419813 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,6 @@ dependencies = [ "dill", "sympy", "scikit-learn", - "smcpy", "pybind11[global]", ] classifiers = [ @@ -43,6 +42,7 @@ Documentation = "https://nasa.github.io/bingo/" Repository = "https://github.com/nasa/bingo" [project.optional-dependencies] +evidence = ["smcpy"] MPI = ["mpi4py>=4.0"] ONNX = ["onnx"] TESTS = [ diff --git a/requirements.txt b/requirements.txt index 7db40728..0ade81f0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,7 +5,6 @@ scipy>=1.6.2 dill>=0.2.9 sympy>=1.0 scikit-learn>=1.1 -smcpy>=0.1.4 pybind11[global] # ONNX onnx diff --git a/setup.py b/setup.py index 87e92798..6a0d9101 100644 --- a/setup.py +++ b/setup.py @@ -184,7 +184,6 @@ def build_extension(self, ext: CMakeExtension) -> None: "bingo.expressions.agraph.pyagraph.evaluation", "bingo.expressions.agraph.pyagraph.simplification", "bingo.expressions.agraph.cppagraph", - "bingo.local_optimizers", "bingo.selection", "bingo.stats", "bingo.symbolic_regression", diff --git a/tests/integration/test_clo_optim.py b/tests/integration/test_clo_optim.py deleted file mode 100644 index c9fd9654..00000000 --- a/tests/integration/test_clo_optim.py +++ /dev/null @@ -1,116 +0,0 @@ -# Ignoring some linting rules in tests -# pylint: disable=redefined-outer-name -# pylint: disable=missing-docstring -import pytest -import numpy as np - -from bingo.evaluation.fitness_function \ - import FitnessFunction, VectorBasedFunction -from bingo.evaluation.gradient_mixin import GradientMixin, VectorGradientMixin -from bingo.local_optimizers.scipy_optimizer import ScipyOptimizer, \ - MINIMIZE_SET, ROOT_SET -from bingo.local_optimizers.local_opt_fitness \ - import LocalOptFitnessFunction -from bingo.chromosomes.multiple_floats import MultipleFloatChromosome - -NUM_VALS = 10 -NUM_OPT = 3 - - -class MultipleFloatValueFitnessFunction(FitnessFunction): - def __call__(self, individual): - return np.linalg.norm(individual.values) - - -class MultipleFloatValueFitnessFunctionWithGradient( - GradientMixin, MultipleFloatValueFitnessFunction): - def get_fitness_and_gradient(self, individual): - full_gradient = individual.values / np.linalg.norm(individual.values) - return self.__call__(individual), \ - [full_gradient[i] for i in individual._needs_opt_list] - - -class FloatVectorFitnessFunction(VectorBasedFunction): - def evaluate_fitness_vector(self, individual): - vals = individual.values - return [x - 0 for x in vals] - - -class FloatVectorFitnessFunctionWithJacobian(VectorGradientMixin, - FloatVectorFitnessFunction): - def get_fitness_vector_and_jacobian(self, individual): - jacobian = np.zeros((len(individual.values), - len(individual._needs_opt_list))) - for i, optimize_i in enumerate(individual._needs_opt_list): - jacobian[optimize_i][i] = 1 - return self.evaluate_fitness_vector(individual), jacobian - - -@pytest.fixture -def opt_individual(): - vals = [1. for _ in range(NUM_VALS)] - return MultipleFloatChromosome(vals, [1, 3, 4]) - - -@pytest.fixture -def reg_individual(): - vals = [1. for _ in range(NUM_VALS)] - return MultipleFloatChromosome(vals) - - -@pytest.mark.parametrize("method", MINIMIZE_SET) -def test_optimize_params_without_gradient( - opt_individual, reg_individual, method): - np.random.seed(2) - fitness_function = MultipleFloatValueFitnessFunction() - local_opt_fitness_function = LocalOptFitnessFunction( - fitness_function, ScipyOptimizer(fitness_function, method=method)) - opt_indv_fitness = local_opt_fitness_function(opt_individual) - reg_indv_fitness = local_opt_fitness_function(reg_individual) - assert opt_indv_fitness == pytest.approx(np.sqrt(NUM_VALS - NUM_OPT), - rel=5.e-6) - assert reg_indv_fitness == pytest.approx(np.sqrt(NUM_VALS)) - - -@pytest.mark.parametrize("method", MINIMIZE_SET) -def test_optimize_params_with_gradient(opt_individual, reg_individual, method): - np.random.seed(2) - fitness_function = MultipleFloatValueFitnessFunctionWithGradient() - local_opt_fitness_function = LocalOptFitnessFunction( - fitness_function, ScipyOptimizer(fitness_function, method=method)) - opt_indv_fitness = local_opt_fitness_function(opt_individual) - reg_indv_fitness = local_opt_fitness_function(reg_individual) - assert opt_indv_fitness == pytest.approx(np.sqrt(NUM_VALS - NUM_OPT), - rel=5.e-6) - assert reg_indv_fitness == pytest.approx(np.sqrt(NUM_VALS)) - - -@pytest.mark.parametrize("method", ROOT_SET) -def test_optimize_fitness_vector_without_jacobian(opt_individual, - reg_individual, method): - reg_list = [1. for _ in range(NUM_VALS)] - opt_list = [1. for _ in range(NUM_VALS)] - opt_list[:3] = [0., 0., 0.] - fitness_function = FloatVectorFitnessFunction() - local_opt_fitness_function = LocalOptFitnessFunction( - fitness_function, ScipyOptimizer(fitness_function, method=method)) - opt_indv_fitness = local_opt_fitness_function(opt_individual) - reg_indv_fitness = local_opt_fitness_function(reg_individual) - assert opt_indv_fitness == pytest.approx(np.mean(opt_list)) - assert reg_indv_fitness == pytest.approx(np.mean(reg_list)) - - -@pytest.mark.parametrize("method", ROOT_SET) -def test_optimize_fitness_vector_with_jacobian(opt_individual, reg_individual, - method): - np.random.seed(0) - reg_list = [1. for _ in range(NUM_VALS)] - opt_list = [1. for _ in range(NUM_VALS)] - opt_list[:3] = [0., 0., 0.] - fitness_function = FloatVectorFitnessFunctionWithJacobian() - local_opt_fitness_function = LocalOptFitnessFunction( - fitness_function, ScipyOptimizer(fitness_function, method=method)) - opt_indv_fitness = local_opt_fitness_function(opt_individual) - reg_indv_fitness = local_opt_fitness_function(reg_individual) - assert opt_indv_fitness == pytest.approx(np.mean(opt_list)) - assert reg_indv_fitness == pytest.approx(np.mean(reg_list)) diff --git a/tests/unit/chromosomes/test_chromosome.py b/tests/unit/chromosomes/test_chromosome.py index 1bb627f1..fe209724 100644 --- a/tests/unit/chromosomes/test_chromosome.py +++ b/tests/unit/chromosomes/test_chromosome.py @@ -32,23 +32,3 @@ def test_genetic_age_starts_at_zero(individual): assert individual.genetic_age == 0 individual.genetic_age = 10 assert individual.genetic_age == 10 - - -def test_optimization_interface_methods_raise_not_implemented(mocker): - mocker.patch.object(Chromosome, "__abstractmethods__", new_callable=set) - - expected_exception_str = ( - "This Chromosome cannot be used in local " - "optimization until its local optimization " - "interface has been implemented" - ) - - with pytest.raises(NotImplementedError) as exc_info: - Chromosome().needs_local_optimization() - assert expected_exception_str == str(exc_info.value) - - assert Chromosome().get_number_local_optimization_params() == 0 - - with pytest.raises(NotImplementedError) as exc_info: - Chromosome().set_local_optimization_params(mocker.Mock()) - assert expected_exception_str == str(exc_info.value) diff --git a/tests/unit/chromosomes/test_multiple_float_chromosome.py b/tests/unit/chromosomes/test_multiple_float_chromosome.py index 2edd3300..7356f748 100644 --- a/tests/unit/chromosomes/test_multiple_float_chromosome.py +++ b/tests/unit/chromosomes/test_multiple_float_chromosome.py @@ -1,65 +1,18 @@ # Ignoring some linting rules in tests # pylint: disable=redefined-outer-name # pylint: disable=missing-docstring -import pytest - from bingo.chromosomes.multiple_floats import MultipleFloatChromosome,\ - MultipleFloatChromosomeGenerator + MultipleFloatChromosomeGenerator DUMMY_VALUE = 999 -def test_multiple_float_needs_local_optimization(): - chromosome_with_opt = MultipleFloatChromosome([1, 2, 3], [1]) - chromosome_without_opt = MultipleFloatChromosome([1, 2, 3]) - - assert chromosome_with_opt. needs_local_optimization() - assert not chromosome_without_opt. needs_local_optimization() - - -@pytest.mark.parametrize("num_params", range(3)) -def test_getting_number_of_optimization_params(num_params): - needs_opt_list = list(range(num_params)) - chromosome = MultipleFloatChromosome([1, 2, 3], needs_opt_list) - assert chromosome.get_number_local_optimization_params() == num_params - - -def test_setting_optimization_params(): - needs_opt_list = [1, 3, 5] - chromosome = MultipleFloatChromosome([0] * 6, needs_opt_list) - chromosome.set_local_optimization_params([1, 1, 1]) - for i in needs_opt_list: - assert chromosome.values[i] == 1 - - -@pytest.mark.parametrize("bad_opt_list", [[-1, 0], - [0, 4], - [0, 0.5]]) -def test_generator_errors_with_bad_opt_list(mocker, bad_opt_list): - with pytest.raises(ValueError): - _ = MultipleFloatChromosomeGenerator(mocker.Mock(), - values_per_chromosome=4, - needs_opt_list=bad_opt_list) - - def test_generator(): def dummy_function(): return DUMMY_VALUE - needs_opt_list = [1, 3, 5] - generator = MultipleFloatChromosomeGenerator(dummy_function, - values_per_chromosome=6, - needs_opt_list=needs_opt_list) - chromosome = generator() - chromosome.set_local_optimization_params([1, 1, 1]) - for i in needs_opt_list: - assert chromosome.values[i] == 1 - -def test_generator_default(): - def dummy_function(): - return DUMMY_VALUE generator = MultipleFloatChromosomeGenerator(dummy_function, - values_per_chromosome=6) + values_per_chromosome=6) chromosome = generator() - assert chromosome.get_number_local_optimization_params() == 0 + assert chromosome.values == [DUMMY_VALUE] * 6 diff --git a/tests/unit/evaluation/test_fitness_function.py b/tests/unit/evaluation/test_fitness_function.py index fa580c95..0e44c6af 100644 --- a/tests/unit/evaluation/test_fitness_function.py +++ b/tests/unit/evaluation/test_fitness_function.py @@ -1,27 +1,10 @@ -"""Tests for generic Python fitness aggregation.""" - -import numpy as np +"""Tests for the fitness-function base class.""" import pytest -from bingo.evaluation.fitness_function import FitnessFunction, VectorBasedFunction +from bingo.evaluation.fitness_function import FitnessFunction from bingo.evaluation.training_data import TrainingData -class _Individual: - def get_number_local_optimization_params(self): - return 2 - - -class _VectorFitness(VectorBasedFunction): - def evaluate_fitness_vector(self, individual): - return np.array([-2.0, -1.0, 0.0, 1.0, 2.0]) - - -class _NanVectorFitness(VectorBasedFunction): - def evaluate_fitness_vector(self, individual): - return np.array([np.nan, -1.0, 0.0, 1.0, 2.0]) - - def test_fitness_function_cannot_be_instantiated(): with pytest.raises(TypeError): FitnessFunction() @@ -36,42 +19,3 @@ def test_fitness_function_stores_training_data(mocker): assert fitness.eval_count == 0 assert fitness.training_data is training_data - - -@pytest.mark.parametrize( - "metric, expected", - [ - ("mae", 1.2), - ("mean absolute error", 1.2), - ("mse", 2.0), - ("mean squared error", 2.0), - ("rmse", np.sqrt(2.0)), - ("root mean squared error", np.sqrt(2.0)), - ("negative nmll laplace", 6.0868339), - ("bic", 22.483434972148757), - ], -) -def test_vector_fitness_aggregates_its_error_vector(metric, expected): - assert _VectorFitness(metric=metric)(_Individual()) == pytest.approx(expected) - - -def test_vector_fitness_rejects_unknown_metric(): - with pytest.raises(ValueError): - _VectorFitness(metric="unknown") - - -@pytest.mark.parametrize( - "metric", - [ - "mae", - "mean absolute error", - "mse", - "mean squared error", - "rmse", - "root mean squared error", - "negative nmll laplace", - "bic", - ], -) -def test_vector_fitness_propagates_nan(metric): - assert np.isnan(_NanVectorFitness(metric=metric)(_Individual())) diff --git a/tests/unit/evaluation/test_gradient_mixin.py b/tests/unit/evaluation/test_gradient_mixin.py deleted file mode 100644 index d9f871c2..00000000 --- a/tests/unit/evaluation/test_gradient_mixin.py +++ /dev/null @@ -1,65 +0,0 @@ -"""Tests for generic Python gradient-based fitness aggregation.""" - -import numpy as np -import pytest - -from bingo.evaluation.fitness_function import VectorBasedFunction -from bingo.evaluation.gradient_mixin import GradientMixin, VectorGradientMixin - - -class _Individual: - def get_number_local_optimization_params(self): - return 2 - - -class _GradientFitness(VectorGradientMixin, VectorBasedFunction): - def get_fitness_vector_and_jacobian(self, individual): - return np.array([-2.0, 0.0, 2.0]), np.array( - [[0.5, 1.0], [1.0, 2.0], [-0.5, 3.0]] - ) - - def evaluate_fitness_vector(self, individual): - return self.get_fitness_vector_and_jacobian(individual)[0] - - -def test_gradient_mixin_cannot_be_instantiated(): - with pytest.raises(TypeError): - GradientMixin() - - -def test_vector_gradient_mixin_requires_vector_fitness_base(): - class _InvalidGradientFitness(VectorGradientMixin): - def get_fitness_vector_and_jacobian(self, individual): - return None - - with pytest.raises(TypeError): - _InvalidGradientFitness() - - -@pytest.mark.parametrize( - "metric, expected_fitness, expected_gradient", - [ - ("mae", 4 / 3, [-1 / 3, 2 / 3]), - ("mean absolute error", 4 / 3, [-1 / 3, 2 / 3]), - ("mse", 8 / 3, [-4 / 3, 8 / 3]), - ("mean squared error", 8 / 3, [-4 / 3, 8 / 3]), - ("rmse", np.sqrt(8 / 3), [-np.sqrt(3 / 8) * 2 / 3, np.sqrt(3 / 8) * 4 / 3]), - ("root mean squared error", np.sqrt(8 / 3), [-np.sqrt(3 / 8) * 2 / 3, np.sqrt(3 / 8) * 4 / 3]), - ("negative nmll laplace", 3.244922013421868, [-0.3169873, 0.6339746]), - ("bic", 14.751955824267544, [-1.5, 3.0]), - ], -) -def test_vector_gradient_fitness_aggregates_vector_and_jacobian( - metric, expected_fitness, expected_gradient -): - fitness, gradient = _GradientFitness(metric=metric).get_fitness_and_gradient( - _Individual() - ) - - assert fitness == pytest.approx(expected_fitness) - np.testing.assert_allclose(gradient, expected_gradient) - - -def test_vector_gradient_mixin_rejects_unknown_metric(): - with pytest.raises(ValueError): - _GradientFitness(metric="unknown") diff --git a/tests/unit/expressions/agraph/test_fit_lifecycle.py b/tests/unit/expressions/agraph/test_fit_lifecycle.py new file mode 100644 index 00000000..6a6523a1 --- /dev/null +++ b/tests/unit/expressions/agraph/test_fit_lifecycle.py @@ -0,0 +1,61 @@ +import copy +import pickle + +import numpy as np +import pytest + +from bingo.expressions.agraph.pyagraph import AGraphExpression as PyAGraphExpression + +try: + from bingo.expressions.agraph.cppagraph import AGraphExpression as CppAGraphExpression +except ImportError: + EXPRESSION_TYPES = [PyAGraphExpression] +else: + EXPRESSION_TYPES = [PyAGraphExpression, CppAGraphExpression] + + +@pytest.mark.parametrize("expression_type", EXPRESSION_TYPES) +def test_commit_fit_validates_atomically_and_clear_fit_resets(expression_type): + expression = expression_type(equation="X0 + 1.0") + + with pytest.raises(ValueError, match="one entry per simplified expression constant"): + expression.commit_fit([]) + assert expression.constants == (1.0,) + assert not expression.is_fitted + + with pytest.raises(ValueError, match="finite"): + expression.commit_fit([np.nan]) + assert expression.constants == (1.0,) + assert not expression.is_fitted + + assert expression.commit_fit([2.0]) is expression + assert expression.constants == (2.0,) + assert expression.is_fitted + + assert expression.clear_fit() is expression + assert not expression.is_fitted + + +@pytest.mark.parametrize("expression_type", EXPRESSION_TYPES) +def test_fit_commit_lifecycle_survives_copy_and_pickle(expression_type): + expression = expression_type(equation="X0 + 1.0") + expression.commit_fit([2.0]) + + copied = copy.deepcopy(expression) + restored = pickle.loads(pickle.dumps(expression)) + + for clone in (copied, restored): + assert clone.constants == (2.0,) + assert clone.is_fitted + + expression.raw_constants = (3.0,) + assert not expression.is_fitted + + +@pytest.mark.parametrize("expression_type", EXPRESSION_TYPES) +def test_clear_fit_is_a_noop_for_constant_free_expressions(expression_type): + expression = expression_type(equation="X0") + + assert expression.is_fitted + assert expression.clear_fit() is expression + assert expression.is_fitted diff --git a/tests/unit/expressions/test_evolvable.py b/tests/unit/expressions/test_evolvable.py index 116895fe..de0537ce 100644 --- a/tests/unit/expressions/test_evolvable.py +++ b/tests/unit/expressions/test_evolvable.py @@ -107,13 +107,6 @@ def test_different_arrays_positive_distance(self): assert a.distance(b) > 0 -class TestNoLocalOptimizationAdapter: - def test_does_not_implement_local_optimization_methods(self): - assert "needs_local_optimization" not in EvolvableExpression.__dict__ - assert "get_number_local_optimization_params" not in EvolvableExpression.__dict__ - assert "set_local_optimization_params" not in EvolvableExpression.__dict__ - - class TestDelegation: def test_user_facing_facade(self, backend): expression = agraph.get_expression_class()( diff --git a/tests/unit/local_optimizers/test_local_opt.py b/tests/unit/local_optimizers/test_local_opt.py deleted file mode 100644 index d0a61b65..00000000 --- a/tests/unit/local_optimizers/test_local_opt.py +++ /dev/null @@ -1,126 +0,0 @@ -from bingo.local_optimizers.local_opt_fitness \ - import LocalOptFitnessFunction -from bingo.chromosomes.chromosome import Chromosome - - -class DummyLocalOptimizationIndividual(Chromosome): - def __init__(self): - super().__init__() - self._params = [1, 2, 3] - self._needs_opt = True - - def needs_local_optimization(self): - return self._needs_opt - - def get_params(self): - return self._params - - def get_number_local_optimization_params(self): - return len(self._params) - - def __str__(self): - pass - - def distance(self, chromosome): - pass - - def set_local_optimization_params(self, params): - self._params = params - self._needs_opt = False - - -def test_get_eval_count_pass_through(mocker): - fitness_function = mocker.Mock() - fitness_function.eval_count = 123 - optimizer = mocker.Mock() - local_opt_fitness_function = \ - LocalOptFitnessFunction(fitness_function, optimizer) - assert local_opt_fitness_function.eval_count == 123 - - -def test_set_eval_count_pass_through(mocker): - fitness_function = mocker.Mock() - optimizer = mocker.Mock() - local_opt_fitness_function = \ - LocalOptFitnessFunction(fitness_function, optimizer) - local_opt_fitness_function.eval_count = 123 - assert fitness_function.eval_count == 123 - - -def test_get_training_data_pass_through(mocker): - fitness_function = mocker.Mock() - fitness_function.training_data = 123 - optimizer = mocker.Mock() - local_opt_fitness_function = \ - LocalOptFitnessFunction(fitness_function, optimizer) - assert local_opt_fitness_function.training_data == 123 - - -def test_set_training_data_pass_through(mocker): - fitness_function = mocker.Mock() - optimizer = mocker.Mock() - local_opt_fitness_function = \ - LocalOptFitnessFunction(fitness_function, optimizer) - local_opt_fitness_function.training_data = 123 - assert fitness_function.training_data == 123 - - -def test_get_and_set_optimizer(mocker): - fitness_function = mocker.Mock() - opt_1 = mocker.Mock() - opt_2 = mocker.Mock() - local_opt_fitness_function = \ - LocalOptFitnessFunction(fitness_function, opt_1) - assert local_opt_fitness_function.optimizer == opt_1 - - local_opt_fitness_function.optimizer = opt_2 - assert local_opt_fitness_function.optimizer == opt_2 - - -def test_call_optimizes_when_necessary(mocker): - fitness_function = mocker.Mock( - side_effect=lambda ind: sum(ind.get_params()) - ) - optimizer = mocker.Mock( - side_effect=lambda ind: ind.set_local_optimization_params([4, 5, 6]) - ) - - individual = DummyLocalOptimizationIndividual() - - local_opt_fitness_function = \ - LocalOptFitnessFunction(fitness_function, optimizer) - - returned_fitness = local_opt_fitness_function(individual) - - # make sure optimizer was called with individual - optimizer.assert_called_once_with(individual) - assert individual.get_params() == [4, 5, 6] - - # make sure that fitness function was called with individual - fitness_function.assert_called_once_with(individual) - assert returned_fitness == fitness_function(individual) - - -def test_call_doesnt_optimize_when_not_needed(mocker): - fitness_function = mocker.Mock( - side_effect=lambda ind: sum(ind.get_params()) - ) - optimizer = mocker.Mock( - side_effect=lambda ind: ind.set_local_optimization_params([4, 5, 6]) - ) - - individual = DummyLocalOptimizationIndividual() - initial_params = [1, 2, 3] - individual.set_local_optimization_params(initial_params) - - local_opt_fitness_function = \ - LocalOptFitnessFunction(fitness_function, optimizer) - - returned_fitness = local_opt_fitness_function(individual) - - assert not optimizer.called - assert individual.get_params() == initial_params - - # make sure that fitness function was called with individual - fitness_function.assert_called_once_with(individual) - assert returned_fitness == fitness_function(individual) diff --git a/tests/unit/local_optimizers/test_normalized_marginal_likelihood.py b/tests/unit/local_optimizers/test_normalized_marginal_likelihood.py deleted file mode 100644 index d5841525..00000000 --- a/tests/unit/local_optimizers/test_normalized_marginal_likelihood.py +++ /dev/null @@ -1,49 +0,0 @@ -import pytest -from bingo.local_optimizers.normalized_marginal_likelihood import ( - NormalizedMarginalLikelihood, -) - - -def test_nml_sets_up_smcpy_optimizer(mocker): - objective_fn = mocker.MagicMock() - deterministic_optimizer = mocker.Mock() - smcpy_optimizer = mocker.patch( - "bingo.local_optimizers.normalized_marginal_likelihood.SmcpyOptimizer", - autospec=True, - ) - - nml = NormalizedMarginalLikelihood( - objective_fn, deterministic_optimizer, log_scale=True, mcmc_steps=999 - ) - nml.training_data = 10 - nml.eval_count = 13 - - smcpy_optimizer.assert_called_once_with( - objective_fn, deterministic_optimizer, mcmc_steps=999 - ) - - assert nml.optimizer.training_data == 10 - assert nml.optimizer.eval_count == 13 - assert nml.training_data == 10 - assert nml.eval_count == 13 - - -@pytest.mark.parametrize( - "log_scale,expected_nml", [(True, -1), (False, -2.718281828459045)] -) -def test_logscale(mocker, log_scale, expected_nml): - individual = mocker.Mock() - objective_fn = mocker.MagicMock() - deterministic_optimizer = mocker.Mock() - mocker.patch( - "bingo.local_optimizers.normalized_marginal_likelihood.SmcpyOptimizer", - autospec=True, - return_value=mocker.Mock(return_value=(1, None, None)), - ) - - nml = NormalizedMarginalLikelihood( - objective_fn, deterministic_optimizer, log_scale=log_scale - ) - - assert nml(individual) == expected_nml - diff --git a/tests/unit/local_optimizers/test_scipy_optimizer.py b/tests/unit/local_optimizers/test_scipy_optimizer.py deleted file mode 100644 index c825523b..00000000 --- a/tests/unit/local_optimizers/test_scipy_optimizer.py +++ /dev/null @@ -1,318 +0,0 @@ -# Ignoring some linting rules in tests -# pylint: disable=redefined-outer-name -# pylint: disable=missing-docstring -from copy import deepcopy - -import pytest -import numpy as np - -from scipy import optimize -from scipy.optimize import OptimizeResult - -from bingo.evaluation.fitness_function \ - import FitnessFunction, VectorBasedFunction -from bingo.evaluation.gradient_mixin import GradientMixin, VectorGradientMixin -from bingo.chromosomes.chromosome import Chromosome -from bingo.local_optimizers.scipy_optimizer import ScipyOptimizer, \ - ROOT_SET, MINIMIZE_SET, JACOBIAN_SET - - -class DummyLocalOptIndividual(Chromosome): - def needs_local_optimization(self): - return True - - def get_number_local_optimization_params(self): - return 1 - - def set_local_optimization_params(self, params): - self.param = params[0] - - def __str__(self): - pass - - def distance(self, chromosome): - pass - - -class BloatedOptIndividual(Chromosome): - def __init__(self): - super().__init__() - self.param = [1, 2, 3] - self._fitness = None - self._fit_set = False - - def needs_local_optimization(self): - return True - - def get_number_local_optimization_params(self): - return 3 - - def set_local_optimization_params(self, params): - self.param = params - - def __str__(self): - pass - - def distance(self, chromosome): - pass - - -@pytest.mark.parametrize("obj_func_type, raises_error", - [(FitnessFunction, True), - (VectorBasedFunction, False)]) -def test_valid_objective_function_init(mocker, obj_func_type, raises_error): - mocked_obj_fn = mocker.create_autospec(obj_func_type) - if raises_error: - with pytest.raises(TypeError): - _ = ScipyOptimizer(mocked_obj_fn, method="lm") - else: - _ = ScipyOptimizer(mocked_obj_fn, method="lm") - - -@pytest.mark.parametrize("obj_func_type, raises_error", - [(FitnessFunction, True), - (VectorBasedFunction, False)]) -def test_valid_objective_function_property(mocker, obj_func_type, raises_error): - mocked_obj_fn = mocker.create_autospec(obj_func_type) - - # construct obj with valid objective fn - valid_obj_fn = mocker.create_autospec(VectorBasedFunction) - opt = ScipyOptimizer(valid_obj_fn, method="lm") - - # test setting obj fn with property - if raises_error: - with pytest.raises(TypeError): - opt.objective_fn = mocked_obj_fn - else: - opt.objective_fn = mocked_obj_fn - - -def test_invalid_method(mocker): - mocked_objective_function = mocker.Mock() - with pytest.raises(KeyError): - ScipyOptimizer(mocked_objective_function, - method="Dwayne - The Rock - Johnson") - - -def get_expected_options(**additional_options): - default_options = {"method": "BFGS", - "tol": 1e-6, - "param_init_bounds": [-10000, 10000]} - default_options.update(additional_options) - return default_options - - -def test_can_set_options_via_constructor(mocker): - mock_obj_fn = \ - mocker.Mock(side_effect=lambda individual: individual.param) - - # testing default options - opt = ScipyOptimizer(mock_obj_fn) - assert opt.options == get_expected_options() - - # testing adding extra options - opt = ScipyOptimizer(mock_obj_fn, - options={"maxiter": 0}) - assert opt.options == get_expected_options(options={"maxiter": 0}) - - # testing setting default options - opt = ScipyOptimizer(mock_obj_fn, - method="lm", - param_init_bounds=[-1, 1], - tol=1e-8) - assert opt.options == get_expected_options( - method="lm", - param_init_bounds=[-1, 1], - tol=1e-8 - ) - - -def test_can_set_options_via_property(mocker): - mock_obj_fn = \ - mocker.Mock(side_effect=lambda individual: individual.param) - opt = ScipyOptimizer(mock_obj_fn) - - # testing default options - assert opt.options == get_expected_options() - - # testing adding extra options - opt_options = {"options": {"maxiter": 0}} - opt.options = opt_options - assert opt.options == get_expected_options(options={"maxiter": 0}) - - # testing removing extra options/setting default options - opt_options = {"method": "lm", - "param_init_bounds": [-1, 1], - "tol": 1e-8} - opt.options = opt_options - assert opt.options == get_expected_options( - method="lm", - param_init_bounds=[-1, 1], - tol=1e-8 - ) - - -@pytest.mark.parametrize("method", ["Nelder-Mead", "lm"]) -# using Nelder-Mead and lm to test minimize and root respectively -def test_set_param_bounds_and_clo_options_affect_clo(mocker, method): - mocked_fitness_function = \ - mocker.Mock(side_effect=lambda individual: individual.param) - - dummy_individual = DummyLocalOptIndividual() - - opt_options = {"method": method, # TODO have to set this or else it will get overwritten by BFGS on update, I don't like this - "tol": 1e-8, - "options": {"maxiter": 1000, - "fatol": 1e-8, - "xatol": 1e-8, - "adaptive": False}} - - expected_options = deepcopy(opt_options) - expected_options["method"] = method - expected_options["args"] = dummy_individual - expected_options["jac"] = False - - # returns x=0 if kwargs != expected_options and x=1 vice versa - def mocked_optimize(*_, **kwargs): - return OptimizeResult(x=[int(kwargs == expected_options)]) - - mocker.patch.object(optimize, "minimize", side_effect=mocked_optimize) - mocker.patch.object(optimize, "root", side_effect=mocked_optimize) - - opt = ScipyOptimizer(mocked_fitness_function, - method=method) - - opt(dummy_individual) - # default options should != opt_options - assert dummy_individual.param == 0 - - opt.options = opt_options - opt(dummy_individual) - # we set opt.options, so options should == opt_options - assert dummy_individual.param == 1 - - -@pytest.mark.parametrize("method", MINIMIZE_SET) -def test_optimize_params_minimize_without_gradient(mocker, method): - fitness_function = mocker.create_autospec(FitnessFunction) - fitness_function.side_effect = lambda individual: 1 + individual.param ** 2 - - opt = ScipyOptimizer(fitness_function, - param_init_bounds=[5, 5], - method=method) - - individual = DummyLocalOptIndividual() - opt(individual) - assert fitness_function(individual) == pytest.approx(1, rel=0.05) - - -class GradientFitnessFunction(GradientMixin, FitnessFunction): - def __call__(self, individual): - pass - - def get_fitness_and_gradient(self, individual): - pass - - -@pytest.mark.parametrize("method", MINIMIZE_SET) -def test_optimize_params_minimize_with_gradient(mocker, method): - fitness_function = mocker.create_autospec(GradientFitnessFunction) - - mocked_fitness = mocker.Mock(side_effect=lambda x: 1 + x.param ** 2) - fitness_function.side_effect = mocked_fitness - - mocked_gradient = \ - mocker.Mock(side_effect=lambda x: (1 + x.param ** 2, 2 * x.param)) - fitness_function.get_fitness_and_gradient = mocked_gradient - - opt = ScipyOptimizer(fitness_function, - param_init_bounds=[5, 5], - method=method) - - individual = DummyLocalOptIndividual() - opt(individual) - assert mocked_fitness.called - if method in JACOBIAN_SET: - assert mocked_gradient.called - assert fitness_function(individual) == pytest.approx(1, rel=0.05) - - -class JacobianVectorFitnessFunction(VectorGradientMixin, VectorBasedFunction): - def evaluate_fitness_vector(self, individual): - pass - - def get_fitness_vector_and_jacobian(self, individual): - pass - - -@pytest.mark.parametrize("method", ROOT_SET) -def test_optimize_params_root_without_jacobian(mocker, method): - fitness_function = mocker.create_autospec(VectorBasedFunction) - fitness_function.evaluate_fitness_vector = lambda x: 1 + np.abs([x.param]) - - opt = ScipyOptimizer(fitness_function, - param_init_bounds=[5, 5], - method=method) - - individual = DummyLocalOptIndividual() - opt(individual) - opt_indv_fitness = fitness_function.evaluate_fitness_vector(individual) - assert opt_indv_fitness[0] == pytest.approx(1, rel=0.05) - - -@pytest.mark.parametrize("method", ROOT_SET) -def test_optimize_params_root_with_jacobian(mocker, method): - fitness_function = mocker.create_autospec(JacobianVectorFitnessFunction) - - mocked_fitness_vector = \ - mocker.Mock(side_effect=lambda x: 1 + np.abs([x.param])) - fitness_function.evaluate_fitness_vector = mocked_fitness_vector - - mocked_jacobian = mocker.Mock( - side_effect=lambda x: (1 + np.abs([x.param]), np.sign([x.param]))) - fitness_function.get_fitness_vector_and_jacobian = mocked_jacobian - - - opt = ScipyOptimizer(fitness_function, - param_init_bounds=[5, 5], - method=method) - - individual = DummyLocalOptIndividual() - opt(individual) - assert mocked_fitness_vector.called - if method in JACOBIAN_SET: - assert mocked_jacobian.called - opt_indv_fitness = fitness_function.evaluate_fitness_vector(individual) - assert opt_indv_fitness[0] == pytest.approx(1, rel=0.05) - - -# easier to do this than to mock a VectorBasedFunction -# because you can't patch __abstractmethods__ with the cpp version -class ReductionFunction(VectorBasedFunction): - def evaluate_fitness_vector(self, individual): - return np.array([1.0 + sum(np.square(individual.param))]) - - -@pytest.mark.parametrize("method", ROOT_SET) -def test_optimize_params_too_many_params(method): - # need a fitness function that has less entries than params - fitness_function = ReductionFunction() - - opt = ScipyOptimizer(fitness_function, - method=method, - param_init_bounds=[-1, -1]) - - individual = BloatedOptIndividual() - - assert opt.options["method"] == method - - # root method will error out, - # should revert to minimize method and optimize - opt(individual) - np.testing.assert_array_almost_equal( - fitness_function.evaluate_fitness_vector(individual), - np.array([1.0]) - ) - - # make sure method didn't stay as minimize method - assert opt.options["method"] == method diff --git a/tests/unit/symbolic_regression/test_custom_regression.py b/tests/unit/symbolic_regression/test_custom_regression.py new file mode 100644 index 00000000..e622aa0a --- /dev/null +++ b/tests/unit/symbolic_regression/test_custom_regression.py @@ -0,0 +1,168 @@ +# pylint: disable=missing-docstring +import numpy as np + +from bingo.expressions.agraph.evolvable import EvolvableExpression +from bingo.expressions.agraph.pyagraph.expression import AGraphExpression +from bingo.expressions.agraph.pyagraph.operators import CONSTANT, VARIABLE +from bingo.symbolic_regression import ( + CustomRegression, + FitResult, + ObjectiveData, + ScipyFitter, + explicit_residuals, + expression_loss, +) + + +def _constant_expression(value=0.0): + expression = AGraphExpression() + expression.raw_command_array = np.array([[CONSTANT, 0, 0]], dtype=np.uint8) + expression.raw_constants = (value,) + return expression + + +def test_custom_regression_commits_fitter_constants_then_ranks_expression(): + expression = _constant_expression() + data = ObjectiveData(np.arange(3.0), np.full(3, 2.0)) + + def fitter(fitted_expression, objective_data, fitting_measure): + assert fitted_expression is expression + assert objective_data is data + assert fitting_measure(fitted_expression, objective_data, [2.0]) == 0.0 + return FitResult([2.0], success=True) + + def measure(fitted_expression, objective_data, constants): + return constants[0] - objective_data.arrays[1][0] + + def loss(fitted_expression, objective_data): + return np.mean((fitted_expression.predict(objective_data.arrays[0]) - 2.0) ** 2) + + objective = CustomRegression(data, fitter, measure, loss) + + assert objective(EvolvableExpression(expression)) == 0.0 + assert expression.constants == (2.0,) + assert expression.is_fitted + + +def test_explicit_measure_evaluates_candidate_constants_without_mutating_expression(): + expression = _constant_expression(1.0) + data = ObjectiveData(np.arange(3.0).reshape(-1, 1), np.full(3, 2.0)) + + residuals = explicit_residuals()(expression, data, [2.0]) + + np.testing.assert_array_equal(residuals, np.zeros(3)) + assert expression.constants == (1.0,) + assert not expression.is_fitted + + +def test_scipy_root_fitter_uses_residual_measure_to_fit_constants(): + expression = _constant_expression() + data = ObjectiveData(np.arange(3.0).reshape(-1, 1), np.full(3, 2.0)) + objective = CustomRegression( + data, + ScipyFitter("lm"), + explicit_residuals(), + lambda fitted_expression, objective_data: fitted_expression.loss( + objective_data.arrays[0], objective_data.arrays[1] + ), + ) + + assert objective(EvolvableExpression(expression)) == 0.0 + assert expression.constants == (2.0,) + + +def test_scipy_minimize_fitter_uses_scalar_measure_to_fit_constants(): + expression = _constant_expression() + data = ObjectiveData(np.arange(3.0).reshape(-1, 1), np.full(3, 2.0)) + objective = CustomRegression( + data, + ScipyFitter("BFGS"), + expression_loss(), + lambda fitted_expression, objective_data: fitted_expression.loss( + objective_data.arrays[0], objective_data.arrays[1] + ), + ) + + assert objective(EvolvableExpression(expression)) < 1e-10 + assert expression.is_fitted + + +def test_invalid_fit_result_leaves_expression_unfitted(): + expression = _constant_expression() + data = ObjectiveData(np.arange(3.0), np.full(3, 2.0)) + objective = CustomRegression( + data, + lambda *_: FitResult([np.nan], success=False), + lambda *_: 0.0, + lambda *_: 0.0, + ) + + with np.testing.assert_raises_regex(ValueError, "fitted constants must be finite"): + objective(EvolvableExpression(expression)) + assert expression.constants == (0.0,) + assert not expression.is_fitted + + +def test_invalid_fit_result_shape_leaves_expression_unfitted(): + expression = _constant_expression() + objective = CustomRegression( + ObjectiveData(np.arange(3.0)), + lambda *_: FitResult([[2.0]], success=True), + lambda *_: 0.0, + lambda *_: 0.0, + ) + + with np.testing.assert_raises_regex(ValueError, "one-dimensional"): + objective(EvolvableExpression(expression)) + assert expression.constants == (0.0,) + assert not expression.is_fitted + + +def test_loss_exception_leaves_newly_fitted_expression_unfitted(): + expression = _constant_expression() + objective = CustomRegression( + ObjectiveData(np.arange(3.0)), + lambda *_: FitResult([2.0], success=True), + lambda *_: 0.0, + lambda *_: (_ for _ in ()).throw(RuntimeError("loss failed")), + ) + + with np.testing.assert_raises_regex(RuntimeError, "loss failed"): + objective(EvolvableExpression(expression)) + assert expression.constants == (0.0,) + assert not expression.is_fitted + + +def test_nonconverged_finite_fit_result_commits_best_constants(): + expression = _constant_expression() + objective = CustomRegression( + ObjectiveData(np.arange(3.0)), + lambda *_: FitResult([3.0], success=False, message="maximum iterations"), + lambda *_: 0.0, + lambda fitted_expression, _: fitted_expression.constants[0], + ) + + assert objective(EvolvableExpression(expression)) == 3.0 + assert expression.constants == (3.0,) + assert expression.is_fitted + + +def test_constant_free_expression_skips_fitter_and_ranks_directly(): + expression = AGraphExpression() + expression.raw_command_array = np.array([[VARIABLE, 0, 0]], dtype=np.uint8) + objective = CustomRegression( + ObjectiveData(np.arange(3.0).reshape(-1, 1)), + lambda *_: (_ for _ in ()).throw(AssertionError("fitter should not run")), + lambda *_: 0.0, + lambda fitted_expression, data: fitted_expression.predict(data.arrays[0]).mean(), + ) + + assert objective(EvolvableExpression(expression)) == 1.0 + + +def test_scipy_fitter_rejects_incompatible_measure_contract(): + expression = _constant_expression() + data = ObjectiveData(np.arange(3.0).reshape(-1, 1), np.full(3, 2.0)) + + with np.testing.assert_raises_regex(TypeError, "residual fitting measure"): + ScipyFitter("lm")(expression, data, expression_loss()) diff --git a/tests/unit/symbolic_regression/test_evidence.py b/tests/unit/symbolic_regression/test_evidence.py new file mode 100644 index 00000000..0a031812 --- /dev/null +++ b/tests/unit/symbolic_regression/test_evidence.py @@ -0,0 +1,266 @@ +# pylint: disable=missing-docstring +import builtins +import importlib.util +import multiprocessing + +import numpy as np +import pytest + +from bingo.expressions.agraph.pyagraph.expression import AGraphExpression +from bingo.expressions.agraph.pyagraph.operators import CONSTANT +from bingo.symbolic_regression import ( + EvidenceResult, + FitResult, + ObjectiveData, + ResidualMeasure, + SmcEvidenceEstimator, + explicit_residuals, +) +from bingo.symbolic_regression.evidence import smc_nmll_loss + + +requires_smcpy = pytest.mark.skipif( + importlib.util.find_spec("smcpy") is None, + reason="requires the optional evidence dependency", +) + +try: + from bingo.expressions.agraph.cppagraph import AGraphExpression as CppAGraphExpression +except ImportError: + CppAGraphExpression = None + + +def _constant_expression(value=0.0): + expression = AGraphExpression() + expression.raw_command_array = np.array([[CONSTANT, 0, 0]], dtype=np.uint8) + expression.raw_constants = (value,) + return expression + + +def _estimate_constant_evidence(_): + result = SmcEvidenceEstimator(num_particles=20, mcmc_steps=2, seed=21).estimate( + _constant_expression(1.0), + ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)), + explicit_residuals(), + ) + return result.smc_nmll, result.map_constants + + +@requires_smcpy +def test_estimate_returns_unsuccessful_result_when_proposal_cannot_be_built(): + expression = _constant_expression(1.0) + estimator = SmcEvidenceEstimator() + + result = estimator.estimate( + expression, + ObjectiveData(np.arange(3.0), np.full(3, 2.0)), + lambda *_: np.full(3, np.inf), + ) + + assert isinstance(result, EvidenceResult) + assert not result.success + assert result.smc_nmll == -np.inf + assert result.map_constants is None + assert expression.constants == (1.0,) + assert not expression.is_fitted + + +@requires_smcpy +def test_estimate_installs_only_posterior_map_constants_after_sampling(): + expression = _constant_expression(1.0) + estimator = SmcEvidenceEstimator(num_particles=20, mcmc_steps=2, seed=4) + + result = estimator.estimate( + expression, + ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)), + explicit_residuals(), + ) + + assert result.success + assert np.isfinite(result.smc_nmll) + assert result.map_constants == expression.constants + assert expression.is_fitted + assert result.posterior is None + + +@requires_smcpy +def test_failed_estimation_maps_to_infinite_ranking_loss(): + result = SmcEvidenceEstimator().estimate( + _constant_expression(), + ObjectiveData(np.arange(3.0), np.full(3, 2.0)), + lambda *_: np.full(3, np.inf), + ) + + assert smc_nmll_loss(result) == np.inf + + +@requires_smcpy +def test_estimate_retries_invalid_proposal_components(): + expression = _constant_expression(1.0) + calls = 0 + + def fitter(*_): + nonlocal calls + calls += 1 + return FitResult([np.nan] if calls == 1 else [2.0], success=True) + + result = SmcEvidenceEstimator( + num_particles=20, mcmc_steps=2, seed=4, fitter=fitter + ).estimate( + expression, + ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)), + explicit_residuals(), + ) + + assert result.success + assert calls == 2 + + +@requires_smcpy +def test_estimate_propagates_user_fitter_exceptions(): + expression = _constant_expression(1.0) + + def failing_fitter(*_): + raise RuntimeError("fitter failed") + + with pytest.raises(RuntimeError, match="fitter failed"): + SmcEvidenceEstimator(fitter=failing_fitter).estimate( + expression, + ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)), + explicit_residuals(), + ) + + assert expression.constants == (1.0,) + assert not expression.is_fitted + + +@requires_smcpy +def test_estimate_propagates_user_measure_exceptions(): + expression = _constant_expression(1.0) + + def failing_measure(*_): + raise RuntimeError("measure failed") + + with pytest.raises(RuntimeError, match="measure failed"): + SmcEvidenceEstimator().estimate( + expression, + ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)), + failing_measure, + ) + + assert expression.constants == (1.0,) + assert not expression.is_fitted + + +@requires_smcpy +def test_estimate_regularizes_indefinite_laplace_curvature(): + measure = ResidualMeasure( + lambda *_: np.ones(5), + jacobian=lambda *_: np.zeros((5, 1)), + residual_hessian=lambda *_: -np.ones((5, 1, 1)), + ) + result = SmcEvidenceEstimator( + num_particles=20, + mcmc_steps=2, + seed=4, + fitter=lambda *_: FitResult([2.0], success=True), + ).estimate( + _constant_expression(1.0), + ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)), + measure, + ) + + assert result.success + + +def test_missing_smcpy_returns_unsuccessful_result(mocker): + real_import = builtins.__import__ + + def unavailable(name, *args, **kwargs): + if name == "smcpy" or name.startswith("smcpy."): + raise ImportError("SMCPy unavailable") + return real_import(name, *args, **kwargs) + + mocker.patch("builtins.__import__", side_effect=unavailable) + result = SmcEvidenceEstimator().estimate( + _constant_expression(), + ObjectiveData(np.arange(3.0), np.full(3, 2.0)), + explicit_residuals(), + ) + + assert not result.success + assert result.smc_nmll == -np.inf + + +@requires_smcpy +def test_fixed_seed_reproduces_evidence_and_can_return_posterior(): + data = ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)) + first = _constant_expression(1.0) + second = _constant_expression(8.0) + estimator = SmcEvidenceEstimator( + num_particles=20, + mcmc_steps=2, + seed=11, + return_posterior=True, + fitter=lambda *_: FitResult([2.0], success=True), + ) + + first_result = estimator.estimate(first, data, explicit_residuals()) + second_result = estimator.estimate(second, data, explicit_residuals()) + + assert first_result.success + assert second_result.success + assert first_result.smc_nmll == second_result.smc_nmll + assert first_result.map_constants == second_result.map_constants + assert first_result.posterior is not None + assert not hasattr(first, "posterior") + + +def test_fixed_seed_ignores_raw_constants(): + first = _constant_expression(1.0) + second = _constant_expression(1.0) + second.raw_constants = (1.0, 2.0) + estimator = SmcEvidenceEstimator(seed=11) + + np.testing.assert_array_equal( + estimator._generator(first).integers(2**32, size=10), + estimator._generator(second).integers(2**32, size=10), + ) + + +@pytest.mark.skipif(CppAGraphExpression is None, reason="C++ expression unavailable") +@requires_smcpy +def test_fixed_seed_is_backend_neutral(): + data = ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)) + python_result = SmcEvidenceEstimator(num_particles=20, mcmc_steps=2, seed=8).estimate( + _constant_expression(1.0), data, explicit_residuals() + ) + cpp_result = SmcEvidenceEstimator(num_particles=20, mcmc_steps=2, seed=8).estimate( + CppAGraphExpression(equation="1.0"), data, explicit_residuals() + ) + + assert python_result.success + assert cpp_result.success + assert python_result.smc_nmll == cpp_result.smc_nmll + assert python_result.map_constants == cpp_result.map_constants + + +@requires_smcpy +def test_fixed_seed_reproduces_evidence_in_multiprocessing(): + with multiprocessing.get_context("spawn").Pool(2) as pool: + results = pool.map(_estimate_constant_evidence, range(2)) + + assert results[0] == results[1] + + +@requires_smcpy +def test_none_seed_permits_nonreproducible_sampling(): + data = ObjectiveData(np.arange(5.0).reshape(-1, 1), np.full(5, 2.0)) + first = SmcEvidenceEstimator(num_particles=20, mcmc_steps=2, seed=None).estimate( + _constant_expression(1.0), data, explicit_residuals() + ) + second = SmcEvidenceEstimator(num_particles=20, mcmc_steps=2, seed=None).estimate( + _constant_expression(1.0), data, explicit_residuals() + ) + + assert (first.smc_nmll, first.map_constants) != (second.smc_nmll, second.map_constants) diff --git a/tests/unit/symbolic_regression/test_explicit_regression.py b/tests/unit/symbolic_regression/test_explicit_regression.py index 1a39d01b..0ec11f0f 100644 --- a/tests/unit/symbolic_regression/test_explicit_regression.py +++ b/tests/unit/symbolic_regression/test_explicit_regression.py @@ -17,15 +17,18 @@ def _constant_expression(value=0.0): def test_explicit_objective_fits_an_unfitted_expression(mocker): expression = _constant_expression() individual = EvolvableExpression(expression) - objective = ExplicitRegression(np.arange(4.0).reshape(-1, 1), np.full(4, 2.0)) + x = np.arange(4.0).reshape(-1, 1) + y = np.full(4, 2.0) + objective = ExplicitRegression(x, y) fit = mocker.spy(expression, "fit") assert objective(individual) == 0.0 fit.assert_called_once() fit_x, fit_y = fit.call_args.args - np.testing.assert_array_equal(fit_x, np.arange(4.0).reshape(-1, 1)) - np.testing.assert_array_equal(fit_y, np.full(4, 2.0)) + np.testing.assert_array_equal(fit_x, x) + np.testing.assert_array_equal(fit_y, y) assert fit.call_args.kwargs == {"tolerance": 1e-5} + assert expression.constants == (2.0,) assert expression.is_fitted assert objective.eval_count == 1 diff --git a/tests/unit/symbolic_regression/test_implicit_regression.py b/tests/unit/symbolic_regression/test_implicit_regression.py index ac22e855..01ad6ae2 100644 --- a/tests/unit/symbolic_regression/test_implicit_regression.py +++ b/tests/unit/symbolic_regression/test_implicit_regression.py @@ -32,14 +32,13 @@ def test_implicit_objective_fits_and_returns_expression_loss(mocker): x = np.arange(4.0).reshape(-1, 1) dx_dt = np.ones((4, 1)) objective = ImplicitRegression(x, dx_dt) - fit = mocker.spy(expression, "fit_implicit") + fit_implicit = mocker.spy(expression, "fit_implicit") loss = objective(individual) - fit.assert_called_once() - fit_x, fit_dx_dt = fit.call_args.args + fit_implicit.assert_called_once() + fit_x, fit_dx_dt = fit_implicit.call_args.args np.testing.assert_array_equal(fit_x, x) np.testing.assert_array_equal(fit_dx_dt, dx_dt) - assert fit.call_args.kwargs == {"tolerance": 1e-5} assert loss == expression.implicit_loss(x, dx_dt) assert expression.is_fitted diff --git a/tests/unit/symbolic_regression/test_objective_data.py b/tests/unit/symbolic_regression/test_objective_data.py new file mode 100644 index 00000000..d99e78dc --- /dev/null +++ b/tests/unit/symbolic_regression/test_objective_data.py @@ -0,0 +1,32 @@ +import numpy as np +import pytest + +from bingo.symbolic_regression import ObjectiveData + + +def test_objective_data_exposes_aligned_arrays_and_indexes_them_together(): + data = ObjectiveData(np.arange(4), np.arange(8).reshape(4, 2)) + + assert isinstance(data.arrays, tuple) + np.testing.assert_array_equal(data.arrays[0], np.arange(4)) + np.testing.assert_array_equal(data.arrays[1], np.arange(8).reshape(4, 2)) + + subset = data[[3, 1]] + + assert isinstance(subset, ObjectiveData) + np.testing.assert_array_equal(subset.arrays[0], [3, 1]) + np.testing.assert_array_equal(subset.arrays[1], [[6, 7], [2, 3]]) + + +def test_objective_data_rejects_misaligned_arrays(): + with pytest.raises(ValueError, match="equal first-axis length"): + ObjectiveData(np.arange(2), np.arange(3)) + + +def test_objective_data_scalar_index_retains_the_sample_axis(): + data = ObjectiveData(np.arange(4), np.arange(8).reshape(4, 2)) + + subset = data[2] + + np.testing.assert_array_equal(subset.arrays[0], [2]) + np.testing.assert_array_equal(subset.arrays[1], [[4, 5]])