r"""Test problems for testing non-linear solvers.
This module contains test problems with known solutions for testing
non-linear solvers.
"""
from __future__ import division
import numpy as np
from numpy import cos, sin
#from sympy import sin, cos
from mmf.objects import StateVars, process_vars
from mmf.objects import Computed
import mmf.math.multigrid.stencil
import mmf.solve.solver_interface as solver_interface
_MESGS = solver_interface.Mesg()
[docs]class TwoDimensional(object):
r"""These are problems of two variables `x` and `y`."""
[docs] def test(self):
"""Test the functions."""
import inspect
from mmf.math.differentiate import differentiate as diff
classes = [cls
for cls in self.__class__.__dict__.values()
if inspect.isclass(cls)]
for cls in classes:
print("Testing %s..." % cls.__name__)
c = cls()
x = 0.12
y = 0.34
J = c.J(x,y)
assert np.allclose(diff(lambda x:c.f(x,y), x, 1, 0.1),
c.f_x(x,y))
assert np.allclose(diff(lambda y:c.f(x,y), y, 1, 0.1),
c.f_y(x,y))
assert np.allclose(diff(lambda x:c.G([x,y])[0], x, 1, 0.1),
J[0,0])
assert np.allclose(diff(lambda y:c.G([x,y])[0], y, 1, 0.1),
J[0,1])
assert np.allclose(diff(lambda x:c.G([x,y])[1], x, 1, 0.1),
J[1,0])
assert np.allclose(diff(lambda y:c.G([x,y])[1], y, 1, 0.1),
J[1,1])
print("Testing %s. Passed." % cls.__name__)
[docs] class MexicanHat(object):
r"""Classical Mexican hat potential with global minimum at
$(0,0)$, saddle point at $(0, y_{+})$ and local maximum at
$(0, y_{-})$ where $2y_{\pm} = 3 \pm \sqrt{1 - b}$:
.. math::
f(x,y) = (x^2 + y^2 - 2y)^2 + \frac{bx^2 + by^2}{2}
The latter two features disappear for $b\geq 1$ while for
$a=b=0$ the hat is completely degenerate. If $a$ is
sufficiently large, then there are additional saddle points
off the $x=0$ axis.
This becomes a very difficult problem for efficient solution
when $a=0$ because the Jacobin becomes singular at the
minimum, spoiling the quadratic convergence. This is due to
the fact that function is quartic along the valley at this
point as can be seen by looking at the second derivative:
.. math::
f_{,xx} = 12x^2 + 4(y^2 - 2y) + a """
def __init__(self, a=0.5, b=0.5):
self.a = a
self.b = b
[docs] def f(self, x, y):
return (x*x + y*(y - 2))**2 + (self.a*x*x + self.b*y*y)/2
[docs] def f_x(self, x, y):
return 4*x*(x*x + y*(y - 2)) + self.a*x
[docs] def f_y(self, x, y):
return 4*(y - 1)*(x*x + y*(y - 2)) + self.b*y
[docs] def F(self, x, y):
return self.f_x(x,y)**2 + self.f_y(x,y)**2
[docs] def G(self, x_):
x, y = x_
return np.array([self.f_x(x,y), self.f_y(x,y)])
[docs] def J(self, x, y):
return np.array([
[12*x*x + 4*y*(y - 2) + self.a,
8*x*(y - 1)],
[8*x*(y - 1),
4*(x*x + y*(y - 2)) + 8*(y - 1)**2 + self.b]])
[docs] def Jinv(self, x, y):
return np.linalg.inv(self.J(x,y))
[docs] class Himmelblau(object):
r"""Himmelblau's problem 2. Find the minimum of
.. math::
f(x, y) = 100(y-x^2)^2 + (1-x)^2
starting from $x=-1.2$, $y=1$. Solution $x=y=1$.
"""
[docs] def f(self, x, y):
return 100*(y - x*x)**2 + (1-x)**2
[docs] def f_x(self, x, y):
return -400*x*(y - x*x) - 2*(1-x)
[docs] def f_y(self, x, y):
return 200*(y - x*x)
[docs] def F(self, x, y):
return self.f_x(x,y)**2 + self.f_y(x,y)**2
[docs] def G(self, x_):
x, y = x_
return np.array([self.f_x(x,y), self.f_y(x,y)])
[docs] def J(self, x, y):
return np.array([
[2 - 400*y + 1200*x*x, -400*x],
[-400*x, 200]])
[docs] def Jinv(self, x, y):
return np.linalg.inv(self.J(x,y))
[docs] class SensitiveParameter(object):
r"""In this potential, `x` is a "sensitive parameter" in the
sense that a good strategy is to define $g(x) =
\min_{y}f(x,y)$ and minimize this. Trying to minimize both
`x` and `y` simultaneously will fail unless close to the
solution.
The solution is `x=y=0`.
.. math::
f(x,y) = (a +\cos(bx))(y - cx^2)^2 + dx^2
with default parameters `a=1.1`, `b=200`, `c=0.1`, `d=0.01`.
"""
def __init__(self, a=1.1, b=100, c=0.1, d=0.01):
self.a = a
self.b = b
self.c = c
self.d = d
[docs] def f(self, x, y):
a, b, c, d = self.a, self.b, self.c, self.d
return (a + cos(b*x))*(y - c*x**2)**2 + d*x**2
[docs] def f_x(self, x, y):
a, b, c, d = self.a, self.b, self.c, self.d
return (2*d*x - b*(y - c*x**2)**2*sin(b*x)
- 4*c*x*(y - c*x**2)*(a + cos(b*x)))
[docs] def f_y(self, x, y):
a, b, c, d = self.a, self.b, self.c, self.d
return (2*y - 2*c*x**2)*(a + cos(b*x))
[docs] def F(self, x, y):
return self.f_x(x,y)**2 + self.f_y(x,y)**2
[docs] def G(self, x_):
x, y = x_
return np.array([self.f_x(x,y), self.f_y(x,y)])
[docs] def J(self, x, y):
a, b, c, d = self.a, self.b, self.c, self.d
return np.array([
[2*d - b**2*(y - c*x**2)**2*cos(b*x)
- 4*c*(y - c*x**2)*(a + cos(b*x))
+ 8*c**2*x**2*(a + cos(b*x))
+ 8*b*c*x*(y - c*x**2)*sin(b*x),
-b*(2*y - 2*c*x**2)*sin(b*x) - 4*c*x*(a + cos(b*x))],
[-b*(2*y - 2*c*x**2)*sin(b*x) - 4*c*x*(a + cos(b*x)),
2*a + 2*cos(b*x)]])
[docs] def Jinv(self, x, y):
return np.linalg.inv(self.J(x,y))
[docs]class MediumScale(object):
"""Medium scale root-finding problems."""
[docs] class GomexRuggiero(solver_interface.Problem):
r"""These are the problems defined in Computers & Mathematics
with Applications Volume 32, Issue 3, August 1996, Pages 1-13
"A numerical study on large-scale nonlinear solvers",
M. A. Gomes-Ruggiero, D. N. Kozakevich and J. M. Martinez.
doi:10.1016/0898-1221(96)00109-5.
These represent a discretization of a function $u:
[0,1]\times[0,1] \mapsto \field{R}$ at $N^2$ equally spaced
points.
The equations are the discretized form of
.. math::
G_c(u) = f(x,y)
where $G_c$ is one of the following:
.. math::
(0)\qquad G_{c}(u) &= -\nabla^2 u
+ c\frac{u^3}{1 + x^2 + y^2}\\
(1)\qquad G_{c}(u) &= -\nabla^2 u
+ c e^{u}\\
(2)\qquad G_{c}(u) &= -\nabla^2 u
+ c u(u_{,x} + u_{,y})
We generate the right-hand side $f(x,y)$ dynamically so that
the solution is exact.
"""
_state_vars = [
('problem', 0, "Problem number"),
('c', 0, r"Parameter $c$. Consider $\abs{c}\leq 2000$."),
('n', 6,
r"Number of points in discretization is $2^{n-1}n_0 + 1$."),
('n_0', 2,
r"Number of points in discretization is $2^{n-1}n_0 + 1$."),
('f', Computed, "Right-hand size"),
('A', Computed, "Discrete negative Laplacian operator."),
('x', Computed, "Abscissa `(x, y, z, ...)`"),
('u_exact', Computed, "Exact solution."),
('state0', Computed),
('mg', Computed, "Multigrid solver for Poisson equation."),
]
process_vars()
def __init__(self, *v, **kw):
d = 2
self.mg = mmf.math.multigrid.MG(n_0=self.n_0, n=self.n,
d=d, L=1)
self.x = self.mg.x()
self.A = self.mg.get_mats()[0][-1]
self.x = (self.x[0].ravel(), self.x[1].ravel())
self.u_exact = self.get_u()
self.f = self.get_f(self.u_exact)
self.state0 = solver_interface.BareState(
x_=0*self.u_exact)
def get_u(self):
r"""Return exact solution"""
x, y = self.x
u = 10*x*y*(1-x)*(1-y)*np.exp(x**4.5)
return u
def get_f(self, u):
x, y = self.x
if 0 == self.problem:
tmp = u**3/(1 + x**2 + y**2)
elif 1 == self.problem:
tmp = np.exp(u)
else:
raise NotImplementedError(
"Problem %i not define" % (self.problem,))
f = self.A*u + self.c*tmp
return f
############## IProblem interface.
def get_initial_state(self):
return self.state0.copy_from()
def F(self, x, iter_data=None, abs_tol=None, rel_tol=None,
mesgs=_MESGS):
"""Fixed-point function."""
rhs = self.f - self.get_f(u=0*x)
f = self.mg.solve(rhs)
return f
def G(self, x, compute_dx=False, iter_data=None, abs_tol=None,
rel_tol=None, mesgs=_MESGS):
"""Objective function."""
g = self.get_f(u=x) - self.f
dx = None
return g, dx
[docs]class DennisSchnabel(object):
"""Test problems from Dennis and Schnabel.
Examples
--------
Check that the exact answers `ans` are solutions:
>>> DennisSchnabel.test_ans() # doctest: +NORMALIZE_WHITESPACE
{'ExtendedPowell': True,
'ExtendedRosenbrock': True,
'HelicalValley': True}
"""
@staticmethod
def test_ans():
"""Test the classes that the exact answer is a solution."""
results = {}
for name, Prob in DennisSchnabel.__dict__.items():
if (name.startswith('_') or
type(Prob) is not type):
# These are not problems.
continue
if (not hasattr(Prob, 'ans')):
# These have no exact solution.
continue
if not Prob.has_G:
continue
if hasattr(Prob, 'm'):
# These have variable length
res = []
for m in xrange(1,5):
p = Prob(m=m)
res.append(np.allclose(p.G(x=p.ans), 0))
else:
# These have fixed length
p = Prob()
res = np.allclose(p.G(x=p.ans)[0], 0)
results[name] = np.all(res)
return results
class ExtendedRosenbrock(solver_interface.Problem):
"""The Extended Rosenbrock Function"""
_state_vars = [
('m', 1, "Integer specifying size of problem `n=2m`."),
('start_factor', 1, "Distance of x0 from origin: 1, 10, or 100"),
('state0', Computed),
('has_G', True),
]
process_vars()
def __init__(self, *v, **kw):
self.state0 = solver_interface.BareState(x_=np.zeros(len(self)))
def __len__(self):
"""Dimension of problem"""
return 2*self.m
def G(self, x, compute_dx=False, iter_data=None, abs_tol=None,
rel_tol=None, mesgs=_MESGS):
"""Objective function"""
i = np.arange(0, self.m)
g = 0*x
g[2*i] = 10*(x[2*i+1] - x[2*i]**2)
g[2*i+1] = 1 - x[2*i]
return g, None
def get_initial_state(self):
x0 = np.ones(len(self), dtype=float)
x0[0::2] *= -1.2
return self.state0.copy_from_x_(x0*self.start_factor)
@property
def ans(self):
"""The solution"""
x_sol = np.ones(len(self), dtype=float)
return x_sol
[docs] class ExtendedPowell(solver_interface.Problem):
"""The Extended Powell Singular Function.
Note that both $G$ and $\nabla^2\norm{G}^2` are singular at
the solution.
"""
_state_vars = [
('m', 1, "Integer specifying size of problem `n=4m`."),
('start_factor', 1, "Distance of x0 from origin: 1, 10, or 100"),
('state0', Computed),
]
process_vars()
def __init__(self, *v, **kw):
self.state0 = solver_interface.BareState(x_=np.zeros(len(self)))
def __len__(self):
"""Dimension of problem"""
return 4*self.m
def G(self, x, compute_dx=False, iter_data=None, abs_tol=None,
rel_tol=None, mesgs=_MESGS):
"""Objective function"""
i = np.arange(0, self.m - 1)
g = 0*x
g[4*i] = x[4*i] + 10*x[4*i+1]
g[4*i+1] = np.sqrt(5)*(x[4*i+2] - x[4*i+3])
g[4*i+2] = (x[4*i+1] - 2*x[4*i+2])**2
g[4*i+3] = np.sqrt(10)*(x[4*i] - x[4*i+3])**2
return g, None
def get_initial_state(self):
x0 = np.ones(len(self), dtype=float)
x0[0::4] = 3
x0[1::4] = -1
x0[2::4] = 0
x0[3::4] = 1
return self.state0.copy_from_x_(x0*self.start_factor)
@property
def ans(self):
"""The solution"""
x_sol = np.zeros(len(self), dtype=float)
return x_sol
class Trigonometric(solver_interface.Problem):
"""A Trigonometric Function."""
_state_vars = [
('m', 1, "Integer specifying size of problem `n=m`."),
('start_factor', 1, "Distance of x0 from origin: 1, 10, or 100"),
('state0', Computed),
]
process_vars()
def __init__(self, *v, **kw):
self.state0 = solver_interface.BareState(x_=np.zeros(len(self)))
def __len__(self):
"""Dimension of problem"""
return self.m
def G(self, x, compute_dx=False, iter_data=None, abs_tol=None,
rel_tol=None, mesgs=_MESGS):
"""Objective function"""
x = np.asarray(x)
xi = x[None, :]
xj = x[:, None]
i = np.arange(0, self.m) + 1
n = self.m + 1
g = n - (np.cos(xj) + i*(1- np.cos(xi)) - np.sin(xi)).sum(axis=0)
return g
def get_initial_state(self):
n = self.m + 1
x0 = np.ones(len(self), dtype=float)/n
return self.state0.copy_from_x_(x0*self.start_factor)
class HelicalValley(solver_interface.Problem):
"""The Helical Valley Function.
"""
_state_vars = [
('start_factor', 1, "Distance of x0 from origin: 1, 10, or 100"),
('state0', Computed),
]
process_vars()
def __init__(self, *v, **kw):
self.state0 = solver_interface.BareState(x_=np.zeros(len(self)))
def __len__(self):
"""Dimension of problem"""
return 3
def _Theta(self, x1, x2):
return (np.arctan2(x2, x1)/2/np.pi
+ np.where(x1 >= 0, 0, 0.5))
def G(self, x, compute_dx=False, iter_data=None, abs_tol=None,
rel_tol=None, mesgs=_MESGS):
"""Objective function"""
g = np.array([10*(x[2] - 10*self._Theta(x[0], x[1])),
10*(np.linalg.norm(x[:2]) -1),
x[2]])
return g, None
def get_initial_state(self):
x0 = np.array([-1.0, 0.0, 0.0])
return self.state0.copy_from_x_(x0*self.start_factor)
@property
def ans(self):
"""The solution"""
x_sol = np.array([1.0, 0.0, 0.0])
return x_sol
class Wood(solver_interface.Problem):
"""The Wood Function. This is an optimization problem. The
objective function here is the gradient.
"""
_state_vars = [
('start_factor', 1, "Distance of x0 from origin: 1, 10, or 100"),
('state0', Computed),
]
process_vars()
def __init__(self, *v, **kw):
self.state0 = solver_interface.BareState(x_=np.zeros(len(self)))
def __len__(self):
"""Dimension of problem"""
return 4
def f(self, x):
"""The objective function"""
return (100*(x[1]**2 - x[1]**2)
+ (1 - x[0])**2
+ 90*(x[2]**2 - x[3])**2
+ (1- x[3])**2
+ 10.1*((1 - x[1])**2 + (1 - x[4])**2)
+ 19.8*(1 - x[1])*(1 - x[3]))
def get_initial_state(self):
x0 = np.array([-3.0, -1.0, -3.0, -1.0])
return self.state0.copy_from_x_(x0*self.start_factor)
@property
def ans(self):
"""The solution"""
x_sol = np.array([1.0, 1.0, 1.0, 1.0])
return x_sol
[docs]class Mine(object):
def test_ans():
"""Test the classes that the exact answer is a solution."""
results = {}
for name, Prob in Mine.__dict__.items():
if (name.startswith('_') or
type(Prob) is not type):
# These are not problems.
continue
if (not hasattr(Prob, 'ans')):
# These have no exact solution.
continue
if hasattr(Prob, 'm'):
# These have variable length
res = []
for m in xrange(1,5):
p = Prob(m=m)
res.append(np.allclose(p.G(x=p.ans), 0))
else:
# These have fixed length
p = Prob()
res = np.allclose(p.G(x=p.ans)[0], 0)
results[name] = np.all(res)
return results
[docs] class Sqrt(solver_interface.Problem):
r"""Computes the square root of two numbers with very
different magnitudes. This is a simple diagonal problem, but
sometimes solvers mix the two solutions causing problems. A
straigh-forward application of Broyden, for example, performs
very badly.
Examples
--------
>>> from mmf.solve.broyden import JinvBroyden
>>> B = JinvBroyden()
>>> p = Mine.Sqrt()
>>> x0 = p.get_initial_state().x_
>>> for n in xrange(100):
... incompatibility = B.update(x=x0, g=p.G(x0)[0])
... x0 = B.get_x()
>>> abs(x0 - p.ans).max() < 10
False
"""
_state_vars = [
('n', np.array([0.01,100.0]), "Computes the square root of this."),
('start_factor', 1, "Distance of x0 from origin: 1, 10, or 100"),
('state0', Computed),
]
process_vars()
def __init__(self, *v, **kw):
self.state0 = solver_interface.BareState(x_=self.n)
def __len__(self):
"""Dimension of problem"""
return 2
def get_initial_state(self):
return self.state0.copy_from_x_(self.n*self.start_factor)
def G(self, x, compute_dx=False, iter_data=None, abs_tol=None,
rel_tol=None, mesgs=_MESGS):
"""Objective function"""
g = x*x - self.n
dx = None
if compute_dx:
dx = x - g/2/x
return g, dx
@property
def ans(self):
"""The solution"""
x_sol = np.sqrt(self.n)
return x_sol