From 49a51a82b115be6344bd17e01cd2a200371eefca Mon Sep 17 00:00:00 2001 From: Maximilian Dinkel Date: Tue, 25 Feb 2025 09:19:15 +0100 Subject: [PATCH] feat: least-squares iterator --- queens/iterators/__init__.py | 4 +- queens/iterators/least_squares_iterator.py | 145 ++++++ queens/iterators/lm_iterator.py | 415 ----------------- queens/iterators/optimization_iterator.py | 41 +- .../python/test_lm_rosenbrock_res.py | 38 +- .../python/test_optimization_lsq_parabola.py | 6 +- .../test_optimization_lsq_rosenbrock.py | 10 +- .../test_optimization_lsq_rosenbrock_1d.py | 6 +- .../unit_tests/iterators/test_lm_iterator.py | 425 ------------------ 9 files changed, 181 insertions(+), 909 deletions(-) create mode 100644 queens/iterators/least_squares_iterator.py delete mode 100644 queens/iterators/lm_iterator.py delete mode 100644 tests/unit_tests/iterators/test_lm_iterator.py diff --git a/queens/iterators/__init__.py b/queens/iterators/__init__.py index 7d9f17769..b0cccf697 100644 --- a/queens/iterators/__init__.py +++ b/queens/iterators/__init__.py @@ -26,8 +26,8 @@ from queens.iterators.elementary_effects_iterator import ElementaryEffectsIterator from queens.iterators.grid_iterator import GridIterator from queens.iterators.hmc_iterator import HMCIterator +from queens.iterators.least_squares_iterator import LeastSquaresIterator from queens.iterators.lhs_iterator import LHSIterator -from queens.iterators.lm_iterator import LMIterator from queens.iterators.metropolis_hastings_iterator import MetropolisHastingsIterator from queens.iterators.metropolis_hastings_pymc_iterator import MetropolisHastingsPyMCIterator from queens.iterators.monte_carlo_iterator import MonteCarloIterator @@ -61,7 +61,7 @@ "points": PointsIterator, "bmfmc": BMFMCIterator, "grid": GridIterator, - "lm": LMIterator, + "least_squares": LeastSquaresIterator, "bbvi": BBVIIterator, "bmfia": BMFIAIterator, "rpvi": RPVIIterator, diff --git a/queens/iterators/least_squares_iterator.py b/queens/iterators/least_squares_iterator.py new file mode 100644 index 000000000..f7ea3dd5d --- /dev/null +++ b/queens/iterators/least_squares_iterator.py @@ -0,0 +1,145 @@ +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (c) 2024-2025, QUEENS contributors. +# +# This file is part of QUEENS. +# +# QUEENS is free software: you can redistribute it and/or modify it under the terms of the GNU +# Lesser General Public License as published by the Free Software Foundation, either version 3 of +# the License, or (at your option) any later version. QUEENS is distributed in the hope that it will +# be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You +# should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, +# see . +# +"""Least squares iterator.""" + +import logging + +import numpy as np +from scipy.optimize import Bounds, least_squares + +from queens.iterators.optimization_iterator import OptimizationIterator +from queens.utils.logger_settings import log_init_args + +_logger = logging.getLogger(__name__) + + +class LeastSquaresIterator(OptimizationIterator): + """Iterator for least-squares optimization. + + Based on the *scipy.optimize.least_squares* optimization toolbox [1]. + + References: + [1]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html + + Attributes: + algorithm (str): Algorithm to perform minimization: + + - trf : Trust Region Reflective algorithm, particularly suitable for large + sparse problems with bounds. Generally robust method. + + - dogbox : dogleg algorithm with rectangular trust regions, typical use + case is small problems with bounds. Not recommended for problems + with rank-deficient Jacobian. + - lm : Levenberg-Marquardt algorithm as implemented in MINPACK. Doesn’t + handle bounds and sparse Jacobians. Usually the most efficient + method for small unconstrained problems. + """ + + @log_init_args + def __init__( + self, + model, + parameters, + global_settings, + initial_guess, + result_description, + verbose_output=False, + bounds=Bounds(lb=-np.inf, ub=np.inf), + max_feval=None, + algorithm="lm", + jac_method="2-point", + jac_rel_step=None, + objective_and_jacobian=True, + ): + """Initialize LeastSquaresIterator. + + Args: + model (Model): Model to be evaluated by iterator + parameters (Parameters): Parameters object + global_settings (GlobalSettings): settings of the QUEENS experiment including its name + and the output directory + initial_guess (array like): initial position at which the optimization starts + result_description (dict): Description of desired post-processing. + verbose_output (int): Integer encoding which kind of verbose information should be + printed by the optimizers. + bounds (sequence, Bounds): Bounds on variables for Nelder-Mead, L-BFGS-B, TNC, SLSQP, + Powell, and trust-constr methods. + There are two ways to specify the bounds: + + 1. Instance of `Bounds` class. + 2. A sequence with 2 elements. The first element corresponds + to a sequence of lower bounds and the second element to + sequence of upper bounds. The length of each of the two + subsequences must be equal to the number of variables. + max_feval (int): Maximum number of function evaluations. + algorithm (str): Algorithm to perform minimization: + + - trf : Trust Region Reflective algorithm, particularly suitable for + large sparse problems with bounds. Generally robust method. + + - dogbox : dogleg algorithm with rectangular trust regions, typical use + case is small problems with bounds. Not recommended for + problems with rank-deficient Jacobian. + - lm : Levenberg-Marquardt algorithm as implemented in MINPACK. + Doesn’t handle bounds and sparse Jacobians. Usually the most + efficient method for small unconstrained problems. + jac_method (str): Method to calculate a finite difference based approximation of the + Jacobian matrix: + + - '2-point': a one-sided scheme by definition + - '3-point': more exact but needs twice as many function evaluations + jac_rel_step (array_like): Relative step size to use for finite difference approximation + of Jacobian matrix. If None (default) then it is selected + automatically. (see SciPy documentation for details) + objective_and_jacobian (bool, opt): If true, every time the objective is evaluated also + the jacobian is evaluated. This leads to improved + batching, but can lead to unnecessary evaluations of + the jacobian during line-search. + Default for is true. + """ + super().__init__( + model=model, + parameters=parameters, + global_settings=global_settings, + initial_guess=initial_guess, + result_description=result_description, + verbose_output=verbose_output, + bounds=bounds, + max_feval=max_feval, + algorithm=algorithm, + jac_method=jac_method, + jac_rel_step=jac_rel_step, + objective_and_jacobian=objective_and_jacobian, + ) + self.algorithm = algorithm # We don't want algorithm.upper() here + + def core_run(self): + """Core run of Optimization iterator.""" + self.solution = least_squares( + self.objective, + self.initial_guess, + method=self.algorithm, + jac=self.jacobian, + bounds=self.bounds, + max_nfev=self.max_feval, + verbose=int(self.verbose_output), + ) + + def post_run(self): + """Analyze the resulting optimum.""" + _logger.info("Optimality:\n\t%s", self.solution.optimality) + _logger.info("Cost:\n\t%s", self.solution.cost) + + super().post_run() diff --git a/queens/iterators/lm_iterator.py b/queens/iterators/lm_iterator.py deleted file mode 100644 index 9714e1dc1..000000000 --- a/queens/iterators/lm_iterator.py +++ /dev/null @@ -1,415 +0,0 @@ -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# Copyright (c) 2024-2025, QUEENS contributors. -# -# This file is part of QUEENS. -# -# QUEENS is free software: you can redistribute it and/or modify it under the terms of the GNU -# Lesser General Public License as published by the Free Software Foundation, either version 3 of -# the License, or (at your option) any later version. QUEENS is distributed in the hope that it will -# be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You -# should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, -# see . -# -"""Levenberg Marquardt iterator.""" - -import logging - -import numpy as np -import pandas as pd -import plotly.express as px - -from queens.iterators.iterator import Iterator -from queens.utils.fd_jacobian import fd_jacobian -from queens.utils.logger_settings import log_init_args - -_logger = logging.getLogger(__name__) - - -class LMIterator(Iterator): - """Iterator for Levenberg-Marquardt deterministic optimization problems. - - Implements the Levenberg-Marquardt algorithm for optimization, adapted from the - 4C `gen_inv_analysis` approach. This class focuses on a simplified yet controlled - optimization process. - - Attributes: - initial_guess (np.ndarray): Initial guess for the optimization parameters. - bounds (tuple): Tuple specifying the lower and upper bounds for parameters. If None, no - bounds are applied. - havebounds (bool): Flag indicating if bounds are provided. - param_current (np.ndarray): Current parameter values being optimized. - jac_rel_step (float): Relative step size used for finite difference approximation of the - Jacobian. - max_feval (int): Maximum number of allowed function evaluations. - result_description (dict): Configuration for result file handling and plotting. - jac_abs_step (float): Absolute step size used for finite difference approximation of the - Jacobian. - reg_param (float): Regularization parameter for the Levenberg-Marquardt algorithm. - init_reg (float): Initial value for the regularization parameter. - update_reg (str): Strategy for updating the regularization parameter ("res" for - residual-based, "grad" for gradient-based). - tolerance (float): Convergence tolerance for the optimization process. - verbose_output (bool): If True, provides detailed output during optimization. - iter_opt (int): The iteration number at which the lowest error was achieved. - lowesterror (float or None): The minimum error encountered during the optimization process. - param_opt (np.ndarray): Parameter values corresponding to the minimum error. - solution (np.ndarray): The final optimized parameter values. - """ - - @log_init_args - def __init__( - self, - model, - parameters, - global_settings, - result_description, - initial_guess=None, - bounds=None, - jac_rel_step=1e-4, - jac_abs_step=0.0, - init_reg=1.0, - update_reg="grad", - convergence_tolerance=1e-6, - max_feval=1, - verbose_output=False, - ): - """Initializes the Levenberg-Marquardt iterator. - - Args: - model (Model): Model to be optimized. - parameters (Parameters): Object containing parameter settings. - global_settings (GlobalSettings): Configuration settings for the QUEENS experiment. - result_description (dict): Description and configuration for result handling and - plotting. - initial_guess (np.ndarray, optional): Initial guess for optimization parameters. - bounds (tuple, optional): Tuple of lower and upper bounds for parameters. If None, no - bounds are used. - jac_rel_step (float, optional): Relative step size for finite difference approximation. - jac_abs_step (float, optional): Absolute step size for finite difference approximation. - init_reg (float, optional): Initial regularization parameter. - update_reg (str, optional): Strategy for updating the regularization parameter ("res" or - "grad"). - convergence_tolerance (float, optional): Convergence tolerance for the optimization. - max_feval (int, optional): Maximum number of function evaluations allowed. - verbose_output (bool, optional): If True, enables verbose output during optimization. - """ - super().__init__(model, parameters, global_settings) - - self.havebounds = True - if bounds is None: - self.havebounds = False - - self.initial_guess = np.array(initial_guess, dtype=float) - self.bounds = bounds - self.param_current = self.initial_guess - self.jac_rel_step = jac_rel_step - self.max_feval = max_feval - self.result_description = result_description - - self.jac_abs_step = jac_abs_step - self.reg_param = init_reg - self.init_reg = init_reg - self.update_reg = update_reg - self.tolerance = convergence_tolerance - - self.verbose_output = verbose_output - self.iter_opt = 0 - self.lowesterror = None - self.param_opt = None - self.solution = None - - def jacobian_and_residual(self, x0): - """Evaluate Jacobian and residual of objective function at *x0*. - - For LM we can restrict to "2-point". - - Args: - x0 (numpy.ndarray): Vector with current parameters - - Returns: - jacobian_matrix (numpy.ndarray): Jacobian Matrix approximation from finite differences - f0 (numpy.ndarray): Residual of objective function at `x0` - """ - positions, delta_positions = self.get_positions_raw_2pointperturb(x0) - - self.model.evaluate(positions) - - f = self.model.response["result"] - f_batch = f[-(len(positions)) :] - - f0 = f_batch[0] # first entry corresponds to f(x0) - f_perturbed = np.delete(f_batch, 0, 0) # delete the first entry - jacobian_matrix = fd_jacobian(f0, f_perturbed, delta_positions, True, "2-point") - # sanity checks: - # in the case of LSQ, the number of residuals needs to be - # greater or equal to the number of parameters to be fitted - num_res, num_par = jacobian_matrix.shape - if num_res < num_par: - raise ValueError( - f"Number of residuals (={num_res}) has to be greater or equal to" - f" number of parameters (={num_par})." - f" You have {num_res}<{num_par}." - ) - - return jacobian_matrix, f0 - - def pre_run(self): - """Initialize run. - - Print console output and optionally open .csv file for results - and write header. - """ - _logger.info("Initialize Levenberg-Marquardt run.") - - # produce .csv file and write header - if self.result_description: - if self.result_description["write_results"]: - df = pd.DataFrame( - columns=["iter", "resnorm", "gradnorm", "params", "delta_params", "mu"], - ) - csv_file = self.global_settings.result_file(".csv") - df.to_csv(csv_file, mode="w", sep="\t", index=None) - - def core_run(self): - """Core run of Levenberg Marquardt iterator.""" - resnorm = np.inf - gradnorm = np.inf - - converged = False - - i = 0 - - # Levenberg Marquardt iterations - while not converged: - if i > self.max_feval: - converged = True - _logger.info("Maximum number of steps max_feval= %d reached.", self.max_feval) - break - - _logger.info( - "iteration: %d reg_param: %s current_parameters: %s", - i, - self.reg_param, - self.param_current, - ) - jacobian, residual = self.jacobian_and_residual(self.param_current) - - # store terms for repeated usage - jacobian_inner_product = (jacobian.T).dot(jacobian) - jacobian_transpose_residual_product = (jacobian.T).dot(residual) - - resnorm_o = resnorm - resnorm = np.linalg.norm(residual) / np.sqrt(residual.size) - - gradnorm_o = gradnorm - gradnorm = np.linalg.norm(jacobian_transpose_residual_product) - - if i == 0: - resnorm_o = resnorm - gradnorm_o = gradnorm - self.lowesterror = resnorm - self.param_opt = self.param_current - - # save most optimal step in resnorm - if resnorm < self.lowesterror: - self.lowesterror = resnorm - self.param_opt = self.param_current - self.iter_opt = i - - # compute update step. Not necessary for last iteration! Should usually be low - # dimension matrix inversion. We need J and r anyway for convergence check. - param_delta = np.linalg.solve( - jacobian_inner_product + self.reg_param * np.diag(np.diag(jacobian_inner_product)), - -(jacobian_transpose_residual_product), - ) - - # output for iteration. Before update including next step - self.printstep(i, resnorm, gradnorm, param_delta) - - # evaluate bounds - if self.havebounds and self.checkbounds(param_delta, i): - if self.reg_param > 1.0e6 * self.init_reg: - _logger.info( - "WARNING: STEP #%d IS OUT OF BOUNDS and reg_param is unreasonably " - "high. Ending iterations!", - i, - ) - break - continue - - # update param for next step - self.param_current = self.param_current + param_delta - - # update reg_param and check for tolerance - if self.update_reg == "res": - if resnorm < self.tolerance: - converged = True - break - self.reg_param = self.reg_param * resnorm / resnorm_o - elif self.update_reg == "grad": - if gradnorm < self.tolerance: - converged = True - break - if resnorm < resnorm_o and gradnorm < gradnorm_o: - self.reg_param = self.reg_param * gradnorm / gradnorm_o - else: - raise ValueError("update_reg unknown") - i += 1 - - # store set of parameters which leads to the lowest residual as solution - self.solution = self.param_opt - - def post_run(self): - """Post run. - - Write solution to the console and optionally create .html plot - from result file. - """ - _logger.info("The optimum:\t%s occurred in iteration #%s.", self.solution, self.iter_opt) - if self.result_description: - if self.result_description["plot_results"] and self.result_description["write_results"]: - data = pd.read_csv( - self.global_settings.result_file(".csv"), - sep="\t", - ) - xydata = data["params"] - xydata = xydata.str.extractall(r"([+-]?\d+\.\d*e?[+-]?\d*)") - xydata = xydata.unstack() - data = data.drop(columns="params") - i = 0 - for column in xydata: - data[self.parameters.names[i]] = xydata[column].astype(float) - i = i + 1 - - if i > 2: - _logger.warning( - "write_results for more than 2 parameters not implemented, " - "because we are limited to 3 dimensions. " - "You have: %d. Plotting is skipped.", - i, - ) - return - if i == 2: - fig = px.line_3d( - data, - x=self.parameters.names[0], - y=self.parameters.names[1], - z="resnorm", - hover_data=[ - "iter", - "resnorm", - "gradnorm", - "delta_params", - "mu", - self.parameters.names[0], - self.parameters.names[1], - ], - ) - fig.update_traces(mode="lines+markers", marker={"size": 2}, line={"width": 4}) - elif i == 1: - fig = px.line( - data, - x=self.parameters.names[0], - y="resnorm", - hover_data=[ - "iter", - "resnorm", - "gradnorm", - "delta_params", - "mu", - self.parameters.names[0], - ], - ) - fig.update_traces(mode="lines+markers", marker={"size": 7}, line={"width": 3}) - else: - raise ValueError("You shouldn't be here without parameters.") - - fig.write_html(self.global_settings.result_file(".html")) - - def get_positions_raw_2pointperturb(self, x0): - """Get parameter sets for objective function evaluations. - - Args: - x0 (numpy.ndarray): Vector with current parameters - - Returns: - positions (numpy.ndarray): Parameter batch for function evaluation - delta_positions (numpy.ndarray): Parameter perturbations for finite difference - scheme - """ - delta_positions = self.jac_abs_step + self.jac_rel_step * np.abs(x0) - # if bounds do not allow finite difference step, use opposite direction - if self.havebounds: - perturbed = x0 + delta_positions - for i, current_p in enumerate(perturbed): - if (current_p < self.bounds[0][i]) or (current_p > self.bounds[1][i]): - delta_positions[i] = -delta_positions[i] - - positions = np.zeros((x0.size + 1, x0.size)) - positions[0] = x0 - - for i in range(x0.size): - positions[i + 1] = x0 - positions[i + 1][i] = positions[i + 1][i] + delta_positions[i] - - delta_positions.shape = (delta_positions.size, 1) - - return positions, delta_positions - - def printstep(self, i, resnorm, gradnorm, param_delta): - """Print iteration data to console and optionally to file. - - Opens file in append mode, so that file is updated frequently. - - Args: - i (int): Iteration number - resnorm (float): Residual norm - gradnorm (float): Gradient norm - param_delta (numpy.ndarray): Parameter step - """ - # write iteration to file - if self.result_description: - if self.result_description["write_results"]: - df = pd.DataFrame( - { - "iter": i, - "resnorm": np.format_float_scientific(resnorm, precision=8), - "gradnorm": np.format_float_scientific(gradnorm, precision=8), - "params": [np.array2string(self.param_current, precision=8)], - "delta_params": [np.array2string(param_delta, precision=8)], - "mu": np.format_float_scientific(self.reg_param, precision=8), - } - ) - csv_file = self.global_settings.result_file(".csv") - df.to_csv( - csv_file, mode="a", sep="\t", header=None, index=None, float_format="%.8f" - ) - - def checkbounds(self, param_delta, i): - """Check if proposed step is in bounds. - - Otherwise double regularization and compute new step. - - Args: - param_delta (numpy.ndarray): Parameter step - i (int): Iteration number - - Returns: - stepisoutside (bool): Flag if proposed step is out of bounds - """ - stepisoutside = False - nextstep = self.param_current + param_delta - if np.any(nextstep < self.bounds[0]) or np.any(nextstep > self.bounds[1]): - stepisoutside = True - _logger.warning( - "WARNING: STEP #%d IS OUT OF BOUNDS; double reg_param and compute new iteration." - "\n declined step was: %s", - i, - nextstep, - ) - - self.reg_param *= 2 - - return stepisoutside diff --git a/queens/iterators/optimization_iterator.py b/queens/iterators/optimization_iterator.py index adcff5241..02aba9b30 100644 --- a/queens/iterators/optimization_iterator.py +++ b/queens/iterators/optimization_iterator.py @@ -18,7 +18,7 @@ import time import numpy as np -from scipy.optimize import Bounds, least_squares, minimize +from scipy.optimize import Bounds, minimize from scipy.optimize._numdiff import _prepare_bounds from queens.iterators.iterator import Iterator @@ -32,10 +32,10 @@ class OptimizationIterator(Iterator): """Iterator for deterministic optimization problems. - Based on the *scipy.optimize* optimization toolbox [1]. + Based on the *scipy.optimize.minimize* optimization toolbox [1]. References: - [1]: https://docs.scipy.org/doc/scipy/reference/optimize.html + [1]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html Attributes: algorithm (str): String that defines the optimization algorithm to be used: @@ -51,7 +51,6 @@ class OptimizationIterator(Iterator): Jacobian necessary - SLSQP: Sequential Least Squares Programming minimization with bounds and constraints using Jacobian - - LSQ: Nonlinear least squares with bounds using Jacobian - COBYLA: Constrained Optimization BY Linear Approximation (no Jacobian) - NELDER-MEAD: Downhill-simplex search method @@ -92,12 +91,8 @@ class OptimizationIterator(Iterator): objective_and_jacobian (bool): If true, every time the objective is evaluated also the jacobian is evaluated. This leads to improved batching, but can lead to unnecessary evaluations of the jacobian during - line-search. - This option is only available for gradient methods. Default - for 'LSQ' is true, and false for remaining gradient methods. - - Returns: - OptimizationIterator (obj): Instance of the OptimizationIterator + line-search. This option is only available for gradient + methods. Default is false. """ @log_init_args @@ -115,7 +110,7 @@ def __init__( algorithm="L-BFGS-B", jac_method="2-point", jac_rel_step=None, - objective_and_jacobian=None, + objective_and_jacobian=False, ): """Initialize an OptimizationIterator. @@ -176,8 +171,7 @@ def __init__( batching, but can lead to unnecessary evaluations of the jacobian during line-search. This option is only available for gradient methods. - Default for 'LSQ' is true, and false for remaining - gradient methods. + Default is false. """ super().__init__(model, parameters, global_settings) @@ -239,11 +233,6 @@ def __init__( self.objective_and_jacobian = objective_and_jacobian if self.algorithm in ["COBYLA", "NELDER-MEAD", "POWELL"]: self.objective_and_jacobian = False - if self.objective_and_jacobian is None: - if self.algorithm == "LSQ": - self.objective_and_jacobian = True - else: - self.objective_and_jacobian = False def objective(self, x0): """Evaluate objective function at *x0*. @@ -325,18 +314,9 @@ def core_run(self): """Core run of Optimization iterator.""" _logger.info("Welcome to Optimization core run.") start = time.time() - # nonlinear least squares with bounds using Jacobian - if self.algorithm == "LSQ": - self.solution = least_squares( - self.objective, - self.initial_guess, - jac=self.jacobian, - bounds=self.bounds, - max_nfev=self.max_feval, - verbose=int(self.verbose_output), - ) + # minimization with bounds using Jacobian - elif self.algorithm in {"L-BFGS-B", "TNC"}: + if self.algorithm in {"L-BFGS-B", "TNC"}: self.solution = minimize( self.objective, self.initial_guess, @@ -390,9 +370,6 @@ def core_run(self): def post_run(self): """Analyze the resulting optimum.""" _logger.info("The optimum:\n\t%s", self.solution.x) - if self.algorithm == "LSQ": - _logger.info("Optimality:\n\t%s", self.solution.optimality) - _logger.info("Cost:\n\t%s", self.solution.cost) if self.result_description: if self.result_description["write_results"]: diff --git a/tests/integration_tests/python/test_lm_rosenbrock_res.py b/tests/integration_tests/python/test_lm_rosenbrock_res.py index c9c550e1b..1919248f5 100644 --- a/tests/integration_tests/python/test_lm_rosenbrock_res.py +++ b/tests/integration_tests/python/test_lm_rosenbrock_res.py @@ -15,15 +15,15 @@ """Integration test for the Levenberg Marquardt iterator.""" import numpy as np -import pandas as pd from queens.distributions.free import FreeVariable from queens.drivers.function_driver import FunctionDriver -from queens.iterators.lm_iterator import LMIterator +from queens.iterators.least_squares_iterator import LeastSquaresIterator from queens.main import run_iterator from queens.models.simulation_model import SimulationModel from queens.parameters.parameters import Parameters from queens.schedulers.pool_scheduler import PoolScheduler +from queens.utils.io_utils import load_result def test_lm_rosenbrock_res(global_settings): @@ -37,35 +37,25 @@ def test_lm_rosenbrock_res(global_settings): driver = FunctionDriver(parameters=parameters, function="rosenbrock60_residual") scheduler = PoolScheduler(experiment_name=global_settings.experiment_name) model = SimulationModel(scheduler=scheduler, driver=driver) - iterator = LMIterator( - jac_rel_step=1e-05, - jac_abs_step=0.001, - max_feval=99, - init_reg=0.01, - update_reg="grad", - convergence_tolerance=1e-06, - initial_guess=[0.9, 0.9], - result_description={"write_results": True, "plot_results": True}, + iterator = LeastSquaresIterator( model=model, parameters=parameters, global_settings=global_settings, + initial_guess=[0.9, 0.9], + result_description={"write_results": True, "plot_results": True}, + verbose_output=False, + max_feval=99, + algorithm="lm", + jac_method="2-point", + jac_rel_step=1e-05, + objective_and_jacobian=True, ) # Actual analysis run_iterator(iterator, global_settings=global_settings) # Load results - result_file = global_settings.result_file(".csv") - data = pd.read_csv( - result_file, - sep="\t", - ) - - params = data.get("params").tail(1) - dfparams = params.str.extractall(r"([+-]?\d+\.\d*e?[+-]?\d*)") - dfparams = dfparams.astype(float) - numpyparams = dfparams.to_numpy() - - np.testing.assert_allclose(numpyparams, np.array([[+1.0], [+1.0]]), rtol=1.0e-5) + results = load_result(global_settings.result_file(".pickle")) - assert global_settings.result_file(".html").is_file() + np.testing.assert_allclose(results.x, np.array([+1.0, +1.0]), rtol=1.0e-5) + np.testing.assert_allclose(results.fun, np.array([+0.0, +0.0])) diff --git a/tests/integration_tests/python/test_optimization_lsq_parabola.py b/tests/integration_tests/python/test_optimization_lsq_parabola.py index 4aa13baf8..67fe74ff6 100644 --- a/tests/integration_tests/python/test_optimization_lsq_parabola.py +++ b/tests/integration_tests/python/test_optimization_lsq_parabola.py @@ -21,7 +21,7 @@ from queens.distributions.free import FreeVariable from queens.drivers.function_driver import FunctionDriver -from queens.iterators.optimization_iterator import OptimizationIterator +from queens.iterators.least_squares_iterator import LeastSquaresIterator from queens.main import run_iterator from queens.models.simulation_model import SimulationModel from queens.parameters.parameters import Parameters @@ -42,8 +42,8 @@ def test_optimization_lsq_parabola(global_settings): driver = FunctionDriver(parameters=parameters, function="parabola_residual") scheduler = PoolScheduler(experiment_name=global_settings.experiment_name) model = SimulationModel(scheduler=scheduler, driver=driver) - iterator = OptimizationIterator( - algorithm="LSQ", + iterator = LeastSquaresIterator( + algorithm="lm", initial_guess=[0.75], result_description={"write_results": True}, bounds=[float("-inf"), float("inf")], diff --git a/tests/integration_tests/python/test_optimization_lsq_rosenbrock.py b/tests/integration_tests/python/test_optimization_lsq_rosenbrock.py index f127ad0a0..871fef6b2 100644 --- a/tests/integration_tests/python/test_optimization_lsq_rosenbrock.py +++ b/tests/integration_tests/python/test_optimization_lsq_rosenbrock.py @@ -21,7 +21,7 @@ from queens.distributions.free import FreeVariable from queens.drivers.function_driver import FunctionDriver -from queens.iterators.optimization_iterator import OptimizationIterator +from queens.iterators.least_squares_iterator import LeastSquaresIterator from queens.main import run_iterator from queens.models.simulation_model import SimulationModel from queens.parameters.parameters import Parameters @@ -40,8 +40,8 @@ def test_optimization_lsq_rosenbrock(global_settings): driver = FunctionDriver(parameters=parameters, function="rosenbrock60_residual") scheduler = PoolScheduler(experiment_name=global_settings.experiment_name) model = SimulationModel(scheduler=scheduler, driver=driver) - iterator = OptimizationIterator( - algorithm="LSQ", + iterator = LeastSquaresIterator( + algorithm="dogbox", initial_guess=[-3.0, -4.0], result_description={"write_results": True}, bounds=[float("-inf"), float("inf")], @@ -72,8 +72,8 @@ def test_optimization_lsq_rosenbrock_error(global_settings): driver = FunctionDriver(parameters=parameters, function="rosenbrock60_residual_3d") scheduler = PoolScheduler(experiment_name=global_settings.experiment_name) model = SimulationModel(scheduler=scheduler, driver=driver) - iterator = OptimizationIterator( - algorithm="LSQ", + iterator = LeastSquaresIterator( + algorithm="trf", initial_guess=[-3.0, -4.0, -5.0], result_description={"write_results": True}, bounds=[float("-inf"), float("inf")], diff --git a/tests/integration_tests/python/test_optimization_lsq_rosenbrock_1d.py b/tests/integration_tests/python/test_optimization_lsq_rosenbrock_1d.py index 201ebda93..ebf8adf33 100644 --- a/tests/integration_tests/python/test_optimization_lsq_rosenbrock_1d.py +++ b/tests/integration_tests/python/test_optimization_lsq_rosenbrock_1d.py @@ -21,7 +21,7 @@ from queens.distributions.free import FreeVariable from queens.drivers.function_driver import FunctionDriver -from queens.iterators.optimization_iterator import OptimizationIterator +from queens.iterators.least_squares_iterator import LeastSquaresIterator from queens.main import run_iterator from queens.models.simulation_model import SimulationModel from queens.parameters.parameters import Parameters @@ -42,8 +42,8 @@ def test_optimization_lsq_rosenbrock_1d(global_settings): driver = FunctionDriver(parameters=parameters, function="rosenbrock60_residual_1d") scheduler = PoolScheduler(experiment_name=global_settings.experiment_name) model = SimulationModel(scheduler=scheduler, driver=driver) - iterator = OptimizationIterator( - algorithm="LSQ", + iterator = LeastSquaresIterator( + algorithm="trf", initial_guess=[3.0], result_description={"write_results": True}, bounds=[float("-inf"), float("inf")], diff --git a/tests/unit_tests/iterators/test_lm_iterator.py b/tests/unit_tests/iterators/test_lm_iterator.py deleted file mode 100644 index 91cc2a4bf..000000000 --- a/tests/unit_tests/iterators/test_lm_iterator.py +++ /dev/null @@ -1,425 +0,0 @@ -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# Copyright (c) 2024-2025, QUEENS contributors. -# -# This file is part of QUEENS. -# -# QUEENS is free software: you can redistribute it and/or modify it under the terms of the GNU -# Lesser General Public License as published by the Free Software Foundation, either version 3 of -# the License, or (at your option) any later version. QUEENS is distributed in the hope that it will -# be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You -# should have received a copy of the GNU Lesser General Public License along with QUEENS. If not, -# see . -# -"""Test for LM iterator.""" - -import logging - -import numpy as np -import pandas as pd -import plotly.express as px -import pytest -from mock import Mock - -from queens.distributions.free import FreeVariable -from queens.iterators.lm_iterator import LMIterator -from queens.models.simulation_model import SimulationModel -from queens.parameters.parameters import Parameters - -_logger = logging.getLogger(__name__) - - -@pytest.fixture(name="iterator_name_cases", scope="module", params=["method"]) -def fixture_iterator_name_cases(request): - """Iterator name cases.""" - return request.param - - -@pytest.fixture(name="model_cases", scope="module", params=[None, "dummy_model"]) -def fixture_model_cases(request): - """Model cases.""" - return request.param - - -@pytest.fixture(name="fix_update_reg", scope="module", params=["grad", "res", "not_valid"]) -def fixture_fix_update_reg(request): - """Update reg cases.""" - return request.param - - -@pytest.fixture(name="fix_tolerance", scope="module", params=[1e-6, 1e0]) -def fixture_fix_tolerance(request): - """Tolerance cases.""" - return request.param - - -@pytest.fixture(name="output_csv") -def fixture_output_csv(global_settings): - """The absolute path to output csv file.""" - return global_settings.result_file(".csv") - - -@pytest.fixture(name="output_html") -def fixture_output_html(global_settings): - """The absolute path to output html file.""" - return global_settings.result_file(".html") - - -@pytest.fixture(name="default_lm_iterator") -def fixture_default_lm_iterator(global_settings): - """A default LMIterator instance.""" - parameters = Parameters(x1=FreeVariable(1), x2=FreeVariable(1)) - model = SimulationModel(scheduler=Mock(), driver=Mock()) - - my_lm_iterator = LMIterator( - model=model, - parameters=parameters, - global_settings=global_settings, - result_description={"write_results": True, "plot_results": True}, - initial_guess=[0.1, 0.2], - jac_rel_step=1e-05, - jac_abs_step=0.001, - init_reg=1.0, - update_reg="grad", - convergence_tolerance=1e-06, - max_feval=99, - ) - - return my_lm_iterator - - -@pytest.fixture(name="fix_true_false_param", scope="module", params=[True, False]) -def fixture_fix_true_false_param(request): - """A boolean parameter.""" - return request.param - - -@pytest.fixture(name="fix_plotly_fig", scope="module") -def fixture_fix_plotly_fig(): - """A Plotly figure.""" - data = pd.DataFrame({"x": [1.0, 2.0], "y": [1.1, 2.1], "z": [1.2, 2.2]}) - fig = px.line_3d( - data, - x="x", - y="y", - z="z", - ) - return fig - - -def test_init(global_settings): - """Test LMIterator initialization.""" - initial_guess = np.array([1, 2.2]) - bounds = np.array([[0.0, 1.0], [1.0, 2.0]]) - jac_rel_step = 1e-3 - jac_abs_step = 1e-2 - init_reg = 1.0 - update_reg = "grad" - tolerance = 1e-8 - max_feval = 99 - model = "dummy_model" - result_description = (True,) - verbose_output = (True,) - - my_lm_iterator = LMIterator( - model=model, - parameters="dummy_parameters", - global_settings=global_settings, - result_description=result_description, - initial_guess=initial_guess, - bounds=bounds, - jac_rel_step=jac_rel_step, - jac_abs_step=jac_abs_step, - init_reg=init_reg, - update_reg=update_reg, - convergence_tolerance=tolerance, - max_feval=max_feval, - verbose_output=verbose_output, - ) - - np.testing.assert_equal(my_lm_iterator.initial_guess, initial_guess) - np.testing.assert_equal(my_lm_iterator.param_current, initial_guess) - assert my_lm_iterator.model == model - assert my_lm_iterator.jac_rel_step == jac_rel_step - assert my_lm_iterator.max_feval == max_feval - assert my_lm_iterator.result_description == result_description - assert my_lm_iterator.jac_abs_step == jac_abs_step - assert my_lm_iterator.reg_param == init_reg - assert my_lm_iterator.update_reg == update_reg - assert my_lm_iterator.verbose_output == verbose_output - - -def test_model_evaluate(default_lm_iterator, mocker): - """Test model evaluation in LMIterator.""" - mp = mocker.patch("queens.models.simulation_model.SimulationModel.evaluate", return_value=None) - default_lm_iterator.model.evaluate(None) - mp.assert_called_once() - - -def test_residual(default_lm_iterator, mocker): - """Test residual calculation in LMIterator.""" - mocker.patch( - "queens.iterators.lm_iterator.LMIterator.get_positions_raw_2pointperturb", - return_value=[np.array([[1.0, 2.2], [1.00101, 2.2], [1.0, 2.201022]]), 1], - ) - - m2 = mocker.patch( - "queens.models.simulation_model.SimulationModel.evaluate", - return_value=None, - ) - - default_lm_iterator.model.response = {"result": np.array([[3.0, 4.2], [99.9, 99.9]])} - - _, result = default_lm_iterator.jacobian_and_residual(np.array([1.0, 2.2])) - - np.testing.assert_equal(result, np.array([3.0, 4.2])) - m2.assert_called_once() - - -def test_jacobian(default_lm_iterator, fix_true_false_param, mocker): - """Test Jacobian calculation in LMIterator.""" - mocker.patch( - "queens.iterators.lm_iterator.LMIterator.get_positions_raw_2pointperturb", - return_value=[ - np.array([[1.0, 2.2], [1.00101, 2.2], [1.0, 2.201022]]), - np.array([0.00101, 0.201022]), - ], - ) - - m3 = mocker.patch( - "queens.models.simulation_model.SimulationModel.evaluate", - return_value=None, - ) - - default_lm_iterator.model.response = {"result": np.array([[3.0, 4.2], [99.9, 99.9]])} - - m5 = mocker.patch( - "queens.iterators.lm_iterator.fd_jacobian", - return_value=np.array([[1.0, 0.0], [0.0, 1.0]]), - ) - - jacobian, _ = default_lm_iterator.jacobian_and_residual(np.array([1.0, 2.2])) - - m3.assert_called_once() - - m5.assert_called_once() - - np.testing.assert_equal(jacobian, np.array([[1.0, 0.0], [0.0, 1.0]])) - - if fix_true_false_param: - with pytest.raises(ValueError): - m5.return_value = np.array([[1.1, 2.2]]) - default_lm_iterator.jacobian_and_residual(np.array([0.1])) - - -def test_pre_run(mocker, fix_true_false_param, default_lm_iterator, output_csv): - """Test pre-run setup in LMIterator.""" - default_lm_iterator.result_description["write_results"] = fix_true_false_param - - mock_pandas_dataframe_to_csv = mocker.patch("pandas.core.generic.NDFrame.to_csv") - default_lm_iterator.pre_run() - if fix_true_false_param: - mock_pandas_dataframe_to_csv.assert_called_once_with( - output_csv, mode="w", sep="\t", index=None - ) - else: - assert not mock_pandas_dataframe_to_csv.called - default_lm_iterator.result_description = None - default_lm_iterator.pre_run() - - -def test_core_run(default_lm_iterator, mocker, fix_update_reg, fix_tolerance): - """Test core run routine in LMIterator.""" - m1 = mocker.patch( - "queens.iterators.lm_iterator.LMIterator.jacobian_and_residual", - return_value=(np.array([[1.0, 2.0], [0.0, 1.0]]), np.array([0.1, 0.01])), - ) - - m3 = mocker.patch("queens.iterators.lm_iterator.LMIterator.printstep") - default_lm_iterator.update_reg = fix_update_reg - default_lm_iterator.max_feval = 2 - default_lm_iterator.tolerance = fix_tolerance - - if fix_update_reg not in ["grad", "res"]: - with pytest.raises(ValueError): - default_lm_iterator.core_run() - else: - default_lm_iterator.core_run() - if fix_tolerance == 1.0: - assert m1.call_count == 1 - assert m3.call_count == 1 - else: - assert m1.call_count == 3 - assert m3.call_count == 3 - np.testing.assert_almost_equal( - default_lm_iterator.param_current, np.array([-0.00875, 0.15875]), 8 - ) - np.testing.assert_almost_equal(default_lm_iterator.param_opt, np.array([0.1, 0.2]), 8) - - -def test_post_run_2param( - mocker, fix_true_false_param, default_lm_iterator, fix_plotly_fig, output_csv, output_html -): - """Test post-run operations in LMIterator with 2 parameters.""" - default_lm_iterator.solution = np.array([1.1, 2.2]) - default_lm_iterator.iter_opt = 3 - - pdata = pd.DataFrame({"params": ["[1.0e3 2.0e-2]", "[1.1 2.1]"], "resnorm": [1.2, 2.2]}) - checkdata = pd.DataFrame({"resnorm": [1.2, 2.2], "x1": [1000.0, 1.1], "x2": [0.02, 2.1]}) - - default_lm_iterator.result_description["plot_results"] = fix_true_false_param - m1 = mocker.patch("pandas.read_csv", return_value=pdata) - m2 = mocker.patch("plotly.express.line_3d", return_value=fix_plotly_fig) - m3 = mocker.patch("plotly.basedatatypes.BaseFigure.update_traces", return_value=None) - m4 = mocker.patch("plotly.basedatatypes.BaseFigure.write_html", return_value=None) - - default_lm_iterator.post_run() - - if fix_true_false_param: - m1.assert_called_once_with(output_csv, sep="\t") - callargs = m2.call_args - pd.testing.assert_frame_equal(callargs[0][0], checkdata) - assert callargs[1]["x"] == "x1" - assert callargs[1]["y"] == "x2" - assert callargs[1]["z"] == "resnorm" - assert callargs[1]["hover_data"] == [ - "iter", - "resnorm", - "gradnorm", - "delta_params", - "mu", - "x1", - "x2", - ] - m4.assert_called_once_with(output_html) - m2.assert_called_once() - - else: - default_lm_iterator.result_description = None - default_lm_iterator.post_run() - m1.assert_not_called() - m2.assert_not_called() - m3.assert_not_called() - m4.assert_not_called() - - -def test_post_run_1param(mocker, default_lm_iterator, fix_plotly_fig, output_html): - """Test post-run operations in LMIterator with 1 parameter.""" - default_lm_iterator.solution = np.array([1.1, 2.2]) - default_lm_iterator.iter_opt = 3 - - pdata = pd.DataFrame({"params": ["[1.0e3]", "[1.1]"], "resnorm": [1.2, 2.2]}) - mocker.patch("pandas.read_csv", return_value=pdata) - mocker.patch("plotly.basedatatypes.BaseFigure.update_traces", return_value=None) - m4 = mocker.patch("plotly.basedatatypes.BaseFigure.write_html", return_value=None) - m6 = mocker.patch("plotly.express.line", return_value=fix_plotly_fig) - - checkdata = pd.DataFrame({"resnorm": [1.2, 2.2], "x1": [1000.0, 1.1]}) - - default_lm_iterator.post_run() - callargs = m6.call_args - pd.testing.assert_frame_equal(callargs[0][0], checkdata) - assert callargs[1]["x"] == "x1" - assert callargs[1]["y"] == "resnorm" - assert callargs[1]["hover_data"] == [ - "iter", - "resnorm", - "gradnorm", - "delta_params", - "mu", - "x1", - ] - m4.assert_called_once_with(output_html) - m6.assert_called_once() - - -def test_post_run_3param(mocker, default_lm_iterator, caplog): - """Test post-run operations in LMIterator with 3 parameters.""" - default_lm_iterator.solution = np.array([1.1, 2.2]) - default_lm_iterator.iter_opt = 3 - - mocker.patch("plotly.basedatatypes.BaseFigure.update_traces", return_value=None) - m4 = mocker.patch("plotly.basedatatypes.BaseFigure.write_html", return_value=None) - pdata = pd.DataFrame({"params": ["[1.0e3 2.0e-2 3.]", "[1.1 2.1 3.1]"], "resnorm": [1.2, 2.2]}) - mocker.patch("pandas.read_csv", return_value=pdata) - - parameters = Parameters(x1=FreeVariable(1), x2=FreeVariable(1), x3=FreeVariable(1)) - default_lm_iterator.parameters = parameters - - with caplog.at_level(logging.WARNING): - default_lm_iterator.post_run() - - expected_warning = ( - "write_results for more than 2 parameters not implemented, " - "because we are limited to 3 dimensions. " - "You have: 3. Plotting is skipped." - ) - assert expected_warning in caplog.text - - m4.assert_not_called() - - -def test_post_run_0param(mocker, default_lm_iterator): - """Test post-run operations in LMIterator with 0 parameters.""" - default_lm_iterator.solution = np.array([1.1, 2.2]) - default_lm_iterator.iter_opt = 3 - - pdata = pd.DataFrame({"params": ["", ""], "resnorm": [1.2, 2.2]}) - mocker.patch("pandas.read_csv", return_value=pdata) - with pytest.raises(ValueError): - default_lm_iterator.post_run() - - -def test_get_positions_raw_2pointperturb(default_lm_iterator): - """Test get_positions_raw_2pointperturb method in LMIterator.""" - x = np.array([1.1, 2.5]) - pos, delta_pos = default_lm_iterator.get_positions_raw_2pointperturb(x) - np.testing.assert_almost_equal(pos, np.array([[1.1, 2.5], [1.101011, 2.5], [1.1, 2.501025]]), 8) - np.testing.assert_almost_equal(delta_pos, np.array([[0.001011], [0.001025]]), 8) - - default_lm_iterator.bounds = [[0.0, 0.0], [np.inf, 2.5]] - default_lm_iterator.havebounds = True - posb, delta_posb = default_lm_iterator.get_positions_raw_2pointperturb(x) - np.testing.assert_almost_equal( - posb, np.array([[1.1, 2.5], [1.101011, 2.5], [1.1, 2.498975]]), 8 - ) - np.testing.assert_almost_equal(delta_posb, np.array([[0.001011], [-0.001025]]), 8) - - -def test_printstep(mocker, default_lm_iterator, fix_true_false_param, output_csv): - """Test print step output in LMIterator.""" - default_lm_iterator.result_description["write_results"] = fix_true_false_param - - mock_pandas_dataframe_to_csv = mocker.patch("pandas.core.generic.NDFrame.to_csv") - default_lm_iterator.printstep(5, 1e-3, 1e-4, np.array([10.1, 11.2])) - if fix_true_false_param: - mock_pandas_dataframe_to_csv.assert_called_once_with( - output_csv, - header=None, - float_format="%.8f", - mode="a", - sep="\t", - index=None, - ) - - else: - assert not mock_pandas_dataframe_to_csv.called - default_lm_iterator.result_description = None - default_lm_iterator.printstep(5, 1e-3, 1e-4, np.array([10.1, 11.2])) - - -def test_checkbounds(default_lm_iterator, caplog): - """Test bound checking.""" - default_lm_iterator.bounds = np.array([[0.0, 0.0], [5.0, 2.0]]) - with caplog.at_level(logging.WARNING): - stepoutside = default_lm_iterator.checkbounds(np.array([1.0, 2.1]), 3) - - assert stepoutside - assert default_lm_iterator.reg_param == 2.0 - - expected_warning = ( - f"WARNING: STEP #{3} IS OUT OF BOUNDS; double reg_param and compute new iteration.\n " - f"declined step was: {np.array([1.1, 2.3])}" - ) - assert expected_warning in caplog.text