test__dual_annealing.py 10.3 KB
Newer Older
xuebingbing's avatar
xuebingbing 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 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
# Dual annealing unit tests implementation.
# Copyright (c) 2018 Sylvain Gubian <sylvain.gubian@pmi.com>,
# Yang Xiang <yang.xiang@pmi.com>
# Author: Sylvain Gubian, PMP S.A.
"""
Unit tests for the dual annealing global optimizer
"""
from scipy.optimize import dual_annealing
from scipy.optimize._dual_annealing import VisitingDistribution
from scipy.optimize._dual_annealing import ObjectiveFunWrapper
from scipy.optimize._dual_annealing import EnergyState
from scipy.optimize._dual_annealing import LocalSearchWrapper
from scipy.optimize import rosen, rosen_der
import pytest
import numpy as np
from numpy.testing import assert_equal, assert_allclose, assert_array_less
from pytest import raises as assert_raises
from scipy._lib._util import check_random_state
from scipy._lib._pep440 import Version


class TestDualAnnealing:

    def setup_method(self):
        # A function that returns always infinity for initialization tests
        self.weirdfunc = lambda x: np.inf
        # 2-D bounds for testing function
        self.ld_bounds = [(-5.12, 5.12)] * 2
        # 4-D bounds for testing function
        self.hd_bounds = self.ld_bounds * 4
        # Number of values to be generated for testing visit function
        self.nbtestvalues = 5000
        self.high_temperature = 5230
        self.low_temperature = 0.1
        self.qv = 2.62
        self.seed = 1234
        self.rs = check_random_state(self.seed)
        self.nb_fun_call = 0
        self.ngev = 0

    def callback(self, x, f, context):
        # For testing callback mechanism. Should stop for e <= 1 as
        # the callback function returns True
        if f <= 1.0:
            return True

    def func(self, x, args=()):
        # Using Rastrigin function for performing tests
        if args:
            shift = args
        else:
            shift = 0
        y = np.sum((x - shift) ** 2 - 10 * np.cos(2 * np.pi * (
            x - shift))) + 10 * np.size(x) + shift
        self.nb_fun_call += 1
        return y

    def rosen_der_wrapper(self, x, args=()):
        self.ngev += 1
        return rosen_der(x, *args)

    def test_visiting_stepping(self):
        lu = list(zip(*self.ld_bounds))
        lower = np.array(lu[0])
        upper = np.array(lu[1])
        dim = lower.size
        vd = VisitingDistribution(lower, upper, self.qv, self.rs)
        values = np.zeros(dim)
        x_step_low = vd.visiting(values, 0, self.high_temperature)
        # Make sure that only the first component is changed
        assert_equal(np.not_equal(x_step_low, 0), True)
        values = np.zeros(dim)
        x_step_high = vd.visiting(values, dim, self.high_temperature)
        # Make sure that component other than at dim has changed
        assert_equal(np.not_equal(x_step_high[0], 0), True)

    def test_visiting_dist_high_temperature(self):
        lu = list(zip(*self.ld_bounds))
        lower = np.array(lu[0])
        upper = np.array(lu[1])
        vd = VisitingDistribution(lower, upper, self.qv, self.rs)
        # values = np.zeros(self.nbtestvalues)
        # for i in np.arange(self.nbtestvalues):
        #     values[i] = vd.visit_fn(self.high_temperature)
        values = vd.visit_fn(self.high_temperature, self.nbtestvalues)

        # Visiting distribution is a distorted version of Cauchy-Lorentz
        # distribution, and as no 1st and higher moments (no mean defined,
        # no variance defined).
        # Check that big tails values are generated
        assert_array_less(np.min(values), 1e-10)
        assert_array_less(1e+10, np.max(values))

    def test_reset(self):
        owf = ObjectiveFunWrapper(self.weirdfunc)
        lu = list(zip(*self.ld_bounds))
        lower = np.array(lu[0])
        upper = np.array(lu[1])
        es = EnergyState(lower, upper)
        assert_raises(ValueError, es.reset, owf, check_random_state(None))

    def test_low_dim(self):
        ret = dual_annealing(
            self.func, self.ld_bounds, seed=self.seed)
        assert_allclose(ret.fun, 0., atol=1e-12)
        assert ret.success

    def test_high_dim(self):
        ret = dual_annealing(self.func, self.hd_bounds, seed=self.seed)
        assert_allclose(ret.fun, 0., atol=1e-12)
        assert ret.success

    def test_low_dim_no_ls(self):
        ret = dual_annealing(self.func, self.ld_bounds,
                             no_local_search=True, seed=self.seed)
        assert_allclose(ret.fun, 0., atol=1e-4)

    def test_high_dim_no_ls(self):
        ret = dual_annealing(self.func, self.hd_bounds,
                             no_local_search=True, seed=self.seed)
        assert_allclose(ret.fun, 0., atol=1e-4)

    def test_nb_fun_call(self):
        ret = dual_annealing(self.func, self.ld_bounds, seed=self.seed)
        assert_equal(self.nb_fun_call, ret.nfev)

    def test_nb_fun_call_no_ls(self):
        ret = dual_annealing(self.func, self.ld_bounds,
                             no_local_search=True, seed=self.seed)
        assert_equal(self.nb_fun_call, ret.nfev)

    def test_max_reinit(self):
        assert_raises(ValueError, dual_annealing, self.weirdfunc,
                      self.ld_bounds)

    def test_reproduce(self):
        res1 = dual_annealing(self.func, self.ld_bounds, seed=self.seed)
        res2 = dual_annealing(self.func, self.ld_bounds, seed=self.seed)
        res3 = dual_annealing(self.func, self.ld_bounds, seed=self.seed)
        # If we have reproducible results, x components found has to
        # be exactly the same, which is not the case with no seeding
        assert_equal(res1.x, res2.x)
        assert_equal(res1.x, res3.x)

    @pytest.mark.skipif(Version(np.__version__) < Version('1.17'),
                        reason='Generator not available for numpy, < 1.17')
    def test_rand_gen(self):
        # check that np.random.Generator can be used (numpy >= 1.17)
        # obtain a np.random.Generator object
        rng = np.random.default_rng(1)

        res1 = dual_annealing(self.func, self.ld_bounds, seed=rng)
        # seed again
        rng = np.random.default_rng(1)
        res2 = dual_annealing(self.func, self.ld_bounds, seed=rng)
        # If we have reproducible results, x components found has to
        # be exactly the same, which is not the case with no seeding
        assert_equal(res1.x, res2.x)

    def test_bounds_integrity(self):
        wrong_bounds = [(-5.12, 5.12), (1, 0), (5.12, 5.12)]
        assert_raises(ValueError, dual_annealing, self.func,
                      wrong_bounds)

    def test_bound_validity(self):
        invalid_bounds = [(-5, 5), (-np.inf, 0), (-5, 5)]
        assert_raises(ValueError, dual_annealing, self.func,
                      invalid_bounds)
        invalid_bounds = [(-5, 5), (0, np.inf), (-5, 5)]
        assert_raises(ValueError, dual_annealing, self.func,
                      invalid_bounds)
        invalid_bounds = [(-5, 5), (0, np.nan), (-5, 5)]
        assert_raises(ValueError, dual_annealing, self.func,
                      invalid_bounds)

    def test_max_fun_ls(self):
        ret = dual_annealing(self.func, self.ld_bounds, maxfun=100,
                             seed=self.seed)

        ls_max_iter = min(max(
            len(self.ld_bounds) * LocalSearchWrapper.LS_MAXITER_RATIO,
            LocalSearchWrapper.LS_MAXITER_MIN),
            LocalSearchWrapper.LS_MAXITER_MAX)
        assert ret.nfev <= 100 + ls_max_iter
        assert not ret.success

    def test_max_fun_no_ls(self):
        ret = dual_annealing(self.func, self.ld_bounds,
                             no_local_search=True, maxfun=500, seed=self.seed)
        assert ret.nfev <= 500
        assert not ret.success

    def test_maxiter(self):
        ret = dual_annealing(self.func, self.ld_bounds, maxiter=700,
                             seed=self.seed)
        assert ret.nit <= 700

    # Testing that args are passed correctly for dual_annealing
    def test_fun_args_ls(self):
        ret = dual_annealing(self.func, self.ld_bounds,
                             args=((3.14159,)), seed=self.seed)
        assert_allclose(ret.fun, 3.14159, atol=1e-6)

    # Testing that args are passed correctly for pure simulated annealing
    def test_fun_args_no_ls(self):
        ret = dual_annealing(self.func, self.ld_bounds,
                             args=((3.14159, )), no_local_search=True,
                             seed=self.seed)
        assert_allclose(ret.fun, 3.14159, atol=1e-4)

    def test_callback_stop(self):
        # Testing that callback make the algorithm stop for
        # fun value <= 1.0 (see callback method)
        ret = dual_annealing(self.func, self.ld_bounds,
                             callback=self.callback, seed=self.seed)
        assert ret.fun <= 1.0
        assert 'stop early' in ret.message[0]
        assert not ret.success

    @pytest.mark.parametrize('method, atol', [
        ('Nelder-Mead', 2e-5),
        ('COBYLA', 1e-5),
        ('Powell', 1e-8),
        ('CG', 1e-8),
        ('BFGS', 1e-8),
        ('TNC', 1e-8),
        ('SLSQP', 2e-7),
    ])
    def test_multi_ls_minimizer(self, method, atol):
        ret = dual_annealing(self.func, self.ld_bounds,
                             local_search_options=dict(method=method),
                             seed=self.seed)
        assert_allclose(ret.fun, 0., atol=atol)

    def test_wrong_restart_temp(self):
        assert_raises(ValueError, dual_annealing, self.func,
                      self.ld_bounds, restart_temp_ratio=1)
        assert_raises(ValueError, dual_annealing, self.func,
                      self.ld_bounds, restart_temp_ratio=0)

    def test_gradient_gnev(self):
        minimizer_opts = {
            'jac': self.rosen_der_wrapper,
        }
        ret = dual_annealing(rosen, self.ld_bounds,
                             local_search_options=minimizer_opts,
                             seed=self.seed)
        assert ret.njev == self.ngev

    def test_from_docstring(self):
        func = lambda x: np.sum(x * x - 10 * np.cos(2 * np.pi * x)) + 10 * np.size(x)
        lw = [-5.12] * 10
        up = [5.12] * 10
        ret = dual_annealing(func, bounds=list(zip(lw, up)), seed=1234)
        assert_allclose(ret.x,
                        [-4.26437714e-09, -3.91699361e-09, -1.86149218e-09,
                         -3.97165720e-09, -6.29151648e-09, -6.53145322e-09,
                         -3.93616815e-09, -6.55623025e-09, -6.05775280e-09,
                         -5.00668935e-09], atol=4e-8)
        assert_allclose(ret.fun, 0.000000, atol=5e-13)