From d602bea5b721bd5ff826bb8380613841b25759a4 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Wed, 2 Feb 2022 10:05:48 +0000
Subject: [PATCH 1/2] refactor: Classification of fit method in
 fits.least_squares simplified, precision of imiunit solver adjusted,
 prefitting for alternative methods removed.

---
 pyerrors/fits.py | 23 +++++++++--------------
 1 file changed, 9 insertions(+), 14 deletions(-)

diff --git a/pyerrors/fits.py b/pyerrors/fits.py
index 311fba89..ea51f6f0 100644
--- a/pyerrors/fits.py
+++ b/pyerrors/fits.py
@@ -96,7 +96,7 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
     initial_guess : list
         can provide an initial guess for the input parameters. Relevant for
                      non-linear fits with many parameters.
-    method : str
+    method : str, optional
         can be used to choose an alternative method for the minimization of chisquare.
         The possible methods are the ones which can be used for scipy.optimize.minimize and
         migrad of iminuit. If no method is specified, Levenberg-Marquard is used.
@@ -487,26 +487,21 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
             chisq = anp.sum(((y_f - model) / dy_f) ** 2)
             return chisq
 
-    if 'method' in kwargs and not (kwargs.get('method', 'Levenberg-Marquardt') == 'Levenberg-Marquardt'):
-        output.method = kwargs.get('method')
-        if not silent:
-            print('Method:', kwargs.get('method'))
-        if kwargs.get('method') == 'migrad':
-            fit_result = iminuit.minimize(chisqfunc, x0)
-            fit_result = iminuit.minimize(chisqfunc, fit_result.x)
+    output.method = kwargs.get('method', 'Levenberg-Marquardt')
+    if not silent:
+        print('Method:', output.method)
+
+    if output.method != 'Levenberg-Marquardt':
+        if output.method == 'migrad':
+            fit_result = iminuit.minimize(chisqfunc, x0, tol=1e-4)  # Stopping crieterion 0.002 * tol * errordef
             output.iterations = fit_result.nfev
         else:
-            fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'))
-            fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=1e-12)
+            fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=1e-12)
             output.iterations = fit_result.nit
 
         chisquare = fit_result.fun
 
     else:
-        output.method = 'Levenberg-Marquardt'
-        if not silent:
-            print('Method: Levenberg-Marquardt')
-
         if kwargs.get('correlated_fit') is True:
             def chisqfunc_residuals(p):
                 model = func(p, x)

From 952c7acef51da03029e515081da046dc1971e4b1 Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Wed, 2 Feb 2022 10:10:20 +0000
Subject: [PATCH 2/2] tests: additional tests for alternative fit methods in
 fits.least_squares added

---
 tests/fits_test.py | 29 +++++++++++++++++++++++++++--
 1 file changed, 27 insertions(+), 2 deletions(-)

diff --git a/tests/fits_test.py b/tests/fits_test.py
index 33a3b958..90a073e7 100644
--- a/tests/fits_test.py
+++ b/tests/fits_test.py
@@ -49,8 +49,6 @@ def test_least_squares():
         y = a[0] * np.exp(-a[1] * x)
         return y
 
-    out = pe.least_squares(x, oy, func, method='migrad')
-    out = pe.least_squares(x, oy, func, method='Powell')
     out = pe.least_squares(x, oy, func, expected_chisquare=True, resplot=True, qqplot=True)
     beta = out.fit_parameters
 
@@ -86,6 +84,33 @@ def test_least_squares():
     assert math.isclose(pe.covariance(betac[0], betac[1]), pcov[0, 1], abs_tol=1e-3)
 
 
+def test_alternative_solvers():
+    dim = 192
+    x = np.arange(dim)
+    y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim)
+    yerr = 0.1 + 0.1 * np.random.rand(dim)
+
+    oy = []
+    for i, item in enumerate(x):
+        oy.append(pe.pseudo_Obs(y[i], yerr[i], 'test'))
+
+    def func(a, x):
+        y = a[0] * np.exp(-a[1] * x)
+        return y
+
+    chisquare_values = []
+    out = pe.least_squares(x, oy, func, method='migrad')
+    chisquare_values.append(out.chisquare)
+    out = pe.least_squares(x, oy, func, method='Powell')
+    chisquare_values.append(out.chisquare)
+    out = pe.least_squares(x, oy, func, method='Nelder-Mead')
+    chisquare_values.append(out.chisquare)
+    out = pe.least_squares(x, oy, func, method='Levenberg-Marquardt')
+    chisquare_values.append(out.chisquare)
+    chisquare_values = np.array(chisquare_values)
+    assert np.all(np.isclose(chisquare_values, chisquare_values[0]))
+
+
 def test_correlated_fit():
     num_samples = 400
     N = 10