New upstream version 0.9.14

1 files changed, 397 insertions, 0 deletions

diff --git a/examples/example_brute.py b/examples/example_brute.py
new file mode 100644
index 0000000..5d2266e
--- /dev/null
+++ b/examples/example_brute.py
@@ -0,0 +1,397 @@
+"""
+Global minimization using the ``brute`` method (a.k.a. grid search)
+===================================================================
+
+"""
+###############################################################################
+# This notebook shows a simple example of using ``lmfit.minimize.brute`` that
+# uses the method with the same name from ``scipy.optimize``.
+#
+# The method computes the function’s value at each point of a multidimensional
+# grid of points, to find the global minimum of the function. It behaves
+# identically to ``scipy.optimize.brute`` in case finite bounds are given on
+# all varying parameters, but will also deal with non-bounded parameters
+# (see below).
+import copy
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lmfit import fit_report, Minimizer, Parameters
+
+###############################################################################
+# Let's start with the example given in the documentation of SciPy:
+#
+# "We illustrate the use of brute to seek the global minimum of a function of
+# two variables that is given as the sum of a positive-definite quadratic and
+# two deep “Gaussian-shaped” craters. Specifically, define the objective
+# function f as the sum of three other functions, ``f = f1 + f2 + f3``. We
+# suppose each of these has a signature ``(z, *params), where z = (x, y)``,
+# and params and the functions are as defined below."
+#
+# First, we create a set of Parameters where all variables except ``x`` and
+# ``y`` are given fixed values.
+params = Parameters()
+params.add_many(
+        ('a', 2, False),
+        ('b', 3, False),
+        ('c', 7, False),
+        ('d', 8, False),
+        ('e', 9, False),
+        ('f', 10, False),
+        ('g', 44, False),
+        ('h', -1, False),
+        ('i', 2, False),
+        ('j', 26, False),
+        ('k', 1, False),
+        ('l', -2, False),
+        ('scale', 0.5, False),
+        ('x', 0.0, True),
+        ('y', 0.0, True))
+
+###############################################################################
+# Second, create the three functions and the objective function:
+
+
+def f1(p):
+    par = p.valuesdict()
+    return (par['a'] * par['x']**2 + par['b'] * par['x'] * par['y'] +
+            par['c'] * par['y']**2 + par['d']*par['x'] + par['e']*par['y'] + par['f'])
+
+
+def f2(p):
+    par = p.valuesdict()
+    return (-1.0*par['g']*np.exp(-((par['x']-par['h'])**2 +
+            (par['y']-par['i'])**2) / par['scale']))
+
+
+def f3(p):
+    par = p.valuesdict()
+    return (-1.0*par['j']*np.exp(-((par['x']-par['k'])**2 +
+            (par['y']-par['l'])**2) / par['scale']))
+
+
+def f(params):
+    return f1(params) + f2(params) + f3(params)
+
+
+###############################################################################
+# Just as in the documentation we will do a grid search between ``-4`` and
+# ``4`` and use a stepsize of ``0.25``. The bounds can be set as usual with
+# the ``min`` and ``max`` attributes, and the stepsize is set using
+# ``brute_step``.
+params['x'].set(min=-4, max=4, brute_step=0.25)
+params['y'].set(min=-4, max=4, brute_step=0.25)
+
+###############################################################################
+# Performing the actual grid search is done with:
+fitter = Minimizer(f, params)
+result = fitter.minimize(method='brute')
+
+###############################################################################
+# , which will increment ``x`` and ``y`` between ``-4`` in increments of
+# ``0.25`` until ``4`` (not inclusive).
+grid_x, grid_y = [np.unique(par.ravel()) for par in result.brute_grid]
+print(grid_x)
+
+###############################################################################
+# The objective function is evaluated on this grid, and the raw output from
+# ``scipy.optimize.brute`` is stored in the MinimizerResult as
+# ``brute_<parname>`` attributes. These attributes are:
+#
+# ``result.brute_x0`` -- A 1-D array containing the coordinates of a point at
+# which the objective function had its minimum value.
+print(result.brute_x0)
+
+###############################################################################
+# ``result.brute_fval`` -- Function value at the point x0.
+print(result.brute_fval)
+
+###############################################################################
+# ``result.brute_grid`` -- Representation of the evaluation grid. It has the
+# same length as x0.
+print(result.brute_grid)
+
+###############################################################################
+# ``result.brute_Jout`` -- Function values at each point of the evaluation
+# grid, i.e., Jout = func(\*grid).
+print(result.brute_Jout)
+
+###############################################################################
+# **Reassuringly, the obtained results are indentical to using the method in
+# SciPy directly!**
+
+###############################################################################
+# Example 2: fit of a decaying sine wave
+#
+# In this example, will explain some of the options ot the algorithm.
+#
+# We start off by generating some synthetic data with noise for a decaying
+# sine wave, define an objective function and create a Parameter set.
+x = np.linspace(0, 15, 301)
+np.random.seed(7)
+noise = np.random.normal(size=x.size, scale=0.2)
+data = (5. * np.sin(2*x - 0.1) * np.exp(-x*x*0.025) + noise)
+plt.plot(x, data, 'b')
+
+
+def fcn2min(params, x, data):
+    """Model decaying sine wave, subtract data."""
+    amp = params['amp']
+    shift = params['shift']
+    omega = params['omega']
+    decay = params['decay']
+    model = amp * np.sin(x*omega + shift) * np.exp(-x*x*decay)
+    return model - data
+
+
+# create a set of Parameters
+params = Parameters()
+params.add('amp', value=7, min=2.5)
+params.add('decay', value=0.05)
+params.add('shift', value=0.0, min=-np.pi/2., max=np.pi/2)
+params.add('omega', value=3, max=5)
+
+###############################################################################
+# In contrast to the implementation in SciPy (as shown in the first example),
+# varying parameters do not need to have finite bounds in lmfit. However, in
+# that case they **do** need the ``brute_step`` attribute specified, so let's
+# do that:
+params['amp'].set(brute_step=0.25)
+params['decay'].set(brute_step=0.005)
+params['omega'].set(brute_step=0.25)
+
+###############################################################################
+# Our initial parameter set is now defined as shown below and this will
+# determine how the grid is set-up.
+params.pretty_print()
+
+###############################################################################
+# First, we initialize a Minimizer and perform the grid search:
+fitter = Minimizer(fcn2min, params, fcn_args=(x, data))
+result_brute = fitter.minimize(method='brute', Ns=25, keep=25)
+
+###############################################################################
+# We used two new parameters here: ``Ns`` and ``keep``. The parameter ``Ns``
+# determines the \'number of grid points along the axes\' similarly to its usage
+# in SciPy. Together with ``brute_step``, ``min`` and ``max`` for a Parameter
+# it will dictate how the grid is set-up:
+#
+# **(1)** finite bounds are specified ("SciPy implementation"): uses
+# ``brute_step`` if present (in the example above) or uses ``Ns`` to generate
+# the grid. The latter scenario that interpolates ``Ns`` points from ``min``
+# to ``max`` (inclusive), is here shown for the parameter ``shift``:
+par_name = 'shift'
+indx_shift = result_brute.var_names.index(par_name)
+grid_shift = np.unique(result_brute.brute_grid[indx_shift].ravel())
+print("parameter = {}\nnumber of steps = {}\ngrid = {}".format(par_name,
+      len(grid_shift), grid_shift))
+
+###############################################################################
+# If finite bounds are not set for a certain parameter then the user **must**
+# specify ``brute_step`` - three more scenarios are considered here:
+#
+# **(2)** lower bound (min) and brute_step are specified:
+# range = (min, min + Ns * brute_step, brute_step)
+par_name = 'amp'
+indx_shift = result_brute.var_names.index(par_name)
+grid_shift = np.unique(result_brute.brute_grid[indx_shift].ravel())
+print("parameter = {}\nnumber of steps = {}\ngrid = {}".format(par_name, len(grid_shift), grid_shift))
+
+###############################################################################
+# **(3)** upper bound (max) and brute_step are specified:
+# range = (max - Ns * brute_step, max, brute_step)
+par_name = 'omega'
+indx_shift = result_brute.var_names.index(par_name)
+grid_shift = np.unique(result_brute.brute_grid[indx_shift].ravel())
+print("parameter = {}\nnumber of steps = {}\ngrid = {}".format(par_name, len(grid_shift), grid_shift))
+
+###############################################################################
+# **(4)** numerical value (value) and brute_step are specified:
+# range = (value - (Ns//2) * brute_step, value + (Ns//2) * brute_step, brute_step)
+par_name = 'decay'
+indx_shift = result_brute.var_names.index(par_name)
+grid_shift = np.unique(result_brute.brute_grid[indx_shift].ravel())
+print("parameter = {}\nnumber of steps = {}\ngrid = {}".format(par_name, len(grid_shift), grid_shift))
+
+###############################################################################
+# The ``MinimizerResult`` contains all the usual best-fit parameters and
+# fitting statistics. For example, the optimal solution from the grid search
+# is given below together with a plot:
+print(fit_report(result_brute))
+plt.plot(x, data, 'b')
+plt.plot(x, data + fcn2min(result_brute.params, x, data), 'r--')
+
+###############################################################################
+# We can see that this fit is already very good, which is what we should expect
+# since our ``brute`` force grid is sampled rather finely and encompasses the
+# "correct" values.
+#
+# In a more realistic, complicated example the ``brute`` method will be used
+# to get reasonable values for the parameters and perform another minimization
+# (e.g., using ``leastsq``) using those as starting values. That is where the
+# `keep`` parameter comes into play: it determines the "number of best
+# candidates from the brute force method that are stored in the ``candidates``
+# attribute". In the example above we store the best-ranking 25 solutions (the
+# default value is ``50`` and storing all the grid points can be accomplished
+# by choosing ``all``). The ``candidates`` attribute contains the parameters
+# and ``chisqr`` from the brute force method as a namedtuple,
+# ``(‘Candidate’, [‘params’, ‘score’])``, sorted on the (lowest) ``chisqr``
+# value. To access the values for a particular candidate one can use
+# ``result.candidate[#].params`` or ``result.candidate[#].score``, where a
+# lower # represents a better candidate. The ``show_candidates(#)`` uses the
+# ``pretty_print()`` method to show a specific candidate-# or all candidates
+# when no number is specified.
+#
+# The optimal fit is, as usual, stored in the ``MinimizerResult.params``
+# attribute and is, therefore, identical to ``result_brute.show_candidates(1)``.
+result_brute.show_candidates(1)
+
+###############################################################################
+#  In this case, the next-best scoring candidate has already a ``chisqr`` that
+#  increased quite a bit:
+result_brute.show_candidates(2)
+
+###############################################################################
+# and is, therefore, probably not so likely... However, as said above, in most
+# cases you'll want to do another minimization using the solutions from the
+# ``brute`` method as starting values. That can be easily accomplished as
+# shown in the code below, where we now perform a ``leastsq`` minimization
+# starting from the top-25 solutions and accept the solution if the ``chisqr``
+# is lower than the previously 'optimal' solution:
+best_result = copy.deepcopy(result_brute)
+
+for candidate in result_brute.candidates:
+    trial = fitter.minimize(method='leastsq', params=candidate.params)
+    if trial.chisqr < best_result.chisqr:
+        best_result = trial
+
+###############################################################################
+# From the ``leastsq`` minimization we obtain the following parameters for the
+# most optimal result:
+print(fit_report(best_result))
+
+###############################################################################
+# As expected the parameters have not changed significantly as they were
+# already very close to the "real" values, which can also be appreciated from
+# the plots below.
+plt.plot(x, data, 'b')
+plt.plot(x, data + fcn2min(result_brute.params, x, data), 'r--',
+         label='brute')
+plt.plot(x, data + fcn2min(best_result.params, x, data), 'g--',
+         label='brute followed by leastsq')
+plt.legend()
+
+
+###############################################################################
+# Finally, the results from the ``brute`` force grid-search can be visualized
+# using the rather lengthy Python function below (which might get incorporated
+# in lmfit at some point).
+def plot_results_brute(result, best_vals=True, varlabels=None,
+                       output=None):
+    """Visualize the result of the brute force grid search.
+
+    The output file will display the chi-square value per parameter and contour
+    plots for all combination of two parameters.
+
+    Inspired by the `corner` package (https://github.com/dfm/corner.py).
+
+    Parameters
+    ----------
+    result : :class:`~lmfit.minimizer.MinimizerResult`
+        Contains the results from the :meth:`brute` method.
+
+    best_vals : bool, optional
+        Whether to show the best values from the grid search (default is True).
+
+    varlabels : list, optional
+        If None (default), use `result.var_names` as axis labels, otherwise
+        use the names specified in `varlabels`.
+
+    output : str, optional
+        Name of the output PDF file (default is 'None')
+    """
+    from matplotlib.colors import LogNorm
+
+    npars = len(result.var_names)
+    fig, axes = plt.subplots(npars, npars)
+
+    if not varlabels:
+        varlabels = result.var_names
+    if best_vals and isinstance(best_vals, bool):
+        best_vals = result.params
+
+    for i, par1 in enumerate(result.var_names):
+        for j, par2 in enumerate(result.var_names):
+
+            # parameter vs chi2 in case of only one parameter
+            if npars == 1:
+                axes.plot(result.brute_grid, result.brute_Jout, 'o', ms=3)
+                axes.set_ylabel(r'$\chi^{2}$')
+                axes.set_xlabel(varlabels[i])
+                if best_vals:
+                    axes.axvline(best_vals[par1].value, ls='dashed', color='r')
+
+            # parameter vs chi2 profile on top
+            elif i == j and j < npars-1:
+                if i == 0:
+                    axes[0, 0].axis('off')
+                ax = axes[i, j+1]
+                red_axis = tuple([a for a in range(npars) if a != i])
+                ax.plot(np.unique(result.brute_grid[i]),
+                        np.minimum.reduce(result.brute_Jout, axis=red_axis),
+                        'o', ms=3)
+                ax.set_ylabel(r'$\chi^{2}$')
+                ax.yaxis.set_label_position("right")
+                ax.yaxis.set_ticks_position('right')
+                ax.set_xticks([])
+                if best_vals:
+                    ax.axvline(best_vals[par1].value, ls='dashed', color='r')
+
+            # parameter vs chi2 profile on the left
+            elif j == 0 and i > 0:
+                ax = axes[i, j]
+                red_axis = tuple([a for a in range(npars) if a != i])
+                ax.plot(np.minimum.reduce(result.brute_Jout, axis=red_axis),
+                        np.unique(result.brute_grid[i]), 'o', ms=3)
+                ax.invert_xaxis()
+                ax.set_ylabel(varlabels[i])
+                if i != npars-1:
+                    ax.set_xticks([])
+                elif i == npars-1:
+                    ax.set_xlabel(r'$\chi^{2}$')
+                if best_vals:
+                    ax.axhline(best_vals[par1].value, ls='dashed', color='r')
+
+            # contour plots for all combinations of two parameters
+            elif j > i:
+                ax = axes[j, i+1]
+                red_axis = tuple([a for a in range(npars) if a != i and a != j])
+                X, Y = np.meshgrid(np.unique(result.brute_grid[i]),
+                                   np.unique(result.brute_grid[j]))
+                lvls1 = np.linspace(result.brute_Jout.min(),
+                                    np.median(result.brute_Jout)/2.0, 7, dtype='int')
+                lvls2 = np.linspace(np.median(result.brute_Jout)/2.0,
+                                    np.median(result.brute_Jout), 3, dtype='int')
+                lvls = np.unique(np.concatenate((lvls1, lvls2)))
+                ax.contourf(X.T, Y.T, np.minimum.reduce(result.brute_Jout, axis=red_axis),
+                            lvls, norm=LogNorm())
+                ax.set_yticks([])
+                if best_vals:
+                    ax.axvline(best_vals[par1].value, ls='dashed', color='r')
+                    ax.axhline(best_vals[par2].value, ls='dashed', color='r')
+                    ax.plot(best_vals[par1].value, best_vals[par2].value, 'rs', ms=3)
+                if j != npars-1:
+                    ax.set_xticks([])
+                elif j == npars-1:
+                    ax.set_xlabel(varlabels[i])
+                if j - i >= 2:
+                    axes[i, j].axis('off')
+
+    if output is not None:
+        plt.savefig(output)
+
+
+###############################################################################
+# and finally, to generated the figure:
+plot_results_brute(result_brute, best_vals=True, varlabels=None)