From 73ea9dc7beb621b0bacb2ee9caa0683a1ad2d0ef Mon Sep 17 00:00:00 2001 From: "G. Bomarito" Date: Wed, 18 Sep 2024 12:03:57 -0400 Subject: [PATCH] using warning that can be filtered, adding linear correction option to agraphs, adding negative nmll explicit regression option --- bingo/evaluation/fitness_function.py | 24 +++++- bingo/evaluation/gradient_mixin.py | 49 +++++++++--- bingo/symbolic_regression/agraph/agraph.py | 77 ++++++++++--------- .../explicit_regression.py | 49 +++++++++--- 4 files changed, 137 insertions(+), 62 deletions(-) diff --git a/bingo/evaluation/fitness_function.py b/bingo/evaluation/fitness_function.py index 33ab21fb..886bf7b6 100644 --- a/bingo/evaluation/fitness_function.py +++ b/bingo/evaluation/fitness_function.py @@ -3,27 +3,39 @@ This module defines the basis of fitness evaluation in bingo evolutionary analyses. """ + from abc import ABCMeta, abstractmethod import numpy as np # Fitness metric functions, outside of FitnessFunction for use in GradientMixin -def mean_absolute_error(vector): +def mean_absolute_error(vector, individual=None): """Calculate the mean absolute error of an error vector""" return np.mean(np.abs(vector)) -def root_mean_squared_error(vector): +def root_mean_squared_error(vector, individual=None): """Calculate the root mean squared error of an error vector""" return np.sqrt(np.mean(np.square(vector))) -def mean_squared_error(vector): +def mean_squared_error(vector, individual=None): """Calculate the mean squared error of an error vector""" return np.mean(np.square(vector)) +def negative_nmll_laplace(vector, individual): + """Calculate the nmll squared error of an error vector""" + 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 + + class FitnessFunction(metaclass=ABCMeta): """Fitness evaluation metric for individuals. @@ -42,6 +54,7 @@ class FitnessFunction(metaclass=ABCMeta): training_data : TrainingData (Optional) data that can be used in fitness evaluation """ + def __init__(self, training_data=None): self.eval_count = 0 self.training_data = training_data @@ -80,6 +93,7 @@ class VectorBasedFunction(FitnessFunction, metaclass=ABCMeta): 'mean absolute error'/'mae', 'mean squared error'/'mse', and 'root mean squared error'/'rmse' """ + def __init__(self, training_data=None, metric="mae"): super().__init__(training_data) @@ -89,6 +103,8 @@ def __init__(self, training_data=None, metric="mae"): 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 else: raise ValueError("Invalid metric for Fitness Function") @@ -110,7 +126,7 @@ def __call__(self, individual): fitness of the individual """ fitness_vector = self.evaluate_fitness_vector(individual) - return self._metric(fitness_vector) + return self._metric(fitness_vector, individual) @abstractmethod def evaluate_fitness_vector(self, individual): diff --git a/bingo/evaluation/gradient_mixin.py b/bingo/evaluation/gradient_mixin.py index b575725c..b636af21 100644 --- a/bingo/evaluation/gradient_mixin.py +++ b/bingo/evaluation/gradient_mixin.py @@ -4,11 +4,16 @@ 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 +from .fitness_function import ( + mean_absolute_error, + mean_squared_error, + root_mean_squared_error, + negative_nmll_laplace, +) class GradientMixin(metaclass=ABCMeta): @@ -17,6 +22,7 @@ class GradientMixin(metaclass=ABCMeta): 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 @@ -53,21 +59,28 @@ class VectorGradientMixin(GradientMixin): 'mean absolute error', 'mean squared error', and 'root mean squared error' """ + 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 = \ + 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 + 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 = \ + 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 + ) else: raise ValueError("Invalid metric for vector gradient mixin") @@ -90,10 +103,10 @@ def get_fitness_and_gradient(self, individual): 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), \ - self._metric_derivative(fitness_vector, jacobian.transpose()) + 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): @@ -133,5 +146,17 @@ def _mean_squared_error_derivative(fitness_vector, fitness_partials): @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) + 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) + dll = -0.5 * n / dmse + dnmll = (1 - b) * dll + return -dnmll diff --git a/bingo/symbolic_regression/agraph/agraph.py b/bingo/symbolic_regression/agraph/agraph.py index d0c1ff5c..13f868db 100644 --- a/bingo/symbolic_regression/agraph/agraph.py +++ b/bingo/symbolic_regression/agraph/agraph.py @@ -47,9 +47,11 @@ 15 hyperbolic cosine :math:`cosh(p1)` ======== ======================================= ================= """ + import logging import numpy as np from sympy.core import Expr +import warnings from .string_parsing import eq_string_to_command_array_and_constants from .string_generation import get_formatted_string @@ -63,6 +65,7 @@ from .evaluation_backend import evaluation_backend from .simplification_backend import simplification_backend + LOGGER = logging.getLogger(__name__) USING_PYTHON_SIMPLIFICATION = False @@ -113,6 +116,7 @@ class AGraph(Equation): constants : tuple of numeric numeric constants that are used in the equation """ + def __init__(self, use_simplification=False, equation=None): super().__init__() @@ -139,8 +143,9 @@ def _init_command_array_and_const(self, equation): self._needs_opt = False self._modified = False elif isinstance(equation, (Expr, str)): - command_array, constants = \ - eq_string_to_command_array_and_constants(str(equation)) + command_array, constants = eq_string_to_command_array_and_constants( + str(equation) + ) self.set_local_optimization_params(constants) if len(constants) > 0: @@ -195,11 +200,13 @@ def _notify_modification(self): def _update(self): if self._use_simplification: - self._simplified_command_array = \ - simplification_backend.simplify_stack(self._command_array) + self._simplified_command_array = simplification_backend.simplify_stack( + self._command_array + ) else: - self._simplified_command_array = \ - simplification_backend.reduce_stack(self._command_array) + self._simplified_command_array = simplification_backend.reduce_stack( + self._command_array + ) const_commands = self._simplified_command_array[:, 0] == CONSTANT num_const = np.count_nonzero(const_commands) @@ -207,11 +214,13 @@ def _update(self): self._simplified_command_array[const_commands, 2] = np.arange(num_const) optimization_aggression = 0 - if optimization_aggression == 0 \ - and num_const <= len(self._simplified_constants): + if optimization_aggression == 0 and num_const <= len( + self._simplified_constants + ): self._simplified_constants = self._simplified_constants[:num_const] - elif optimization_aggression == 1 \ - and num_const == len(self._simplified_constants): + elif optimization_aggression == 1 and num_const == len( + self._simplified_constants + ): self._simplified_constants = self._simplified_constants[:num_const] else: self._simplified_constants = (1.0,) * num_const @@ -282,8 +291,7 @@ def get_utilized_commands(self): list of bool of length N Boolean values for whether each command is utilized. """ - return simplification_backend.get_utilized_commands( - self._command_array) + return simplification_backend.get_utilized_commands(self._command_array) def evaluate_equation_at(self, x): """Evaluate the `AGraph` equation. @@ -304,13 +312,12 @@ def evaluate_equation_at(self, x): if self._modified: self._update() try: - f_of_x = \ - evaluation_backend.evaluate(self._simplified_command_array, x, - self._simplified_constants) + f_of_x = evaluation_backend.evaluate( + self._simplified_command_array, x, self._simplified_constants + ) return f_of_x - except (ArithmeticError, OverflowError, ValueError, - FloatingPointError) as err: - LOGGER.warning("%s in stack evaluation", err) + except (ArithmeticError, OverflowError, ValueError, FloatingPointError) as err: + warnings.warn(f"{err} in stack evaluation") return np.full(x.shape, np.nan) def evaluate_equation_with_x_gradient_at(self, x): @@ -333,12 +340,11 @@ def evaluate_equation_with_x_gradient_at(self, x): self._update() try: f_of_x, df_dx = evaluation_backend.evaluate_with_derivative( - self._simplified_command_array, x, - self._simplified_constants, True) + self._simplified_command_array, x, self._simplified_constants, True + ) return f_of_x, df_dx - except (ArithmeticError, OverflowError, ValueError, - FloatingPointError) as err: - LOGGER.warning("%s in stack evaluation/deriv", err) + except (ArithmeticError, OverflowError, ValueError, FloatingPointError) as err: + warnings.warn(f"{err} in stack evaluation/deriv") nan_array = np.full(x.shape, np.nan) return nan_array, np.array(nan_array) @@ -363,14 +369,12 @@ def evaluate_equation_with_local_opt_gradient_at(self, x): self._update() try: f_of_x, df_dc = evaluation_backend.evaluate_with_derivative( - self._simplified_command_array, x, - self._simplified_constants, False) + self._simplified_command_array, x, self._simplified_constants, False + ) return f_of_x, df_dc - except (ArithmeticError, OverflowError, ValueError, - FloatingPointError) as err: - LOGGER.warning("%s in stack evaluation/const-deriv", err) - nan_array = np.full((x.shape[0], len(self._simplified_constants)), - np.nan) + except (ArithmeticError, OverflowError, ValueError, FloatingPointError) as err: + warnings.warn(f"{err} in stack evaluation/const-deriv") + nan_array = np.full((x.shape[0], len(self._simplified_constants)), np.nan) return nan_array, np.array(nan_array) def __str__(self): @@ -404,8 +408,9 @@ def get_formatted_string(self, format_, raw=False): return get_formatted_string(format_, self._command_array, tuple()) if self._modified: self._update() - return get_formatted_string(format_, self._simplified_command_array, - self._simplified_constants) + return get_formatted_string( + format_, self._simplified_command_array, self._simplified_constants + ) def get_complexity(self): """Calculate complexity of agraph equation. @@ -448,10 +453,10 @@ def _copy_agraph_values_to_new_graph(self, agraph_duplicate): agraph_duplicate._fitness = self._fitness agraph_duplicate._fit_set = self._fit_set agraph_duplicate._command_array = np.copy(self.command_array) - agraph_duplicate._simplified_command_array = \ - np.copy(self._simplified_command_array) - agraph_duplicate._simplified_constants = \ - tuple(self._simplified_constants) + agraph_duplicate._simplified_command_array = np.copy( + self._simplified_command_array + ) + agraph_duplicate._simplified_constants = tuple(self._simplified_constants) agraph_duplicate._needs_opt = self._needs_opt agraph_duplicate._modified = self._modified agraph_duplicate._use_simplification = self._use_simplification diff --git a/bingo/symbolic_regression/explicit_regression.py b/bingo/symbolic_regression/explicit_regression.py index a76bd843..ce8cc297 100644 --- a/bingo/symbolic_regression/explicit_regression.py +++ b/bingo/symbolic_regression/explicit_regression.py @@ -7,8 +7,10 @@ that are unique to explicit symbolic regression. Namely, these classes are an appropriate fitness evaluator and a corresponding training data container. """ + import logging import numpy as np +from scipy.stats import linregress from ..evaluation.fitness_function import VectorBasedFunction from ..evaluation.gradient_mixin import VectorGradientMixin @@ -33,13 +35,19 @@ class ExplicitRegression(VectorGradientMixin, VectorBasedFunction): relative : bool Whether to use relative, pointwise normalization of errors. Default: False. + use_linear_correction : bool + Whether to adjust outputs of equations by a least squares linear correction. Default: False. """ - def __init__(self, training_data, metric="mae", relative=False): + + def __init__( + self, training_data, metric="mae", relative=False, use_linear_correction=False + ): super().__init__(training_data, metric) self._relative = relative + self._linear_correction = use_linear_correction def evaluate_fitness_vector(self, individual): - """ Traditional fitness evaluation for symbolic regression + """Traditional fitness evaluation for symbolic regression fitness = y - f(x) where x and y are in the training_data (i.e. training_data.x and training_data.y) and the function f is defined by @@ -57,6 +65,16 @@ def evaluate_fitness_vector(self, individual): """ self.eval_count += 1 f_of_x = individual.evaluate_equation_at(self.training_data.x) + + if self._linear_correction: + try: + slope, intercept, _, _, _ = linregress( + f_of_x.flatten(), self.training_data.y.flatten() + ) + f_of_x = intercept + slope * f_of_x + except ValueError: + pass + error = f_of_x - self.training_data.y if not self._relative: return np.squeeze(error) @@ -90,14 +108,24 @@ def get_fitness_vector_and_jacobian(self, individual): to the individual's constants """ self.eval_count += 1 - f_of_x, df_dc = \ - individual.evaluate_equation_with_local_opt_gradient_at( - self.training_data.x) + f_of_x, df_dc = individual.evaluate_equation_with_local_opt_gradient_at( + self.training_data.x + ) + + if self._linear_correction: + try: + slope, intercept, _, _, _ = linregress( + f_of_x.flatten(), self.training_data.y.flatten() + ) + f_of_x = intercept + slope * f_of_x + df_dc *= slope + except ValueError: + pass + error = f_of_x - self.training_data.y if not self._relative: return np.squeeze(error), df_dc - return np.squeeze(error / self.training_data.y), \ - df_dc / self.training_data.y + return np.squeeze(error / self.training_data.y), df_dc / self.training_data.y class ExplicitTrainingData(TrainingData): @@ -113,20 +141,21 @@ class ExplicitTrainingData(TrainingData): y : 2D numpy array dependent variable """ + def __init__(self, x, y): if x.ndim == 1: # warnings.warn("Explicit training x should be 2 dim array, " + # "reshaping array") x = x.reshape([-1, 1]) if x.ndim > 2: - raise TypeError('Explicit training x should be 2 dim array') + raise TypeError("Explicit training x should be 2 dim array") if y.ndim == 1: # warnings.warn("Explicit training y should be 2 dim array, " + # "reshaping array") y = y.reshape([-1, 1]) if y.ndim > 2: - raise TypeError('Explicit training y should be 2 dim array') + raise TypeError("Explicit training y should be 2 dim array") self._x = x self._y = y @@ -158,7 +187,7 @@ def __getitem__(self, items): return temp def __len__(self): - """ gets the length of the first dimension of the data + """gets the length of the first dimension of the data Returns -------