diff --git a/README.md b/README.md index 6c57878..c993b0b 100644 --- a/README.md +++ b/README.md @@ -35,6 +35,6 @@ Documentation can be found at: We absolutely welcome contributions. **EasyScience** is maintained by the ESS and on a volunteer basis and thus we need to foster a community that can support user questions and develop new features to make this software a useful tool for all users while encouraging every member of the community to share their ideas. ## License -While **EasyScience** is under the BSD-3 license, DFO_LS is subject to the GPL license. +While **EasyScience** is under the BSD-3 license, DFO-LS is subject to the GPL license. diff --git a/examples_old/example1.py b/examples_old/example1.py index 867c288..18ca3b5 100644 --- a/examples_old/example1.py +++ b/examples_old/example1.py @@ -3,7 +3,7 @@ import numpy as np -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example1_dream.py b/examples_old/example1_dream.py index c0e0e2c..f485d05 100644 --- a/examples_old/example1_dream.py +++ b/examples_old/example1_dream.py @@ -3,7 +3,7 @@ import numpy as np -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example2.py b/examples_old/example2.py index 07cc009..db1228b 100644 --- a/examples_old/example2.py +++ b/examples_old/example2.py @@ -3,7 +3,7 @@ import numpy as np -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example3.py b/examples_old/example3.py index ef69429..5191612 100644 --- a/examples_old/example3.py +++ b/examples_old/example3.py @@ -5,7 +5,7 @@ import numpy as np -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example4.py b/examples_old/example4.py index f5858f6..7b6e6c9 100644 --- a/examples_old/example4.py +++ b/examples_old/example4.py @@ -10,7 +10,7 @@ import numpy as np from easyscience import borg -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.core import ComponentSerializer from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example5_broken.py b/examples_old/example5_broken.py index 3ac38a9..17abd5d 100644 --- a/examples_old/example5_broken.py +++ b/examples_old/example5_broken.py @@ -9,7 +9,7 @@ import numpy as np from easyscience import borg -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.Base import BaseObj from easyscience.Objects.Base import Parameter from easyscience.Objects.core import ComponentSerializer diff --git a/examples_old/example6_broken.py b/examples_old/example6_broken.py index e3a1183..0c2b9c0 100644 --- a/examples_old/example6_broken.py +++ b/examples_old/example6_broken.py @@ -9,7 +9,7 @@ import numpy as np from easyscience import borg -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.Base import BaseObj from easyscience.Objects.Base import Parameter from easyscience.Objects.core import ComponentSerializer diff --git a/examples_old/example_dataset2.py b/examples_old/example_dataset2.py index 791f76c..e0667a1 100644 --- a/examples_old/example_dataset2.py +++ b/examples_old/example_dataset2.py @@ -5,7 +5,7 @@ import numpy as np from easyscience.Datasets.xarray import xr -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example_dataset2pt2.py b/examples_old/example_dataset2pt2.py index 49c64c5..c33c9a9 100644 --- a/examples_old/example_dataset2pt2.py +++ b/examples_old/example_dataset2pt2.py @@ -5,7 +5,7 @@ import numpy as np from easyscience.Datasets.xarray import xr -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter @@ -39,7 +39,7 @@ def fit_fun(x, *args, **kwargs): f.initialize(b, fit_fun) fig, ax = plt.subplots(2, 3, sharey='row') -for idx, minimizer in enumerate(['lmfit', 'bumps', 'DFO_LS']): +for idx, minimizer in enumerate(['lmfit', 'bumps', 'dfo_ls']): b.m = m_starting_point b.c = c_starting_point diff --git a/examples_old/example_dataset2pt2_broken.py b/examples_old/example_dataset2pt2_broken.py index 7c569d7..b56ea7a 100644 --- a/examples_old/example_dataset2pt2_broken.py +++ b/examples_old/example_dataset2pt2_broken.py @@ -5,7 +5,7 @@ import numpy as np from easyscience.Datasets.xarray import xr -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.Base import BaseObj from easyscience.Objects.Base import Parameter @@ -39,7 +39,7 @@ def fit_fun(x, *args, **kwargs): f.initialize(b, fit_fun) fig, ax = plt.subplots(2, 3, sharey='row') -for idx, minimizer in enumerate(['lmfit', 'bumps', 'DFO_LS']): +for idx, minimizer in enumerate(['lmfit', 'bumps', 'dfo_ls']): b.m = m_starting_point b.c = c_starting_point diff --git a/examples_old/example_dataset3.py b/examples_old/example_dataset3.py index ad0658f..ec1420d 100644 --- a/examples_old/example_dataset3.py +++ b/examples_old/example_dataset3.py @@ -7,7 +7,7 @@ import numpy as np from easyscience.Datasets.xarray import xr -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example_dataset3pt2.py b/examples_old/example_dataset3pt2.py index 060a0a6..a21ff8e 100644 --- a/examples_old/example_dataset3pt2.py +++ b/examples_old/example_dataset3pt2.py @@ -5,7 +5,7 @@ import numpy as np from easyscience.Datasets.xarray import xr -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter @@ -41,7 +41,7 @@ def fit_fun(x, *args, **kwargs): cbar_ax1 = fig.add_axes([0.85, 0.15, 0.05, 0.3]) cbar_ax2 = fig.add_axes([0.85, 0.60, 0.05, 0.3]) -for idx, minimizer in enumerate(['lmfit', 'bumps', 'DFO_LS']): +for idx, minimizer in enumerate(['lmfit', 'bumps', 'dfo_ls']): b.s_off = s_off_start_point b.c_off = c_off_start_point diff --git a/examples_old/example_dataset4.py b/examples_old/example_dataset4.py index 374ae46..a146501 100644 --- a/examples_old/example_dataset4.py +++ b/examples_old/example_dataset4.py @@ -5,7 +5,7 @@ import numpy as np from easyscience.Datasets.xarray import xr -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/examples_old/example_dataset4_2.py b/examples_old/example_dataset4_2.py index 6ebb4cb..63ecfe4 100644 --- a/examples_old/example_dataset4_2.py +++ b/examples_old/example_dataset4_2.py @@ -5,7 +5,7 @@ import numpy as np from easyscience.Datasets.xarray import xr -from easyscience.Fitting.Fitting import Fitter +from easyscience.Fitting import Fitter from easyscience.Objects.ObjectClasses import BaseObj from easyscience.Objects.ObjectClasses import Parameter diff --git a/src/easyscience/Datasets/xarray.py b/src/easyscience/Datasets/xarray.py index 3971444..3663321 100644 --- a/src/easyscience/Datasets/xarray.py +++ b/src/easyscience/Datasets/xarray.py @@ -20,7 +20,7 @@ import xarray as xr from easyscience import ureg -from easyscience.Fitting.fitting_template import FitResults +from easyscience.Fitting import FitResults T_ = TypeVar('T_') diff --git a/src/easyscience/Fitting/__init__.py b/src/easyscience/Fitting/__init__.py index 4a2a542..4c4f3ac 100644 --- a/src/easyscience/Fitting/__init__.py +++ b/src/easyscience/Fitting/__init__.py @@ -1,35 +1,3 @@ -# SPDX-FileCopyrightText: 2023 EasyScience contributors -# SPDX-License-Identifier: BSD-3-Clause -# © 2021-2023 Contributors to the EasyScience project List[str]: :return: List of available fitting engines :rtype: List[str] """ - if Fitting.engines is None: + if minimizers.engines is None: raise ImportError('There are no available fitting engines. Install `lmfit` and/or `bumps`') - return [engine.name for engine in Fitting.engines] + return [engine.name for engine in minimizers.engines] @property def can_fit(self) -> bool: @@ -198,7 +190,7 @@ def fit_function(self, fit_function: Callable): self.__initialize() @property - def fit_object(self) -> B: + def fit_object(self): """ The EasyScience object which will be used as a model :return: EasyScience Model @@ -206,7 +198,7 @@ def fit_object(self) -> B: return self._fit_object @fit_object.setter - def fit_object(self, fit_object: B): + def fit_object(self, fit_object): """ Set the EasyScience object which wil be used as a model :param fit_object: New EasyScience object @@ -248,7 +240,7 @@ def inner_fit_callable( weights: Optional[np.ndarray] = None, vectorized: bool = False, **kwargs, - ) -> FR: + ) -> FitResults: """ This is a wrapped callable which performs the actual fitting. It is split into 3 sections, PRE/ FIT/ POST. @@ -334,7 +326,7 @@ def _precompute_reshaping( return x_for_fit, x_new, y_new, weights, x_shape, kwargs @staticmethod - def _post_compute_reshaping(fit_result: FR, x: np.ndarray, y: np.ndarray, weights: np.ndarray) -> FR: + def _post_compute_reshaping(fit_result: FitResults, x: np.ndarray, y: np.ndarray, weights: np.ndarray) -> FitResults: """ Reshape the output of the fitter into the correct dimensions. :param fit_result: Output from the fitter @@ -347,127 +339,3 @@ def _post_compute_reshaping(fit_result: FR, x: np.ndarray, y: np.ndarray, weight setattr(fit_result, 'y_calc', np.reshape(fit_result.y_calc, y.shape)) setattr(fit_result, 'y_err', np.reshape(fit_result.y_err, y.shape)) return fit_result - - -class MultiFitter(Fitter): - """ - Extension of Fitter to enable multiple dataset/fit function fitting. We can fit these types of data simultaneously: - - Multiple models on multiple datasets. - """ - - def __init__( - self, - fit_objects: Optional[List[B]] = None, - fit_functions: Optional[List[Callable]] = None, - ): - # Create a dummy core object to hold all the fit objects. - self._fit_objects = BaseCollection('multi', *fit_objects) - self._fit_functions = fit_functions - # Initialize with the first of the fit_functions, without this it is - # not possible to change the fitting engine. - super().__init__(self._fit_objects, self._fit_functions[0]) - - def _fit_function_wrapper(self, real_x=None, flatten: bool = True) -> Callable: - """ - Simple fit function which injects the N real X (independent) values into the - optimizer function. This will also flatten the results if needed. - :param real_x: List of independent x parameters to be injected - :param flatten: Should the result be a flat 1D array? - :return: Wrapped optimizer function. - """ - # Extract of a list of callable functions - wrapped_fns = [] - for this_x, this_fun in zip(real_x, self._fit_functions): - self._fit_function = this_fun - wrapped_fns.append(Fitter._fit_function_wrapper(self, this_x, flatten=flatten)) - - def wrapped_fun(x, **kwargs): - # Generate an empty Y based on x - y = np.zeros_like(x) - i = 0 - # Iterate through wrapped functions, passing the WRONG x, the correct - # x was injected in the step above. - for idx, dim in enumerate(self._dependent_dims): - ep = i + np.prod(dim) - y[i:ep] = wrapped_fns[idx](x, **kwargs) - i = ep - return y - - return wrapped_fun - - @staticmethod - def _precompute_reshaping( - x: List[np.ndarray], - y: List[np.ndarray], - weights: Optional[List[np.ndarray]], - vectorized: bool, - kwargs, - ): - """ - Convert an array of X's and Y's to an acceptable shape for fitting. - :param x: List of independent variables. - :param y: List of dependent variables. - :param vectorized: Is the fn input vectorized or point based? - :param kwargs: Additional kwy words. - :return: Variables for optimization - """ - if weights is None: - weights = [None] * len(x) - _, _x_new, _y_new, _weights, _dims, kwargs = Fitter._precompute_reshaping(x[0], y[0], weights[0], vectorized, kwargs) - x_new = [_x_new] - y_new = [_y_new] - w_new = [_weights] - dims = [_dims] - for _x, _y, _w in zip(x[1::], y[1::], weights[1::]): - _, _x_new, _y_new, _weights, _dims, _ = Fitter._precompute_reshaping(_x, _y, _w, vectorized, kwargs) - x_new.append(_x_new) - y_new.append(_y_new) - w_new.append(_weights) - dims.append(_dims) - y_new = np.hstack(y_new) - if w_new[0] is None: - w_new = None - else: - w_new = np.hstack(w_new) - x_fit = np.linspace(0, y_new.size - 1, y_new.size) - return x_fit, x_new, y_new, w_new, dims, kwargs - - def _post_compute_reshaping( - self, - fit_result_obj: FR, - x: List[np.ndarray], - y: List[np.ndarray], - weights: List[np.ndarray], - ) -> List[FR]: - """ - Take a fit results object and split it into n chuncks based on the size of the x, y inputs - :param fit_result_obj: Result from a multifit - :param x: List of X co-ords - :param y: List of Y co-ords - :return: List of fit results - """ - - cls = fit_result_obj.__class__ - sp = 0 - fit_results_list = [] - for idx, this_x in enumerate(x): - # Create a new Results obj - current_results = cls() - ep = sp + int(np.array(self._dependent_dims[idx]).prod()) - - # Fill out the new result obj (see EasyScience.Fitting.Fitting_template.FitResults) - current_results.success = fit_result_obj.success - current_results.fitting_engine = fit_result_obj.fitting_engine - current_results.p = fit_result_obj.p - current_results.p0 = fit_result_obj.p0 - current_results.x = this_x - current_results.y_obs = y[idx] - current_results.y_calc = np.reshape(fit_result_obj.y_calc[sp:ep], current_results.y_obs.shape) - current_results.y_err = np.reshape(fit_result_obj.y_err[sp:ep], current_results.y_obs.shape) - current_results.engine_result = fit_result_obj.engine_result - - # Attach an additional field for the un-modified results - current_results.total_results = fit_result_obj - fit_results_list.append(current_results) - sp = ep - return fit_results_list diff --git a/src/easyscience/Fitting/minimizers/__init__.py b/src/easyscience/Fitting/minimizers/__init__.py new file mode 100644 index 0000000..6afab1c --- /dev/null +++ b/src/easyscience/Fitting/minimizers/__init__.py @@ -0,0 +1,37 @@ +# SPDX-FileCopyrightText: 2023 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause +# © 2021-2023 Contributors to the EasyScience project # SPDX-License-Identifier: BSD-3-Clause # © 2021-2023 Contributors to the EasyScience project str: - return getattr(self._borg.map.get_item_by_key(item_key), 'name', '') - - def get_item_from_key(self, item_key: int) -> object: - return self._borg.map.get_item_by_key(item_key) - - def get_key(self, item: object) -> int: - return self._borg.map.convert_id_to_key(item) - - -class FitError(Exception): - def __init__(self, e: Exception = None): - self.e = e - - def __str__(self) -> str: - s = '' - if self.e is not None: - s = f'{self.e}\n' - return s + 'Something has gone wrong with the fit' diff --git a/src/easyscience/Fitting/bumps.py b/src/easyscience/Fitting/minimizers/minimizer_bumps.py similarity index 83% rename from src/easyscience/Fitting/bumps.py rename to src/easyscience/Fitting/minimizers/minimizer_bumps.py index b7b19a7..05000fd 100644 --- a/src/easyscience/Fitting/bumps.py +++ b/src/easyscience/Fitting/minimizers/minimizer_bumps.py @@ -2,35 +2,35 @@ # SPDX-License-Identifier: BSD-3-Clause # © 2021-2023 Contributors to the EasyScience project Callable: + def make_model(self, pars: Optional[List[BumpsParameter]] = None) -> Callable: """ Generate a bumps model from the supplied `fit_function` and parameters in the base object. Note that this makes a callable as it needs to be initialized with *x*, *y*, *weights* @@ -62,12 +62,10 @@ def make_func(x, y, weights): par = {} if not pars: for name, item in obj._cached_pars.items(): - par["p" + str(name)] = obj.convert_to_par_object(item) + par['p' + str(name)] = obj.convert_to_par_object(item) else: for item in pars: - par[ - "p" + str(NameConverter().get_key(item)) - ] = obj.convert_to_par_object(item) + par['p' + str(NameConverter().get_key(item))] = obj.convert_to_par_object(item) return Curve(fit_func, x, y, dy=weights, **par) return make_func @@ -125,12 +123,10 @@ def fit_function(x: np.ndarray, **kwargs): self._cached_pars_order = tuple(self._cached_pars.keys()) params = [ - inspect.Parameter( - "x", inspect.Parameter.POSITIONAL_OR_KEYWORD, annotation=inspect._empty - ), + inspect.Parameter('x', inspect.Parameter.POSITIONAL_OR_KEYWORD, annotation=inspect._empty), *[ inspect.Parameter( - "p" + str(name), + 'p' + str(name), inspect.Parameter.POSITIONAL_OR_KEYWORD, annotation=inspect._empty, default=self._cached_pars[name].raw_value, @@ -148,8 +144,8 @@ def fit( x: np.ndarray, y: np.ndarray, weights: Optional[np.ndarray] = None, - model: Optional = None, - parameters: Optional = None, + model=None, + parameters=None, method: Optional[str] = None, minimizer_kwargs: Optional[dict] = None, engine_kwargs: Optional[dict] = None, @@ -167,7 +163,7 @@ def fit( :param model: Optional Model which is being fitted to :type model: lmModel :param parameters: Optional parameters for the fit - :type parameters: List[bumpsParameter] + :type parameters: List[BumpsParameter] :param kwargs: Additional arguments for the fitting function. :param method: Method for minimization :type method: str @@ -177,7 +173,7 @@ def fit( default_method = {} if method is not None and method in self.available_methods(): - default_method["method"] = method + default_method['method'] = method if weights is None: weights = np.sqrt(np.abs(y)) @@ -195,10 +191,7 @@ def fit( model = self.make_model(pars=parameters) model = model(x, y, weights) self._cached_model = model - self.p_0 = { - f"p{key}": self._cached_pars[key].raw_value - for key in self._cached_pars.keys() - } + self._p_0 = {f'p{key}': self._cached_pars[key].raw_value for key in self._cached_pars.keys()} problem = FitProblem(model) # Why do we do this? Because a fitting template has to have borg instantiated outside pre-runtime from easyscience import borg @@ -207,9 +200,7 @@ def fit( borg.stack.enabled = False try: - model_results = bumps_fit( - problem, **default_method, **minimizer_kwargs, **kwargs - ) + model_results = bumps_fit(problem, **default_method, **minimizer_kwargs, **kwargs) self._set_parameter_fit_result(model_results, stack_status) results = self._gen_fit_results(model_results) except Exception as e: @@ -218,16 +209,14 @@ def fit( raise FitError(e) return results - def convert_to_pars_obj( - self, par_list: Optional[List] = None - ) -> List[bumpsParameter]: + def convert_to_pars_obj(self, par_list: Optional[List] = None) -> List[BumpsParameter]: """ Create a container with the `Parameters` converted from the base object. :param par_list: If only a single/selection of parameter is required. Specify as a list :type par_list: List[str] :return: bumps Parameters list - :rtype: List[bumpsParameter] + :rtype: List[BumpsParameter] """ if par_list is None: # Assume that we have a BaseObj for which we can obtain a list @@ -237,15 +226,15 @@ def convert_to_pars_obj( # For some reason I have to double staticmethod :-/ @staticmethod - def convert_to_par_object(obj) -> bumpsParameter: + def convert_to_par_object(obj) -> BumpsParameter: """ Convert an `EasyScience.Objects.Base.Parameter` object to a bumps Parameter object :return: bumps Parameter compatible object. - :rtype: bumpsParameter + :rtype: BumpsParameter """ - return bumpsParameter( - name="p" + str(NameConverter().get_key(obj)), + return BumpsParameter( + name='p' + str(NameConverter().get_key(obj)), value=obj.raw_value, bounds=[obj.min, obj.max], fixed=obj.fixed, @@ -268,7 +257,7 @@ def _set_parameter_fit_result(self, fit_result, stack_status: bool): pars[name].value = self._cached_pars_vals[name][0] pars[name].error = self._cached_pars_vals[name][1] borg.stack.enabled = True - borg.stack.beginMacro("Fitting routine") + borg.stack.beginMacro('Fitting routine') for index, name in enumerate(self._cached_model._pnames): dict_name = int(name[1:]) @@ -296,7 +285,7 @@ def _gen_fit_results(self, fit_results, **kwargs) -> FitResults: for index, name in enumerate(self._cached_model._pnames): dict_name = int(name[1:]) item[name] = pars[dict_name].raw_value - results.p0 = self.p_0 + results.p0 = self._p_0 results.p = item results.x = self._cached_model.x results.y_obs = self._cached_model.y diff --git a/src/easyscience/Fitting/DFO_LS.py b/src/easyscience/Fitting/minimizers/minimizer_dfo.py similarity index 86% rename from src/easyscience/Fitting/DFO_LS.py rename to src/easyscience/Fitting/minimizers/minimizer_dfo.py index 53ab0ff..48afbff 100644 --- a/src/easyscience/Fitting/DFO_LS.py +++ b/src/easyscience/Fitting/minimizers/minimizer_dfo.py @@ -2,31 +2,30 @@ # SPDX-License-Identifier: BSD-3-Clause # © 2021-2023 Contributors to the EasyScience project Callable: """ @@ -57,18 +56,18 @@ def make_func(x, y, weights): par = {} if not pars: for name, item in obj._cached_pars.items(): - par["p" + str(name)] = item.raw_value + par['p' + str(name)] = item.raw_value else: for item in pars: - par["p" + str(NameConverter().get_key(item))] = item.raw_value + par['p' + str(NameConverter().get_key(item))] = item.raw_value def residuals(x0) -> np.ndarray: for idx, par_name in enumerate(par.keys()): par[par_name] = x0[idx] return (y - fit_func(x, **par)) / weights - setattr(residuals, "x", x) - setattr(residuals, "y", y) + setattr(residuals, 'x', x) + setattr(residuals, 'y', y) return residuals return make_func @@ -157,7 +156,7 @@ def fit( default_method = {} if method is not None and method in self.available_methods(): - default_method["method"] = method + default_method['method'] = method if weights is None: weights = np.sqrt(np.abs(y)) @@ -166,10 +165,7 @@ def fit( model = self.make_model(pars=parameters) model = model(x, y, weights) self._cached_model = model - self.p_0 = { - f"p{key}": self._cached_pars[key].raw_value - for key in self._cached_pars.keys() - } + self._p_0 = {f'p{key}': self._cached_pars[key].raw_value for key in self._cached_pars.keys()} # Why do we do this? Because a fitting template has to have borg instantiated outside pre-runtime from easyscience import borg @@ -202,9 +198,7 @@ def convert_to_par_object(obj) -> None: """ pass - def _set_parameter_fit_result( - self, fit_result, stack_status, ci: float = 0.95 - ) -> None: + def _set_parameter_fit_result(self, fit_result, stack_status, ci: float = 0.95) -> None: """ Update parameters to their final values and assign a std error to them. @@ -221,11 +215,9 @@ def _set_parameter_fit_result( pars[name].value = self._cached_pars_vals[name][0] pars[name].error = self._cached_pars_vals[name][1] borg.stack.enabled = True - borg.stack.beginMacro("Fitting routine") + borg.stack.beginMacro('Fitting routine') - error_matrix = self._error_from_jacobian( - fit_result.jacobian, fit_result.resid, ci - ) + error_matrix = self._error_from_jacobian(fit_result.jacobian, fit_result.resid, ci) for idx, par in enumerate(pars.values()): par.value = fit_result.x[idx] par.error = error_matrix[idx, idx] @@ -250,8 +242,8 @@ def _gen_fit_results(self, fit_results, weights, **kwargs) -> FitResults: pars = self._cached_pars item = {} for p_name, par in pars.items(): - item[f"p{p_name}"] = par.raw_value - results.p0 = self.p_0 + item[f'p{p_name}'] = par.raw_value + results.p0 = self._p_0 results.p = item results.x = self._cached_model.x results.y_obs = self._cached_model.y @@ -267,7 +259,7 @@ def _gen_fit_results(self, fit_results, weights, **kwargs) -> FitResults: return results def available_methods(self) -> List[str]: - return ["leastsq"] + return ['leastsq'] def dfols_fit(self, model: Callable, **kwargs): """ diff --git a/src/easyscience/Fitting/lmfit.py b/src/easyscience/Fitting/minimizers/minimizer_lmfit.py similarity index 73% rename from src/easyscience/Fitting/lmfit.py rename to src/easyscience/Fitting/minimizers/minimizer_lmfit.py index c77e603..0adad4a 100644 --- a/src/easyscience/Fitting/lmfit.py +++ b/src/easyscience/Fitting/minimizers/minimizer_lmfit.py @@ -2,76 +2,59 @@ # SPDX-License-Identifier: BSD-3-Clause # © 2021-2023 Contributors to the EasyScience project lmModel: + def make_model(self, pars: Optional[LMParameters] = None) -> LMModel: """ Generate a lmfit model from the supplied `fit_function` and parameters in the base object. :return: Callable lmfit model - :rtype: lmModel + :rtype: LMModel """ # Generate the fitting function fit_func = self._generate_fit_function() if pars is None: pars = self._cached_pars # Create the model - model = lmModel( + model = LMModel( fit_func, - independent_vars=["x"], - param_names=["p" + str(key) for key in pars.keys()], + independent_vars=['x'], + param_names=['p' + str(key) for key in pars.keys()], ) # Assign values from the `Parameter` to the model for name, item in pars.items(): - if isinstance(item, lmParameter): + if isinstance(item, LMParameter): value = item.value else: value = item.raw_value - model.set_param_hint( - "p" + str(name), value=value, min=item.min, max=item.max - ) + model.set_param_hint('p' + str(name), value=value, min=item.min, max=item.max) # Cache the model for later reference self._cached_model = model @@ -127,12 +110,10 @@ def fit_function(x: np.ndarray, **kwargs): # f = (x, a=1, b=2)... # Where we need to be generic. Note that this won't hold for much outside of this scope. params = [ - inspect.Parameter( - "x", inspect.Parameter.POSITIONAL_OR_KEYWORD, annotation=inspect._empty - ), + inspect.Parameter('x', inspect.Parameter.POSITIONAL_OR_KEYWORD, annotation=inspect._empty), *[ inspect.Parameter( - "p" + str(name), + 'p' + str(name), inspect.Parameter.POSITIONAL_OR_KEYWORD, annotation=inspect._empty, default=parameter.raw_value, @@ -150,8 +131,8 @@ def fit( x: np.ndarray, y: np.ndarray, weights: Optional[np.ndarray] = None, - model: Optional[lmModel] = None, - parameters: Optional[lmParameters] = None, + model: Optional[LMModel] = None, + parameters: Optional[LMParameters] = None, method: Optional[str] = None, minimizer_kwargs: Optional[dict] = None, engine_kwargs: Optional[dict] = None, @@ -169,9 +150,9 @@ def fit( :param weights: Weights for supplied measured points :type weights: np.ndarray :param model: Optional Model which is being fitted to - :type model: lmModel + :type model: LMModel :param parameters: Optional parameters for the fit - :type parameters: lmParameters + :type parameters: LMParameters :param minimizer_kwargs: Arguments to be passed directly to the minimizer :type minimizer_kwargs: dict :param kwargs: Additional arguments for the fitting function. @@ -180,7 +161,7 @@ def fit( """ default_method = {} if method is not None and method in self.available_methods(): - default_method["method"] = method + default_method['method'] = method if weights is None: weights = 1 / np.sqrt(np.abs(y)) @@ -191,7 +172,7 @@ def fit( if minimizer_kwargs is None: minimizer_kwargs = {} else: - minimizer_kwargs = {"fit_kws": minimizer_kwargs} + minimizer_kwargs = {'fit_kws': minimizer_kwargs} minimizer_kwargs.update(engine_kwargs) # Why do we do this? Because a fitting template has to have borg instantiated outside pre-runtime @@ -204,9 +185,7 @@ def fit( if model is None: model = self.make_model() - model_results = model.fit( - y, x=x, weights=weights, **default_method, **minimizer_kwargs, **kwargs - ) + model_results = model.fit(y, x=x, weights=weights, **default_method, **minimizer_kwargs, **kwargs) self._set_parameter_fit_result(model_results, stack_status) results = self._gen_fit_results(model_results) except Exception as e: @@ -215,33 +194,31 @@ def fit( raise FitError(e) return results - def convert_to_pars_obj(self, par_list: Optional[List] = None) -> lmParameters: + def convert_to_pars_obj(self, par_list: Optional[List] = None) -> LMParameters: """ Create an lmfit compatible container with the `Parameters` converted from the base object. :param par_list: If only a single/selection of parameter is required. Specify as a list :type par_list: List[str] :return: lmfit Parameters compatible object - :rtype: lmParameters + :rtype: LMParameters """ if par_list is None: # Assume that we have a BaseObj for which we can obtain a list par_list = self._object.get_fit_parameters() - pars_obj = lmParameters().add_many( - [self.__class__.convert_to_par_object(obj) for obj in par_list] - ) + pars_obj = LMParameters().add_many([self.__class__.convert_to_par_object(obj) for obj in par_list]) return pars_obj @staticmethod - def convert_to_par_object(obj) -> lmParameter: + def convert_to_par_object(obj) -> LMParameter: """ Convert an `EasyScience.Objects.Base.Parameter` object to a lmfit Parameter object. :return: lmfit Parameter compatible object. - :rtype: lmParameter + :rtype: LMParameter """ - return lmParameter( - "p" + str(NameConverter().get_key(obj)), + return LMParameter( + 'p' + str(NameConverter().get_key(obj)), value=obj.raw_value, vary=not obj.fixed, min=obj.min, @@ -266,11 +243,11 @@ def _set_parameter_fit_result(self, fit_result: ModelResult, stack_status: bool) pars[name].value = self._cached_pars_vals[name][0] pars[name].error = self._cached_pars_vals[name][1] borg.stack.enabled = True - borg.stack.beginMacro("Fitting routine") + borg.stack.beginMacro('Fitting routine') for name in pars.keys(): - pars[name].value = fit_result.params["p" + str(name)].value + pars[name].value = fit_result.params['p' + str(name)].value if fit_result.errorbars: - pars[name].error = fit_result.params["p" + str(name)].stderr + pars[name].error = fit_result.params['p' + str(name)].stderr else: pars[name].error = 0.0 if stack_status: @@ -294,7 +271,7 @@ def _gen_fit_results(self, fit_results: ModelResult, **kwargs) -> FitResults: results.success = fit_results.success results.y_obs = fit_results.data # results.residual = fit_results.residual - results.x = fit_results.userkws["x"] + results.x = fit_results.userkws['x'] results.p = fit_results.values results.p0 = fit_results.init_values # results.goodness_of_fit = fit_results.chisqr @@ -309,16 +286,16 @@ def _gen_fit_results(self, fit_results: ModelResult, **kwargs) -> FitResults: def available_methods(self) -> List[str]: return [ - "least_squares", - "leastsq", - "differential_evolution", - "basinhopping", - "ampgo", - "nelder", - "lbfgsb", - "powell", - "cg", - "newton", - "cobyla", - "bfgs", + 'least_squares', + 'leastsq', + 'differential_evolution', + 'basinhopping', + 'ampgo', + 'nelder', + 'lbfgsb', + 'powell', + 'cg', + 'newton', + 'cobyla', + 'bfgs', ] diff --git a/src/easyscience/Fitting/minimizers/utils.py b/src/easyscience/Fitting/minimizers/utils.py new file mode 100644 index 0000000..ff56e87 --- /dev/null +++ b/src/easyscience/Fitting/minimizers/utils.py @@ -0,0 +1,79 @@ +import numpy as np + + +class FitResults: + """ + At the moment this is just a dummy way of unifying the returned fit parameters. + """ + + __slots__ = [ + 'success', + 'fitting_engine', + 'fit_args', + 'p', + 'p0', + 'x', + 'x_matrices', + 'y_obs', + 'y_calc', + 'y_err', + 'engine_result', + 'total_results', + ] + + def __init__(self): + self.success = False + self.fitting_engine = None + self.fit_args = {} + self.p = {} + self.p0 = {} + self.x = np.ndarray([]) + self.x_matrices = np.ndarray([]) + self.y_obs = np.ndarray([]) + self.y_calc = np.ndarray([]) + self.y_err = np.ndarray([]) + self.engine_result = None + self.total_results = None + + @property + def n_pars(self): + return len(self.p) + + @property + def residual(self): + return self.y_obs - self.y_calc + + @property + def chi2(self): + return ((self.residual / self.y_err) ** 2).sum() + + @property + def reduced_chi(self): + return self.chi2 / (len(self.x) - self.n_pars) + + +class NameConverter: + def __init__(self): + from easyscience import borg + + self._borg = borg + + def get_name_from_key(self, item_key: int) -> str: + return getattr(self._borg.map.get_item_by_key(item_key), 'name', '') + + def get_item_from_key(self, item_key: int) -> object: + return self._borg.map.get_item_by_key(item_key) + + def get_key(self, item: object) -> int: + return self._borg.map.convert_id_to_key(item) + + +class FitError(Exception): + def __init__(self, e: Exception = None): + self.e = e + + def __str__(self) -> str: + s = '' + if self.e is not None: + s = f'{self.e}\n' + return s + 'Something has gone wrong with the fit' diff --git a/src/easyscience/Fitting/multi_fitter.py b/src/easyscience/Fitting/multi_fitter.py new file mode 100644 index 0000000..f7f89fd --- /dev/null +++ b/src/easyscience/Fitting/multi_fitter.py @@ -0,0 +1,140 @@ +__author__ = 'github.com/wardsimon' +__version__ = '0.0.1' + +# SPDX-FileCopyrightText: 2023 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause +# © 2021-2023 Contributors to the EasyScience project Callable: + """ + Simple fit function which injects the N real X (independent) values into the + optimizer function. This will also flatten the results if needed. + :param real_x: List of independent x parameters to be injected + :param flatten: Should the result be a flat 1D array? + :return: Wrapped optimizer function. + """ + # Extract of a list of callable functions + wrapped_fns = [] + for this_x, this_fun in zip(real_x, self._fit_functions): + self._fit_function = this_fun + wrapped_fns.append(Fitter._fit_function_wrapper(self, this_x, flatten=flatten)) + + def wrapped_fun(x, **kwargs): + # Generate an empty Y based on x + y = np.zeros_like(x) + i = 0 + # Iterate through wrapped functions, passing the WRONG x, the correct + # x was injected in the step above. + for idx, dim in enumerate(self._dependent_dims): + ep = i + np.prod(dim) + y[i:ep] = wrapped_fns[idx](x, **kwargs) + i = ep + return y + + return wrapped_fun + + @staticmethod + def _precompute_reshaping( + x: List[np.ndarray], + y: List[np.ndarray], + weights: Optional[List[np.ndarray]], + vectorized: bool, + kwargs, + ): + """ + Convert an array of X's and Y's to an acceptable shape for fitting. + :param x: List of independent variables. + :param y: List of dependent variables. + :param vectorized: Is the fn input vectorized or point based? + :param kwargs: Additional kwy words. + :return: Variables for optimization + """ + if weights is None: + weights = [None] * len(x) + _, _x_new, _y_new, _weights, _dims, kwargs = Fitter._precompute_reshaping(x[0], y[0], weights[0], vectorized, kwargs) + x_new = [_x_new] + y_new = [_y_new] + w_new = [_weights] + dims = [_dims] + for _x, _y, _w in zip(x[1::], y[1::], weights[1::]): + _, _x_new, _y_new, _weights, _dims, _ = Fitter._precompute_reshaping(_x, _y, _w, vectorized, kwargs) + x_new.append(_x_new) + y_new.append(_y_new) + w_new.append(_weights) + dims.append(_dims) + y_new = np.hstack(y_new) + if w_new[0] is None: + w_new = None + else: + w_new = np.hstack(w_new) + x_fit = np.linspace(0, y_new.size - 1, y_new.size) + return x_fit, x_new, y_new, w_new, dims, kwargs + + def _post_compute_reshaping( + self, + fit_result_obj: FitResults, + x: List[np.ndarray], + y: List[np.ndarray], + weights: List[np.ndarray], + ) -> List[FitResults]: + """ + Take a fit results object and split it into n chuncks based on the size of the x, y inputs + :param fit_result_obj: Result from a multifit + :param x: List of X co-ords + :param y: List of Y co-ords + :return: List of fit results + """ + + cls = fit_result_obj.__class__ + sp = 0 + fit_results_list = [] + for idx, this_x in enumerate(x): + # Create a new Results obj + current_results = cls() + ep = sp + int(np.array(self._dependent_dims[idx]).prod()) + + # Fill out the new result obj (see EasyScience.Fitting.Fitting_template.FitResults) + current_results.success = fit_result_obj.success + current_results.fitting_engine = fit_result_obj.fitting_engine + current_results.p = fit_result_obj.p + current_results.p0 = fit_result_obj.p0 + current_results.x = this_x + current_results.y_obs = y[idx] + current_results.y_calc = np.reshape(fit_result_obj.y_calc[sp:ep], current_results.y_obs.shape) + current_results.y_err = np.reshape(fit_result_obj.y_err[sp:ep], current_results.y_obs.shape) + current_results.engine_result = fit_result_obj.engine_result + + # Attach an additional field for the un-modified results + current_results.total_results = fit_result_obj + fit_results_list.append(current_results) + sp = ep + return fit_results_list diff --git a/src/easyscience/Objects/Inferface.py b/src/easyscience/Objects/Inferface.py index bc9c3eb..2ff89dd 100644 --- a/src/easyscience/Objects/Inferface.py +++ b/src/easyscience/Objects/Inferface.py @@ -4,8 +4,8 @@ # SPDX-License-Identifier: BSD-3-Clause # © 2021-2023 Contributors to the EasyScience project 0: # Fallback name interface_name = self.return_name(self._interfaces[0]) else: raise NotImplementedError else: - interface_name = kwargs.pop("interface_name") + interface_name = kwargs.pop('interface_name') interfaces = self.available_interfaces if interface_name in interfaces: self._current_interface = self._interfaces[interfaces.index(interface_name)] @@ -74,20 +74,20 @@ def switch(self, new_interface: str, fitter: Optional[Type[Fitter]] = None): self._current_interface = self._interfaces[interfaces.index(new_interface)] self.__interface_obj = self._current_interface() else: - raise AttributeError("The user supplied interface is not valid.") + raise AttributeError('The user supplied interface is not valid.') if fitter is not None: - if hasattr(fitter, "_fit_object"): - obj = getattr(fitter, "_fit_object") + if hasattr(fitter, '_fit_object'): + obj = getattr(fitter, '_fit_object') try: - if hasattr(obj, "update_bindings"): + if hasattr(obj, 'update_bindings'): obj.update_bindings() except Exception as e: - print(f"Unable to auto generate bindings.\n{e}") - elif hasattr(fitter, "generate_bindings"): + print(f'Unable to auto generate bindings.\n{e}') + elif hasattr(fitter, 'generate_bindings'): try: fitter.generate_bindings() except Exception as e: - print(f"Unable to auto generate bindings.\n{e}") + print(f'Unable to auto generate bindings.\n{e}') @property def available_interfaces(self) -> List[str]: @@ -189,8 +189,8 @@ def return_name(this_interface) -> str: Return an interfaces name """ interface_name = this_interface.__name__ - if hasattr(this_interface, "name"): - interface_name = getattr(this_interface, "name") + if hasattr(this_interface, 'name'): + interface_name = getattr(this_interface, 'name') return interface_name @@ -225,4 +225,4 @@ def set_value(value): return set_value -iF = TypeVar("iF", bound=InterfaceFactoryTemplate) +iF = TypeVar('iF', bound=InterfaceFactoryTemplate) diff --git a/tests/integration_tests/Fitting/test_fitter.py b/tests/integration_tests/Fitting/test_fitter.py new file mode 100644 index 0000000..0180ef3 --- /dev/null +++ b/tests/integration_tests/Fitting/test_fitter.py @@ -0,0 +1,278 @@ +# SPDX-FileCopyrightText: 2023 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause +# © 2021-2023 Contributors to the EasyScience project 0 % This does not work as some methods don't calculate error + assert item1.error == pytest.approx(0, abs=1e-1) + assert item1.raw_value == pytest.approx(item2.raw_value, abs=5e-3) + y_calc_ref = ref_sin(x) + assert result.y_calc == pytest.approx(y_calc_ref, abs=1e-2) + assert result.residual == pytest.approx(sp_sin(x) - y_calc_ref, abs=1e-2) + + +@pytest.mark.parametrize("with_errors", [False, True]) +@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "dfo_ls"]) +def test_basic_fit(fit_engine, with_errors): + ref_sin = AbsSin(0.2, np.pi) + sp_sin = AbsSin(0.354, 3.05) + + x = np.linspace(0, 5, 200) + y = ref_sin(x) + + sp_sin.offset.fixed = False + sp_sin.phase.fixed = False + + f = Fitter(sp_sin, sp_sin) + if fit_engine is not None: + try: + f.switch_engine(fit_engine) + except AttributeError: + pytest.skip(msg=f"{fit_engine} is not installed") + args = [x, y] + kwargs = {} + if with_errors: + kwargs["weights"] = 1 / np.sqrt(y) + result = f.fit(*args, **kwargs) + + if fit_engine is not None: + assert result.fitting_engine.name == fit_engine + assert sp_sin.phase.raw_value == pytest.approx(ref_sin.phase.raw_value, rel=1e-3) + assert sp_sin.offset.raw_value == pytest.approx(ref_sin.offset.raw_value, rel=1e-3) + + +@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "dfo_ls"]) +def test_fit_result(fit_engine): + ref_sin = AbsSin(0.2, np.pi) + sp_sin = AbsSin(0.354, 3.05) + + x = np.linspace(0, 5, 200) + y = ref_sin(x) + + sp_sin.offset.fixed = False + sp_sin.phase.fixed = False + + sp_ref1 = { + f"p{sp_sin._borg.map.convert_id_to_key(item1)}": item1.raw_value + for item1, item2 in zip(sp_sin._kwargs.values(), ref_sin._kwargs.values()) + } + sp_ref2 = { + f"p{sp_sin._borg.map.convert_id_to_key(item1)}": item2.raw_value + for item1, item2 in zip(sp_sin._kwargs.values(), ref_sin._kwargs.values()) + } + + f = Fitter(sp_sin, sp_sin) + + if fit_engine is not None: + try: + f.switch_engine(fit_engine) + except AttributeError: + pytest.skip(msg=f"{fit_engine} is not installed") + + result = f.fit(x, y) + check_fit_results(result, sp_sin, ref_sin, x, sp_ref1=sp_ref1, sp_ref2=sp_ref2) + + +@pytest.mark.parametrize("fit_method", ["leastsq", "powell", "cobyla"]) +def test_lmfit_methods(fit_method): + ref_sin = AbsSin(0.2, np.pi) + sp_sin = AbsSin(0.354, 3.05) + + x = np.linspace(0, 5, 200) + y = ref_sin(x) + + sp_sin.offset.fixed = False + sp_sin.phase.fixed = False + + f = Fitter(sp_sin, sp_sin) + assert fit_method in f.available_methods() + result = f.fit(x, y, method=fit_method) + check_fit_results(result, sp_sin, ref_sin, x) + + +@pytest.mark.xfail(reason="known bumps issue") +@pytest.mark.parametrize("fit_method", ["newton", "lm"]) +def test_bumps_methods(fit_method): + ref_sin = AbsSin(0.2, np.pi) + sp_sin = AbsSin(0.354, 3.05) + + x = np.linspace(0, 5, 200) + y = ref_sin(x) + + sp_sin.offset.fixed = False + sp_sin.phase.fixed = False + + f = Fitter(sp_sin, sp_sin) + f.switch_engine("bumps") + assert fit_method in f.available_methods() + result = f.fit(x, y, method=fit_method) + check_fit_results(result, sp_sin, ref_sin, x) + + +@pytest.mark.parametrize("fit_engine", ["lmfit", "bumps", "dfo_ls"]) +def test_fit_constraints(fit_engine): + ref_sin = AbsSin(np.pi * 0.45, 0.45 * np.pi * 0.5) + sp_sin = AbsSin(1, 0.5) + + x = np.linspace(0, 5, 200) + y = ref_sin(x) + + sp_sin.phase.fixed = False + + f = Fitter(sp_sin, sp_sin) + + assert len(f.fit_constraints()) == 0 + c = ObjConstraint(sp_sin.offset, "2*", sp_sin.phase) + f.add_fit_constraint(c) + + if fit_engine is not None: + try: + f.switch_engine(fit_engine) + except AttributeError: + pytest.skip(msg=f"{fit_engine} is not installed") + + result = f.fit(x, y) + check_fit_results(result, sp_sin, ref_sin, x) + assert len(f.fit_constraints()) == 1 + f.remove_fit_constraint(0) + assert len(f.fit_constraints()) == 0 + + +@pytest.mark.parametrize("with_errors", [False, True]) +@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "dfo_ls"]) +def test_2D_vectorized(fit_engine, with_errors): + x = np.linspace(0, 5, 200) + mm = AbsSin2D(0.3, 1.6) + m2 = AbsSin2D( + 0.1, 1.8 + ) # The fit is quite sensitive to the initial values :-( + X, Y = np.meshgrid(x, x) + XY = np.stack((X, Y), axis=2) + ff = Fitter(m2, m2) + if fit_engine is not None: + try: + ff.switch_engine(fit_engine) + except AttributeError: + pytest.skip(msg=f"{fit_engine} is not installed") + try: + args = [XY, mm(XY)] + kwargs = {"vectorized": True} + if with_errors: + kwargs["weights"] = 1 / np.sqrt(args[1]) + result = ff.fit(*args, **kwargs) + except FitError as e: + if "Unable to allocate" in str(e): + pytest.skip(msg="MemoryError - Matrix too large") + else: + raise e + assert result.n_pars == len(m2.get_fit_parameters()) + assert result.reduced_chi == pytest.approx(0, abs=1.5e-3) + assert result.success + assert np.all(result.x == XY) + y_calc_ref = m2(XY) + assert result.y_calc == pytest.approx(y_calc_ref, abs=1e-2) + assert result.residual == pytest.approx(mm(XY) - y_calc_ref, abs=1e-2) + + +@pytest.mark.parametrize("with_errors", [False, True]) +@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "dfo_ls"]) +def test_2D_non_vectorized(fit_engine, with_errors): + x = np.linspace(0, 5, 200) + mm = AbsSin2DL(0.3, 1.6) + m2 = AbsSin2DL( + 0.1, 1.8 + ) # The fit is quite sensitive to the initial values :-( + X, Y = np.meshgrid(x, x) + XY = np.stack((X, Y), axis=2) + ff = Fitter(m2, m2) + if fit_engine is not None: + try: + ff.switch_engine(fit_engine) + except AttributeError: + pytest.skip(msg=f"{fit_engine} is not installed") + try: + args = [XY, mm(XY.reshape(-1, 2))] + kwargs = {"vectorized": False} + if with_errors: + kwargs["weights"] = 1 / np.sqrt(args[1]) + result = ff.fit(*args, **kwargs) + except FitError as e: + if "Unable to allocate" in str(e): + pytest.skip(msg="MemoryError - Matrix too large") + else: + raise e + assert result.n_pars == len(m2.get_fit_parameters()) + assert result.reduced_chi == pytest.approx(0, abs=1.5e-3) + assert result.success + assert np.all(result.x == XY) + y_calc_ref = m2(XY.reshape(-1, 2)) + assert result.y_calc == pytest.approx(y_calc_ref, abs=1e-2) + assert result.residual == pytest.approx( + mm(XY.reshape(-1, 2)) - y_calc_ref, abs=1e-2 + ) diff --git a/tests/integration_tests/Fitting/test_multi_fitter.py b/tests/integration_tests/Fitting/test_multi_fitter.py new file mode 100644 index 0000000..5ec8de9 --- /dev/null +++ b/tests/integration_tests/Fitting/test_multi_fitter.py @@ -0,0 +1,279 @@ +# SPDX-FileCopyrightText: 2023 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause +# © 2021-2023 Contributors to the EasyScience project np.ndarray: y = l1(x) + 0.125 * (dy - 0.5) - from easyscience.Fitting.Fitting import Fitter + from easyscience.Fitting import Fitter f = Fitter(l2, l2) try: diff --git a/tests/unit_tests/Fitting/test_fitting.py b/tests/unit_tests/Fitting/test_fitting.py deleted file mode 100644 index fba4fd0..0000000 --- a/tests/unit_tests/Fitting/test_fitting.py +++ /dev/null @@ -1,573 +0,0 @@ -# SPDX-FileCopyrightText: 2023 EasyScience contributors -# SPDX-License-Identifier: BSD-3-Clause -# © 2021-2023 Contributors to the EasyScience project 0 % This does not work as some methods don't calculate error - assert item1.error == pytest.approx(0, abs=1e-1) - assert item1.raw_value == pytest.approx(item2.raw_value, abs=5e-3) - y_calc_ref = ref_sin(x) - assert result.y_calc == pytest.approx(y_calc_ref, abs=1e-2) - assert result.residual == pytest.approx(sp_sin(x) - y_calc_ref, abs=1e-2) - - -@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "DFO_LS"]) -def test_fit_result(genObjs, fit_engine): - ref_sin = genObjs[0] - sp_sin = genObjs[1] - - x = np.linspace(0, 5, 200) - y = ref_sin(x) - - sp_sin.offset.fixed = False - sp_sin.phase.fixed = False - - sp_ref1 = { - f"p{sp_sin._borg.map.convert_id_to_key(item1)}": item1.raw_value - for item1, item2 in zip(sp_sin._kwargs.values(), ref_sin._kwargs.values()) - } - sp_ref2 = { - f"p{sp_sin._borg.map.convert_id_to_key(item1)}": item2.raw_value - for item1, item2 in zip(sp_sin._kwargs.values(), ref_sin._kwargs.values()) - } - - f = Fitter(sp_sin, sp_sin) - - if fit_engine is not None: - try: - f.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - - result = f.fit(x, y) - check_fit_results(result, sp_sin, ref_sin, x, sp_ref1=sp_ref1, sp_ref2=sp_ref2) - - -@pytest.mark.parametrize("fit_method", ["leastsq", "powell", "cobyla"]) -def test_lmfit_methods(genObjs, fit_method): - ref_sin = genObjs[0] - sp_sin = genObjs[1] - - x = np.linspace(0, 5, 200) - y = ref_sin(x) - - sp_sin.offset.fixed = False - sp_sin.phase.fixed = False - - f = Fitter(sp_sin, sp_sin) - assert fit_method in f.available_methods() - result = f.fit(x, y, method=fit_method) - check_fit_results(result, sp_sin, ref_sin, x) - - -@pytest.mark.xfail(reason="known bumps issue") -@pytest.mark.parametrize("fit_method", ["newton", "lm"]) -def test_bumps_methods(genObjs, fit_method): - ref_sin = genObjs[0] - sp_sin = genObjs[1] - - x = np.linspace(0, 5, 200) - y = ref_sin(x) - - sp_sin.offset.fixed = False - sp_sin.phase.fixed = False - - f = Fitter(sp_sin, sp_sin) - f.switch_engine("bumps") - assert fit_method in f.available_methods() - result = f.fit(x, y, method=fit_method) - check_fit_results(result, sp_sin, ref_sin, x) - - -@pytest.mark.parametrize("fit_engine", ["lmfit", "bumps", "DFO_LS"]) -def test_fit_constraints(genObjs2, fit_engine): - ref_sin = genObjs2[0] - sp_sin = genObjs2[1] - - x = np.linspace(0, 5, 200) - y = ref_sin(x) - - sp_sin.phase.fixed = False - - f = Fitter(sp_sin, sp_sin) - - assert len(f.fit_constraints()) == 0 - c = ObjConstraint(sp_sin.offset, "2*", sp_sin.phase) - f.add_fit_constraint(c) - - if fit_engine is not None: - try: - f.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - - result = f.fit(x, y) - check_fit_results(result, sp_sin, ref_sin, x) - assert len(f.fit_constraints()) == 1 - f.remove_fit_constraint(0) - assert len(f.fit_constraints()) == 0 - - -# def test_fit_makeModel(genObjs): -# ref_sin = genObjs[0] -# sp_sin = genObjs[1] -# -# x = np.linspace(0, 5, 200) -# y = ref_sin(x) -# -# sp_sin.offset.fixed = False -# sp_sin.phase.fixed = False -# -# f = Fitter(sp_sin, sp_sin) -# model = f.make_model() -# result = f.fit(x, y, model=model) -# check_fit_results(result, sp_sin, ref_sin, x) - - -@pytest.mark.parametrize("with_errors", [False, True]) -@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "DFO_LS"]) -def test_multi_fit(genObjs, genObjs2, fit_engine, with_errors): - ref_sin1 = genObjs[0] - ref_sin2 = genObjs2[0] - - ref_sin1.offset.user_constraints["ref_sin2"] = ObjConstraint( - ref_sin2.offset, "", ref_sin1.offset - ) - ref_sin1.offset.user_constraints["ref_sin2"]() - - sp_sin1 = genObjs[1] - sp_sin2 = genObjs2[1] - - sp_sin1.offset.user_constraints["sp_sin2"] = ObjConstraint( - sp_sin2.offset, "", sp_sin1.offset - ) - sp_sin1.offset.user_constraints["sp_sin2"]() - - x1 = np.linspace(0, 5, 200) - y1 = ref_sin1(x1) - x2 = np.copy(x1) - y2 = ref_sin2(x2) - - sp_sin1.offset.fixed = False - sp_sin1.phase.fixed = False - sp_sin2.phase.fixed = False - - f = MultiFitter([sp_sin1, sp_sin2], [sp_sin1, sp_sin2]) - if fit_engine is not None: - try: - f.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - - args = [[x1, x2], [y1, y2]] - kwargs = {} - if with_errors: - kwargs["weights"] = [1 / np.sqrt(y1), 1 / np.sqrt(y2)] - - results = f.fit(*args, **kwargs) - X = [x1, x2] - Y = [y1, y2] - F_ref = [ref_sin1, ref_sin2] - F_real = [sp_sin1, sp_sin2] - for idx, result in enumerate(results): - assert result.n_pars == len(sp_sin1.get_fit_parameters()) + len( - sp_sin2.get_fit_parameters() - ) - assert result.chi2 == pytest.approx( - 0, abs=1.5e-3 * (len(result.x) - result.n_pars) - ) - assert result.reduced_chi == pytest.approx(0, abs=1.5e-3) - assert result.success - assert np.all(result.x == X[idx]) - assert np.all(result.y_obs == Y[idx]) - assert result.y_calc == pytest.approx(F_ref[idx](X[idx]), abs=1e-2) - assert result.residual == pytest.approx( - F_real[idx](X[idx]) - F_ref[idx](X[idx]), abs=1e-2 - ) - - -@pytest.mark.parametrize("with_errors", [False, True]) -@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "DFO_LS"]) -def test_multi_fit2(genObjs, genObjs2, fit_engine, with_errors): - ref_sin1 = genObjs[0] - ref_sin2 = genObjs2[0] - ref_line = Line.from_pars(1, 4.6) - - ref_sin1.offset.user_constraints["ref_sin2"] = ObjConstraint( - ref_sin2.offset, "", ref_sin1.offset - ) - ref_sin1.offset.user_constraints["ref_line"] = ObjConstraint( - ref_line.m, "", ref_sin1.offset - ) - ref_sin1.offset.user_constraints["ref_sin2"]() - ref_sin1.offset.user_constraints["ref_line"]() - - sp_sin1 = genObjs[1] - sp_sin2 = genObjs2[1] - sp_line = Line.from_pars(0.43, 6.1) - - sp_sin1.offset.user_constraints["sp_sin2"] = ObjConstraint( - sp_sin2.offset, "", sp_sin1.offset - ) - sp_sin1.offset.user_constraints["sp_line"] = ObjConstraint( - sp_line.m, "", sp_sin1.offset - ) - sp_sin1.offset.user_constraints["sp_sin2"]() - sp_sin1.offset.user_constraints["sp_line"]() - - x1 = np.linspace(0, 5, 200) - y1 = ref_sin1(x1) - x3 = np.copy(x1) - y3 = ref_sin2(x3) - x2 = np.copy(x1) - y2 = ref_line(x2) - - sp_sin1.offset.fixed = False - sp_sin1.phase.fixed = False - sp_sin2.phase.fixed = False - sp_line.c.fixed = False - - f = MultiFitter([sp_sin1, sp_line, sp_sin2], [sp_sin1, sp_line, sp_sin2]) - if fit_engine is not None: - try: - f.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - - args = [[x1, x2, x3], [y1, y2, y3]] - kwargs = {} - if with_errors: - kwargs["weights"] = [1 / np.sqrt(y1), 1 / np.sqrt(y2), 1 / np.sqrt(y3)] - - results = f.fit(*args, **kwargs) - X = [x1, x2, x3] - Y = [y1, y2, y3] - F_ref = [ref_sin1, ref_line, ref_sin2] - F_real = [sp_sin1, sp_line, sp_sin2] - - assert len(results) == len(X) - - for idx, result in enumerate(results): - assert result.n_pars == len(sp_sin1.get_fit_parameters()) + len( - sp_sin2.get_fit_parameters() - ) + len(sp_line.get_fit_parameters()) - assert result.chi2 == pytest.approx( - 0, abs=1.5e-3 * (len(result.x) - result.n_pars) - ) - assert result.reduced_chi == pytest.approx(0, abs=1.5e-3) - assert result.success - assert np.all(result.x == X[idx]) - assert np.all(result.y_obs == Y[idx]) - assert result.y_calc == pytest.approx(F_real[idx](X[idx]), abs=1e-2) - assert result.residual == pytest.approx( - F_ref[idx](X[idx]) - F_real[idx](X[idx]), abs=1e-2 - ) - - -class AbsSin2D(BaseObj): - def __init__(self, offset: Parameter, phase: Parameter): - super(AbsSin2D, self).__init__("sin2D", offset=offset, phase=phase) - - @classmethod - def from_pars(cls, offset, phase): - offset = Parameter("offset", offset) - phase = Parameter("phase", phase) - return cls(offset=offset, phase=phase) - - def __call__(self, x): - X = x[:, :, 0] - Y = x[:, :, 1] - return np.abs( - np.sin(self.phase.raw_value * X + self.offset.raw_value) - ) * np.abs(np.sin(self.phase.raw_value * Y + self.offset.raw_value)) - - -class AbsSin2DL(AbsSin2D): - def __call__(self, x): - X = x[:, 0] - Y = x[:, 1] - return np.abs( - np.sin(self.phase.raw_value * X + self.offset.raw_value) - ) * np.abs(np.sin(self.phase.raw_value * Y + self.offset.raw_value)) - - -@pytest.mark.parametrize("with_errors", [False, True]) -@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "DFO_LS"]) -def test_2D_vectorized(fit_engine, with_errors): - x = np.linspace(0, 5, 200) - mm = AbsSin2D.from_pars(0.3, 1.6) - m2 = AbsSin2D.from_pars( - 0.1, 1.8 - ) # The fit is quite sensitive to the initial values :-( - X, Y = np.meshgrid(x, x) - XY = np.stack((X, Y), axis=2) - ff = Fitter(m2, m2) - if fit_engine is not None: - try: - ff.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - try: - args = [XY, mm(XY)] - kwargs = {"vectorized": True} - if with_errors: - kwargs["weights"] = 1 / np.sqrt(args[1]) - result = ff.fit(*args, **kwargs) - except FitError as e: - if "Unable to allocate" in str(e): - pytest.skip(msg="MemoryError - Matrix too large") - else: - raise e - assert result.n_pars == len(m2.get_fit_parameters()) - assert result.reduced_chi == pytest.approx(0, abs=1.5e-3) - assert result.success - assert np.all(result.x == XY) - y_calc_ref = m2(XY) - assert result.y_calc == pytest.approx(y_calc_ref, abs=1e-2) - assert result.residual == pytest.approx(mm(XY) - y_calc_ref, abs=1e-2) - - -@pytest.mark.parametrize("with_errors", [False, True]) -@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "DFO_LS"]) -def test_2D_non_vectorized(fit_engine, with_errors): - x = np.linspace(0, 5, 200) - mm = AbsSin2DL.from_pars(0.3, 1.6) - m2 = AbsSin2DL.from_pars( - 0.1, 1.8 - ) # The fit is quite sensitive to the initial values :-( - X, Y = np.meshgrid(x, x) - XY = np.stack((X, Y), axis=2) - ff = Fitter(m2, m2) - if fit_engine is not None: - try: - ff.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - try: - args = [XY, mm(XY.reshape(-1, 2))] - kwargs = {"vectorized": False} - if with_errors: - kwargs["weights"] = 1 / np.sqrt(args[1]) - result = ff.fit(*args, **kwargs) - except FitError as e: - if "Unable to allocate" in str(e): - pytest.skip(msg="MemoryError - Matrix too large") - else: - raise e - assert result.n_pars == len(m2.get_fit_parameters()) - assert result.reduced_chi == pytest.approx(0, abs=1.5e-3) - assert result.success - assert np.all(result.x == XY) - y_calc_ref = m2(XY.reshape(-1, 2)) - assert result.y_calc == pytest.approx(y_calc_ref, abs=1e-2) - assert result.residual == pytest.approx( - mm(XY.reshape(-1, 2)) - y_calc_ref, abs=1e-2 - ) - - -@pytest.mark.parametrize("with_errors", [False, True]) -@pytest.mark.parametrize("fit_engine", [None, "lmfit", "bumps", "DFO_LS"]) -def test_multi_fit_1D_2D(genObjs, fit_engine, with_errors): - # Generate fit and reference objects - ref_sin1D = genObjs[0] - sp_sin1D = genObjs[1] - - ref_sin2D = AbsSin2D.from_pars(0.3, 1.6) - sp_sin2D = AbsSin2D.from_pars( - 0.1, 1.75 - ) # The fit is VERY sensitive to the initial values :-( - - # Link the parameters - ref_sin1D.offset.user_constraints["ref_sin2"] = ObjConstraint( - ref_sin2D.offset, "", ref_sin1D.offset - ) - ref_sin1D.offset.user_constraints["ref_sin2"]() - - sp_sin1D.offset.user_constraints["sp_sin2"] = ObjConstraint( - sp_sin2D.offset, "", sp_sin1D.offset - ) - sp_sin1D.offset.user_constraints["sp_sin2"]() - - # Generate data - x1D = np.linspace(0.2, 3.8, 400) - y1D = ref_sin1D(x1D) - - x = np.linspace(0, 5, 200) - X, Y = np.meshgrid(x, x) - x2D = np.stack((X, Y), axis=2) - y2D = ref_sin2D(x2D) - - ff = MultiFitter([sp_sin1D, sp_sin2D], [sp_sin1D, sp_sin2D]) - if fit_engine is not None: - try: - ff.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - - sp_sin1D.offset.fixed = False - sp_sin1D.phase.fixed = False - sp_sin2D.phase.fixed = False - - f = MultiFitter([sp_sin1D, sp_sin2D], [sp_sin1D, sp_sin2D]) - if fit_engine is not None: - try: - f.switch_engine(fit_engine) - except AttributeError: - pytest.skip(msg=f"{fit_engine} is not installed") - try: - args = [[x1D, x2D], [y1D, y2D]] - kwargs = {"vectorized": True} - if with_errors: - kwargs["weights"] = [1 / np.sqrt(y1D), 1 / np.sqrt(y2D)] - results = f.fit(*args, **kwargs) - except FitError as e: - if "Unable to allocate" in str(e): - pytest.skip(msg="MemoryError - Matrix too large") - else: - raise e - - X = [x1D, x2D] - Y = [y1D, y2D] - F_ref = [ref_sin1D, ref_sin2D] - F_real = [sp_sin1D, sp_sin2D] - for idx, result in enumerate(results): - assert result.n_pars == len(sp_sin1D.get_fit_parameters()) + len( - sp_sin2D.get_fit_parameters() - ) - assert result.chi2 == pytest.approx( - 0, abs=1.5e-3 * (len(result.x) - result.n_pars) - ) - assert result.reduced_chi == pytest.approx(0, abs=1.5e-3) - assert result.success - assert np.all(result.x == X[idx]) - assert np.all(result.y_obs == Y[idx]) - assert result.y_calc == pytest.approx(F_ref[idx](X[idx]), abs=1e-2) - assert result.residual == pytest.approx( - F_real[idx](X[idx]) - F_ref[idx](X[idx]), abs=1e-2 - )