From 347a5b772aaebbe85dda7d0dbcd0404f2f6efe7f Mon Sep 17 00:00:00 2001 From: Martin Reinecke <martin@mpa-garching.mpg.de> Date: Fri, 26 Jan 2018 10:28:01 +0100 Subject: [PATCH] more scipy work --- ci/install_basics.sh | 4 +- nifty4/__init__.py | 1 + nifty4/data_objects/distributed_do.py | 7 +- nifty4/data_objects/numpy_do.py | 4 + nifty4/dobj.py | 2 +- nifty4/minimization/scipy_minimizer.py | 109 ++++++++++++++++++++++ test/test_minimization/test_minimizers.py | 25 +++-- 7 files changed, 140 insertions(+), 12 deletions(-) create mode 100644 nifty4/minimization/scipy_minimizer.py diff --git a/ci/install_basics.sh b/ci/install_basics.sh index ba35c23d4..d58cbbe55 100644 --- a/ci/install_basics.sh +++ b/ci/install_basics.sh @@ -1,3 +1,3 @@ apt-get install -y build-essential git autoconf libtool pkg-config libfftw3-dev openmpi-bin libopenmpi-dev \ - python python-pip python-dev python-nose python-numpy python-matplotlib python-future python-mpi4py \ - python3 python3-pip python3-dev python3-nose python3-numpy python3-matplotlib python3-future python3-mpi4py + python python-pip python-dev python-nose python-numpy python-matplotlib python-future python-mpi4py python-scipy \ + python3 python3-pip python3-dev python3-nose python3-numpy python3-matplotlib python3-future python3-mpi4py python3-scipy diff --git a/nifty4/__init__.py b/nifty4/__init__.py index c0df5e303..0cb58e806 100644 --- a/nifty4/__init__.py +++ b/nifty4/__init__.py @@ -44,6 +44,7 @@ from .minimization.descent_minimizer import DescentMinimizer from .minimization.steepest_descent import SteepestDescent from .minimization.vl_bfgs import VL_BFGS from .minimization.relaxed_newton import RelaxedNewton +from .minimization.scipy_minimizer import NewtonCG, L_BFGS_B from .minimization.energy import Energy from .minimization.quadratic_energy import QuadraticEnergy from .minimization.line_energy import LineEnergy diff --git a/nifty4/data_objects/distributed_do.py b/nifty4/data_objects/distributed_do.py index c568207f9..c4dbd3172 100644 --- a/nifty4/data_objects/distributed_do.py +++ b/nifty4/data_objects/distributed_do.py @@ -27,6 +27,10 @@ rank = _comm.Get_rank() master = (rank == 0) +def is_numpy(): + return False + + def mprint(*args): if master: print(*args) @@ -259,6 +263,7 @@ class data_object(object): def fill(self, value): self._data.fill(value) + def full(shape, fill_value, dtype=None, distaxis=0): return data_object(shape, np.full(local_shape(shape, distaxis), fill_value, dtype), distaxis) @@ -346,7 +351,7 @@ def local_data(arr): def ibegin_from_shape(glob_shape, distaxis=0): res = [0] * len(glob_shape) - if distaxis<0: + if distaxis < 0: return res res[distaxis] = _shareRange(glob_shape[distaxis], ntask, rank)[0] return tuple(res) diff --git a/nifty4/data_objects/numpy_do.py b/nifty4/data_objects/numpy_do.py index 189ca40f0..7eeb0eaa1 100644 --- a/nifty4/data_objects/numpy_do.py +++ b/nifty4/data_objects/numpy_do.py @@ -30,6 +30,10 @@ rank = 0 master = True +def is_numpy(): + return True + + def mprint(*args): print(*args) diff --git a/nifty4/dobj.py b/nifty4/dobj.py index 2aa0eebf9..4d1375970 100644 --- a/nifty4/dobj.py +++ b/nifty4/dobj.py @@ -33,4 +33,4 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", "log", "tanh", "sqrt", "from_object", "from_random", "local_data", "ibegin", "ibegin_from_shape", "np_allreduce_sum", "distaxis", "from_local_data", "from_global_data", "to_global_data", - "redistribute", "default_distaxis", "mprint"] + "redistribute", "default_distaxis", "mprint", "is_numpy"] diff --git a/nifty4/minimization/scipy_minimizer.py b/nifty4/minimization/scipy_minimizer.py new file mode 100644 index 000000000..86065c582 --- /dev/null +++ b/nifty4/minimization/scipy_minimizer.py @@ -0,0 +1,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) diff --git a/test/test_minimization/test_minimizers.py b/test/test_minimization/test_minimizers.py index b7e5b1af0..30b54e334 100644 --- a/test/test_minimization/test_minimizers.py +++ b/test/test_minimization/test_minimizers.py @@ -18,14 +18,16 @@ import unittest import numpy as np -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_equal import nifty4 as ift from itertools import product from test.common import expand +from nose.plugins.skip import SkipTest spaces = [ift.RGSpace([1024], distances=0.123), ift.HPSpace(32)] minimizers = [ift.SteepestDescent, ift.RelaxedNewton, ift.VL_BFGS, - ift.ConjugateGradient, ift.NonlinearCG] + ift.ConjugateGradient, ift.NonlinearCG, + ift.NewtonCG, ift.L_BFGS_B] class Test_Minimizers(unittest.TestCase): @@ -39,13 +41,20 @@ class Test_Minimizers(unittest.TestCase): covariance = ift.DiagonalOperator(covariance_diagonal) required_result = ift.Field.ones(space, dtype=np.float64) - IC = ift.GradientNormController(tol_abs_gradnorm=1e-5) - minimizer = minimizer_class(controller=IC) - energy = ift.QuadraticEnergy(A=covariance, b=required_result, - position=starting_point) + IC = ift.GradientNormController(verbose=True,tol_abs_gradnorm=1e-5, iteration_limit=1000) + try: + minimizer = minimizer_class(controller=IC) + energy = ift.QuadraticEnergy(A=covariance, b=required_result, + position=starting_point) - (energy, convergence) = minimizer(energy) - assert convergence == IC.CONVERGED + (energy, convergence) = minimizer(energy) + except NotImplementedError: + raise SkipTest + + assert_equal(convergence, IC.CONVERGED) assert_allclose(ift.dobj.to_global_data(energy.position.val), 1./ift.dobj.to_global_data(covariance_diagonal.val), rtol=1e-3, atol=1e-3) + + +#MR FIXME: add Rosenbrock test -- GitLab