scipy_minimizer.py 3.99 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program 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 General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

from __future__ import division
from .minimizer import Minimizer
from ..field import Field
from .. import dobj


class ScipyMinimizer(Minimizer):
    """Scipy-based minimizer

    Parameters
    ----------
    controller : IterationController
        Object that decides when to terminate the minimization.
    method     : str
        The selected Scipy minimization method.
    options    : dictionary
        A set of custom options for the selected minimizer.
    """

    def __init__(self, controller, method, options, need_hessp):
        super(ScipyMinimizer, self).__init__()
        if not dobj.is_numpy():
            raise NotImplementedError
        self._controller = controller
        self._method = method
        self._options = options
        self._need_hessp = need_hessp

    def __call__(self, energy):
        class _MinimizationDone:
            pass

        class _MinHelper(object):
            def __init__(self, controller, energy):
                self._controller = controller
                self._energy = energy
                self._domain = energy.position.domain

            def _update(self, x):
                pos = Field(self._domain, x.reshape(self._domain.shape))
                if (pos.val != self._energy.position.val).any():
                    self._energy = self._energy.at(pos)
                    status = self._controller.check(self._energy)
                    if status != self._controller.CONTINUE:
                        raise _MinimizationDone

            def fun(self, x):
                self._update(x)
                return self._energy.value

            def jac(self, x):
                self._update(x)
                return self._energy.gradient.val.reshape(-1)

            def hessp(self, x, p):
                self._update(x)
                vec = Field(self._domain, p.reshape(self._domain.shape))
                res = self._energy.curvature(vec)
                return res.val.reshape(-1)

        import scipy.optimize as opt
        hlp = _MinHelper(self._controller, energy)
        energy = None
        status = self._controller.start(hlp._energy)
        if status != self._controller.CONTINUE:
            return hlp._energy, status
        try:
            if self._need_hessp:
                opt.minimize(hlp.fun, hlp._energy.position.val.reshape(-1),
                             method=self._method, jac=hlp.jac,
                             hessp=hlp.hessp,
                             options=self._options)
            else:
                opt.minimize(hlp.fun, hlp._energy.position.val.reshape(-1),
                             method=self._method, jac=hlp.jac,
                             options=self._options)
        except _MinimizationDone:
            status = self._controller.check(hlp._energy)
            return hlp._energy, self._controller.check(hlp._energy)
        return hlp._energy, self._controller.ERROR


def NewtonCG(controller):
    return ScipyMinimizer(controller, "Newton-CG",
                          {"xtol": 1e-20, "maxiter": None}, True)


def L_BFGS_B(controller, maxcor=10):
    return ScipyMinimizer(controller, "L-BFGS-B",
                          {"ftol": 1e-20, "gtol": 1e-20, "maxcor": maxcor},
                          False)