From fc079bfb91aac2084734cd9e52dcb7dd6b3384f4 Mon Sep 17 00:00:00 2001 From: "Taewon D. Kim" Date: Tue, 10 Apr 2018 14:24:36 -0400 Subject: [PATCH] Attempt at new API --- flik/__init__.py | 41 ----- flik/algorithm.py | 11 ++ flik/approx_jacobian.py | 204 --------------------- flik/jacobian.py | 214 ---------------------- flik/model.py | 63 +++++++ flik/nonlinear.py | 294 ------------------------------ flik/subproblem.py | 68 +++++++ flik/test/__init__.py | 25 --- flik/test/test_approx_jacobian.py | 192 ------------------- flik/test/test_gauss_newton.py | 196 -------------------- flik/test/test_jacobian.py | 245 ------------------------- flik/test/test_newton.py | 234 ------------------------ flik/trustregion.py | 13 ++ 13 files changed, 155 insertions(+), 1645 deletions(-) delete mode 100644 flik/__init__.py create mode 100644 flik/algorithm.py delete mode 100644 flik/approx_jacobian.py delete mode 100644 flik/jacobian.py create mode 100644 flik/model.py delete mode 100644 flik/nonlinear.py create mode 100644 flik/subproblem.py delete mode 100644 flik/test/__init__.py delete mode 100644 flik/test/test_approx_jacobian.py delete mode 100644 flik/test/test_gauss_newton.py delete mode 100644 flik/test/test_jacobian.py delete mode 100644 flik/test/test_newton.py create mode 100644 flik/trustregion.py diff --git a/flik/__init__.py b/flik/__init__.py deleted file mode 100644 index 8fbf671..0000000 --- a/flik/__init__.py +++ /dev/null @@ -1,41 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -""" -An experimental local optimization package. - -Copyright (C) 2018 Ayers Lab . - -""" - - -from flik.nonlinear import nonlinear_solve - -from flik.jacobian import Jacobian - -from flik.approx_jacobian import ForwardDiffJacobian -from flik.approx_jacobian import CentralDiffJacobian - - -__all__ = [ - "nonlinear_solve", - "Jacobian", - "ForwardDiffJacobian", - "CentralDiffJacobian", - ] diff --git a/flik/algorithm.py b/flik/algorithm.py new file mode 100644 index 0000000..752c00c --- /dev/null +++ b/flik/algorithm.py @@ -0,0 +1,11 @@ +def algorithm(x, model, subproblem_solver, trustregion, niter): + for i in range(niter): + # solve model (using given algorithm) + step = subproblem_solver(x, model, trustregion) + # update trust region and model + trustregion.update(x, step, model) + model.update(x, step) + # update x + if trustregion.check_step(step): + x += step + return x diff --git a/flik/approx_jacobian.py b/flik/approx_jacobian.py deleted file mode 100644 index fb53317..0000000 --- a/flik/approx_jacobian.py +++ /dev/null @@ -1,204 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -r""" -Classes for numerical approximations for Jacobians of analytical functions. - -Numerical approximations to the Jacobian are useful for optimization purposes -where the analytical Jacobian is unavailable or prohibitively expensive to -compute. - -The forward difference Jacobian approximation uses the formula: -..math:: \frac{\partial f_i(x)}{\partial x_j} - = \frac{f(x + \epsilon e_j) - f(x)}{\epsilon} - -The central difference Jacobian approximation uses the formula: -..math:: \frac{\partial f_i(x)}{\partial x_j} - = \frac{f(x + \epsilon e_j) - f(x - \epsilon e_j)}{2 \epsilon} - -where :math: `e_j` is the unit vector in dimension :math: `j` and :math: -`\epsilon` is a small finite increment over which to approximate the Jacobian. - -""" - - -from numbers import Real, Integral - -import numpy as np - -from flik.jacobian import Jacobian - - -__all__ = [ - "FiniteDiffJacobian", - "ForwardDiffJacobian", - "CentralDiffJacobian", - ] - - -class FiniteDiffJacobian(Jacobian): - r"""Finite difference Jacobian approximation class.""" - - def __init__(self, f, m, n=None, eps=1.0e-4): - r""" - Construct a finite difference approximate Jacobian function. - - Parameters - ---------- - f : callable - The function for which the Jacobian is being approximated. - m : int - Size of the function output vector. - n : int, optional - Size of the function argument vector (default is ``n`` == ``m``). - eps : float or np.ndarray, optional - Increment in the function's argument to use when approximating the - Jacobian. - - Raises - ------ - TypeError - If an argument of an invalid type or shape is passed. - ValueError - If an argument passed has an unreasonable value. - - """ - # Check input types and values - if n is None: - n = m - if not callable(f): - raise TypeError("f must be a callable object") - if not isinstance(m, Integral): - raise TypeError("m must be an integral type") - if not isinstance(n, Integral): - raise TypeError("n must be an integral type") - if m <= 0: - raise ValueError("m must be > 0") - if n <= 0: - raise ValueError("n must be > 0") - if not (isinstance(eps, np.ndarray) and eps.ndim == 1): - if not isinstance(eps, Real): - raise TypeError("eps must be a float or 1-dimensional array") - if isinstance(eps, np.ndarray): - if eps.size != n: - raise ValueError("eps must be of the same length as the input vector") - eps = np.copy(eps) - else: - eps = np.full(int(n), float(eps), dtype=np.float) - if np.any(eps <= 0.0): - raise ValueError("eps must be > 0.0") - # Assign internal attributes - self._function = f - self._m = int(m) - self._n = int(n) - self._eps = eps - - -class ForwardDiffJacobian(FiniteDiffJacobian): - r"""Forward difference Jacobian approximation class.""" - - def __call__(self, x, fx=None): - r""" - Evaluate the approximate Jacobian at position ``x``. - - Parameters - ---------- - x : np.ndarray - Argument vector to the approximate Jacobian function. - fx : np.ndarray, optional - Output vector of the function at position `x` (optional, but avoids - an extra function call). - - Returns - ------- - jacobian : np.ndarray - Value of the approximate Jacobian at position ``x``. - - """ - # Note: In order to stick to row-major iteration, this algorithm - # computes the transpose of the approximate Jacobian into the jac - # vector. This function, being the Jacobian proper, returns the - # transpose of the jac vector. - jac = np.empty((self._n, self._m), dtype=np.float) - # Evaluate function at x (fx = f(x)) if required - if fx is None: - fx = self._function(x) - # Copy x to vector dx - dx = np.copy(x) - # Iterate over elements of `x` to increment - for i in range(self._n): - # Add forward-epsilon increment to dx (dx = x + e_i * eps_i) - dx[i] += self._eps[i] - # Evaluate function at dx (dfx = f(dx)) - dfx = self._function(dx) - # Calculate df[j]/dx[i] = (dfx - fx) / eps_i into dfx vector - dfx -= fx - dfx /= self._eps[i] - # Put result from dfx into the ith row of the jac matrix - jac[i, :] = dfx - # Reset dx = x - dx[i] = x[i] - # df[i]/dx[j] = transpose(jac) - return jac.transpose() - - -class CentralDiffJacobian(FiniteDiffJacobian): - r"""Central difference Jacobian approximation class.""" - - def __call__(self, x): - r""" - Evaluate the approximate Jacobian at position ``x``. - - Parameters - ---------- - x : np.ndarray - Argument vector to the approximate Jacobian function. - - Returns - ------- - jacobian : np.ndarray - Value of the approximate Jacobian at position ``x``. - - """ - # Note: In order to stick to row-major iteration, this algorithm - # computes the transpose of the approximate Jacobian into the jac - # vector. This function, being the Jacobian proper, returns the - # transpose of the jac vector. - jac = np.empty((self._n, self._m), dtype=np.float) - # Copy x to vector dx - dx = np.copy(x) - # Iterate over elements of `x` to increment - for i in range(self._n): - # Add forward-epsilon increment to dx (+dx = x + e_i * eps_i) - dx[i] += self._eps[i] - # Evaluate function at +dx (dfx2 = f(+dx)) - dfx2 = self._function(dx) - # Add backward-epsilon increment to dx (-dx = x - e_i * eps_i) - dx[i] = x[i] - self._eps[i] - # Evaluate function at -dx (dfx1 = f(-dx)) - dfx1 = self._function(dx) - # Calculate df[j]/dx[i] = (dfx2 - dfx1) / (2 * eps_i) into dfx2 vector - dfx2 -= dfx1 - dfx2 /= 2 * self._eps[i] - # Put result from dfx2 into the ith row of the jac matrix - jac[i, :] = dfx2 - # Reset dx = x - dx[i] = x[i] - # df[i]/dx[j] = transpose(jac) - return jac.transpose() diff --git a/flik/jacobian.py b/flik/jacobian.py deleted file mode 100644 index 1f5f07d..0000000 --- a/flik/jacobian.py +++ /dev/null @@ -1,214 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -r"""Base Jacobian class.""" - - -import numpy as np - - -__all__ = [ - "Jacobian", - ] - - -class Jacobian: - r""" - Jacobian class with analytical evaluation by callable Jacobian function. - - The Jacobian class and its subclasses are used for evaluating and updating - Jacobians as part of the Newton iterations. - - """ - - def __init__(self, jac): - r""" - Construct a Jacobian class for a callable analytical jacobian. - - Parameters - ---------- - jac : callable, optional - - Raises - ------ - TypeError - If an argument of an invalid type or shape is passed. - - """ - if not callable(jac): - raise TypeError("J must be a callable object") - self._jac = jac - - def __call__(self, x): - r""" - Compute the Jacobian at position ``x``. - - Parameters - ---------- - x : np.ndarray - - Returns - ------- - y :np.ndarray - - """ - return self._jac(x) - - def update_newton(self, A, new_x, *_): - r""" - Update the Jacobian matrix ``A`` at new solution vector ``x_(k+1)``. - - Parameters - ---------- - new_x : np.ndarray - ``x_(k+1)`` - dx : np.ndarray - ``x_(k+1) - x_k`` - df : np.ndarray - ``f_(k+1) - f_k`` - - """ - A[...] = self(new_x) - - @staticmethod - def update_goodbroyden(A, _, dx, df): - r""" - Update the Jacobian matrix ``A`` at new solution vector ``x_(k+1)``. - - Parameters - ---------- - new_x : np.ndarray - ``x_(k+1)`` - dx : np.ndarray - ``x_(k+1) - x_k`` - df : np.ndarray - ``f_(k+1) - f_k`` - - """ - # Compute Good Broyden right hand side second term numerator - t = df - t -= np.dot(A, dx) - # Divide by dx norm - t /= np.dot(dx, dx) - # Compute matrix from dot product of f and transposed dx - A += np.outer(t, dx.T) - - @staticmethod - def update_badbroyden(A, _, dx, df): - r""" - Update the Jacobian matrix ``A`` at new solution vector ``x_(k+1)``. - - Parameters - ---------- - new_x : np.ndarray - ``x_(k+1)`` - dx : np.ndarray - ``x_(k+1) - x_k`` - df : np.ndarray - ``f_(k+1) - f_k`` - - """ - t2 = np.dot(dx.T, A) - norm = np.dot(t2, df) - t1 = dx - t1 -= np.dot(A, df) - t1 /= norm - A += np.outer(t1, t2) - - @staticmethod - def update_dfp(A, _, dx, df): - r""" - Update the Jacobian matrix ``A`` at new solution vector ``x_(k+1)``. - - Parameters - ---------- - new_x : np.ndarray - ``x_(k+1)`` - dx : np.ndarray - ``x_(k+1) - x_k`` - df : np.ndarray - ``f_(k+1) - f_k`` - - """ - norm = np.dot(df, dx) - t1 = np.outer(df, dx.T) - t1 /= -norm - t1 += np.eye(t1.shape[0]) - t2 = np.outer(dx, df.T) - t2 /= -norm - t2 += np.eye(t2.shape[0]) - A[:] = np.dot(t1, A) - A[:] = np.dot(A, t2) - t1 = np.outer(df, df.T) - t1 /= norm - A += t1 - - @staticmethod - def update_bfgs(A, _, dx, df): - r""" - Update the Jacobian matrix ``A`` at new solution vector ``x_(k+1)``. - - Parameters - ---------- - new_x : np.ndarray - ``x_(k+1)`` - dx : np.ndarray - ``x_(k+1) - x_k`` - df : np.ndarray - ``f_(k+1) - f_k`` - - """ - Jacobian.update_dfp(A, None, df, dx) - - @staticmethod - def update_sr1(A, _, dx, df): - r""" - Update the Jacobian matrix ``A`` at new solution vector ``x_(k+1)``. - - Parameters - ---------- - new_x : np.ndarray - ``x_(k+1)`` - dx : np.ndarray - ``x_(k+1) - x_k`` - df : np.ndarray - ``f_(k+1) - f_k`` - - """ - t1 = df - np.dot(A, dx) - t2 = np.outer(t1, t1.T) - t2 /= np.dot(t1.T, dx) - A += t2 - - @staticmethod - def update_sr1inv(A, _, dx, df): - r""" - Update the Jacobian matrix ``A`` at new solution vector ``x_(k+1)``. - - Parameters - ---------- - new_x : np.ndarray - ``x_(k+1)`` - dx : np.ndarray - ``x_(k+1) - x_k`` - df : np.ndarray - ``f_(k+1) - f_k`` - - """ - Jacobian.update_sr1(A, None, df, dx) diff --git a/flik/model.py b/flik/model.py new file mode 100644 index 0000000..e1d5cec --- /dev/null +++ b/flik/model.py @@ -0,0 +1,63 @@ +import numpy as np + + +class Model: + def __init__(self, func): + self.func = func + # or whatever else structure for storing data + self.cache = {} + + def update(self): + # update cache + pass + + +class LinearModel(Model): + def __init__(self, func, grad): + super().__init__(func) + self.grad = grad + + +class QuadraticModel(Model): + def __init__(self, func, grad, hess=None): + super().__init__(func) + self.grad = grad + self.hess = hess + + +class HessianUpdateModel(QuadraticModel): + def __init__(self, func, grad, hess, initial_hessian): + super().__init__(func, grad, hess) + self._hessian = initial_hessian + + def hessian_dot(self, vec): + return self._hessian.dot(vec) + + def hessian_update(self, newx, newf): + # update self.hessian + pass + + +class InverseHessianUpdateModel(QuadraticModel): + def __init__(self, func, grad, hess, initial_inv_hessian): + super().__init__(func, grad, hess) + self._inv_hessian = initial_inv_hessian + + def inv_hessian_dot(self, vec): + return self._inv_hessian.dot(vec) + + def inv_hessian_update(self, newx, newf): + # update self.inv_hessian + pass + + +class ExactHessianModel(QuadraticModel): + def __init__(self, func, grad, hess, initial_hessian): + super().__init__(func, grad, hess) + + def hessian_dot(self, vec): + return self.hess(vec).dot(vec) + + def inv_hessian_dot(self, vec): + hessian = self.hess(vec) + return np.linalg.vec(hessian).dot(vec) diff --git a/flik/nonlinear.py b/flik/nonlinear.py deleted file mode 100644 index d7d9b88..0000000 --- a/flik/nonlinear.py +++ /dev/null @@ -1,294 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -r""" -Solvers for nonlinear systems using the Newton and Gauss-Newton algorithms. - -Functions used to find the roots of a nonlinear system of equations given the -residual function `f` and its analytical Jacobian `J`. An exactly-determined -system (`m` equations, `m` variables) is best solved with the Newton method. -Over- or under- determined systems (`m` equations, `n` variables) must be -solved in the least-squares sense using the Gauss-Newton method. - -""" - - -from numbers import Integral -from numbers import Real - -import numpy as np - -from flik.jacobian import Jacobian -from flik.approx_jacobian import CentralDiffJacobian - - -__all__ = [ - "nonlinear_solve", - ] - - -def nonlinear_solve(f, x_0, J=None, stepsize=1.0, eps=1.0e-6, maxiter=100, method="newton"): - r""" - Solve a system of nonlinear equations with the Newton method. - - Parameters - ---------- - f : callable - Vector-valued function corresponding to nonlinear system of equations. - Must be of the form f(x), where x is a 1-dimensional array. - x_0 : np.ndarray - Solution initial guess. - J : callable, optional - Jacobian of function f. Must be of the form J(x), where x is a - 1-dimensional array. If none is given, then the Jacobian is calculated - using finite differences. - stepsize : float or np.ndarray, optional - Scaling factor for Newton step. - eps : float, optional - Convergence threshold for vector function f norm. - maxiter : int, optional - Maximum number of iterations to perform. - method : str, optional - Update method for the (approximated) J(x) or the inverse of J(x). The - default uses Newton method. - - Returns - ------- - result : dict - A dictionary with the keys: - success - Boolean variable informing whether the algorithm succeeded or not. - message - Information about the cause of the termination. - niter - Number of actual iterations performed. - x - Nonlinear system of equations solution (Root). - f - Vector function evaluated at solution. - J - Jacobian evaluated at solution. - eps - Convergence threshold for vector function f norm. - - Raises - ------ - TypeError - If an argument of an invalid type or shape is passed. - ValueError - If an argument passed has an unreasonable value. - - """ - # Check input types and values - if not callable(f): - raise TypeError("Argument f should be callable") - if not callable(J) and J is not None: - raise TypeError("Argument J should be callable") - if not (isinstance(x_0, np.ndarray) and x_0.ndim == 1): - raise TypeError("Argument x_0 should be a 1-dimensional numpy array") - if not isinstance(eps, Real): - raise TypeError("Argument eps should be a real number") - if not isinstance(maxiter, Integral): - raise TypeError("Argument maxiter should be an integer number") - if not isinstance(method, str): - raise TypeError("Argument method should be a string") - if eps < 0.0: - raise ValueError("Argument eps should be >= 0.0") - if maxiter < 1: - raise ValueError("Argument maxiter should be >= 1") - eps = float(eps) - maxiter = int(maxiter) - # Check stepsize argument - if isinstance(stepsize, Real): - stepsize = float(stepsize) - elif isinstance(stepsize, np.ndarray): - if stepsize.shape != x_0.shape: - raise TypeError("stepsize and x_0 must have the same shape") - else: - raise TypeError("Argument stepsize should be a float or numpy array") - if np.any(stepsize < 0.0): - raise ValueError("Argument stepsize should be >= 0.0") - # Check J (Jacobian function) argument - if J is None: - m = f(x_0).shape[0] - n = x_0.shape[0] - J = CentralDiffJacobian(f, m, n) - else: - J = J if isinstance(J, Jacobian) else Jacobian(J) - # Choose the step/update function and inverse option - inverse, step, update = _nonlinear_functions(J, method) - # Return result of Newton iterations - return _nonlinear_iterations(f, x_0, J, stepsize, eps, maxiter, inverse, step, update) - - -def _nonlinear_functions(J, method): - r""" - Return the functions used in ``nonlinear_solve`` according to the method. - - Parameters - ---------- - J : Jacobian - method : str - - Returns - ------- - inv : bool - True if inverse Jacobian is used, otherwise False - step : function - Newton step function - update : function - Jacobian update function - - """ - method = method.lower() - if method == "newton": - result = False, _step_linear, J.update_newton - elif method == "goodbroyden": - result = False, _step_linear, J.update_goodbroyden - elif method == "dfp": - result = False, _step_linear, J.update_dfp - elif method == "sr1": - result = False, _step_linear, J.update_sr1 - elif method == "badbroyden": - result = True, _step_inverse, J.update_badbroyden - elif method == "bfgs": - result = True, _step_inverse, J.update_bfgs - elif method == "sr1inv": - result = True, _step_inverse, J.update_sr1inv - elif method == "gaussnewton": - result = False, _step_gauss_newton, J.update_newton - else: - raise ValueError("Argument method is not a valid option") - return result - - -def _nonlinear_iterations(f, x_0, J, stepsize, eps, maxiter, inverse, step, update): - r"""Run the iterations for ``newton_solve``.""" - # Calculate f_0 and J_0 - A = np.linalg.inv(J(x_0)) if inverse else J(x_0) - b = f(x_0) - # Iterations - success = False - message = "Maximum number of iterations reached." - for niter in range(1, maxiter + 1): - b *= -1 - # Calculate step function, take Newton step - try: - dx = step(b, A) - dx *= stepsize - except np.linalg.LinAlgError: - message = "Singular Jacobian; no solution found." - break - x_0 += dx - # Evaluate function and Jacobian for next step or result - df = b - b = f(x_0) - df += b - update(A, x_0, dx, df) - # Check for convergence - if np.linalg.norm(b) < eps: - # If so, we're done (SUCCESS) - success = True - message = "Convergence obtained." - break - # Return result dictionary - return { - "success": success, - "message": message, - "niter": niter, - "x": x_0, - "f": b, - "J": np.linalg.inv(A) if inverse else A, - "eps": eps, - } - - -def _step_linear(b, A): - r""" - Compute the Newton step for the Jacobian. - - Calculate the roots for the next step of a method that updates the - (approximated) Jacobian matrix. - - Parameters - ---------- - b : np.ndarray - 1-dimensional array of the negative of the function evaluated at the - current guess of the roots -f(x_0). - A : np.ndarray - 2-dimensional array of the Jacobian evaluated at the current guess of - the roots J(x_0). - - Returns - ------- - dx : np.ndarray - Step length for the method. - - """ - return np.linalg.solve(A, b) - - -def _step_inverse(b, A): - r""" - Compute the Newton step for the inverse of the Jacobian. - - Calculate the roots for the next step of a method that updates the - inverse of the approximated Jacobian matrix. - - Parameters - ---------- - b : np.ndarray - 1-dimensional array of the negative of the function evaluated at the - current guess of the roots -f(x_0). - A : np.ndarray - 2-dimensional array of the Jacobian evaluated at the current guess of - the roots J(x_0). - - Returns - ------- - dx : np.ndarray - Step length for the method. - - """ - return np.dot(A, b) - - -def _step_gauss_newton(b, A): - r""" - Compute the Gauss-Newton step for the Jacobian. - - Calculate the roots for the next step of a method that updates the - (approximated) Jacobian matrix. - - Parameters - ---------- - b : np.ndarray - 1-dimensional array of the negative of the function evaluated at the - current guess of the roots -f(x_0). - A : np.ndarray - 2-dimensional array of the Jacobian evaluated at the current guess of - the roots J(x_0). - - Returns - ------- - dx : np.ndarray - Step length for the method. - - """ - return np.linalg.solve(np.dot(A.T, A), np.dot(A.T, b)) diff --git a/flik/subproblem.py b/flik/subproblem.py new file mode 100644 index 0000000..d290a29 --- /dev/null +++ b/flik/subproblem.py @@ -0,0 +1,68 @@ +"""Problem of finding the minima in the given trust region.""" +from flik.model import HessianUpdateModel, InverseHessianUpdateModel, ExactHessianModel + + +class HessianDotSubproblem: + def __init__(self, model): + if not isinstance(model, (HessianUpdateModel, ExactHessianModel)): + raise ValueError + self.model = model + + def solve(self, x, trustregion): + raise NotImplementedError + + +class InvHessianDotSubproblem: + def __init__(self, model): + if not isinstance(model, (InverseHessianUpdateModel, ExactHessianModel)): + raise ValueError + self.model = model + + def solve(self, x, trustregion): + raise NotImplementedError + + +class Dogleg(InvHessianDotSubproblem): + pass + + +class Subspace(InvHessianDotSubproblem): + pass + + +class CauchyPoint(InvHessianDotSubproblem): + pass + + +class Steinhaug(HessianDotSubproblem): + pass + + +class Iterative(HessianDotSubproblem): + pass + + +# OR ALTERNATIVELY +def solve_dogleg(x, model, trustregion): + if not isinstance(model, (InverseHessianUpdateModel, ExactHessianModel)): + raise ValueError + + +def solve_subspace(x, model, trustregion): + if not isinstance(model, (InverseHessianUpdateModel, ExactHessianModel)): + raise ValueError + + +def solve_cauchy_point(x, model, trustregion): + if not isinstance(model, (InverseHessianUpdateModel, ExactHessianModel)): + raise ValueError + + +def solve_steinhaug(x, model, trustregion): + if not isinstance(model, (HessianUpdateModel, ExactHessianModel)): + raise ValueError + + +def solve_iterative(x, model, trustregion): + if not isinstance(model, (HessianUpdateModel, ExactHessianModel)): + raise ValueError diff --git a/flik/test/__init__.py b/flik/test/__init__.py deleted file mode 100644 index e5f695d..0000000 --- a/flik/test/__init__.py +++ /dev/null @@ -1,25 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -""" -An experimental local optimization package. - -Copyright (C) 2018 Ayers Lab . - -""" diff --git a/flik/test/test_approx_jacobian.py b/flik/test/test_approx_jacobian.py deleted file mode 100644 index d07a19c..0000000 --- a/flik/test/test_approx_jacobian.py +++ /dev/null @@ -1,192 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -"""Test file for `flik.approx_jacobian`.""" - - -import numpy as np -import numpy.testing as npt - -from flik import ForwardDiffJacobian -from flik import CentralDiffJacobian - - -__all__ = [ - "test_finite_diff_jacobian_inputs", - "test_forward_diff_jacobian_square", - "test_central_diff_jacobian_square", - "test_forward_diff_jacobian_rectangular", - "test_central_diff_jacobian_rectangular", - ] - - -# Seed the numpy rng for consistency -np.random.seed(101010101) - - -# Define some analytical test functions and Jacobians - - -def f1(x): - """Test function.""" - y = np.copy(x) - y **= 2 - y[0] += x[1] - y[1] -= x[0] - return y - - -def j1(x): - """Test function.""" - y = np.empty((2, 2), dtype=x.dtype) - y[0, 0] = 2.0 * x[0] - y[0, 1] = 1.0 - y[1, 0] = -1.0 - y[1, 1] = 2.0 * x[1] - return y - - -def f2(x): - """Test function.""" - y = np.empty((3,), dtype=x.dtype) - y[0] = x[0] * x[1] - y[1] = x[0] ** 3 - np.sqrt(x[1]) - y[2] = 4.0 * x[0] * x[1] ** 2 - return y - - -def j2(x): - """Test function.""" - y = np.empty((3, 2), dtype=x.dtype) - y[0, 0] = x[1] - y[0, 1] = x[0] - y[1, 0] = 3.0 * x[0] ** 2 - y[1, 1] = -0.5 * x[1] ** (-0.5) - y[2, 0] = 4.0 * x[1] ** 2 - y[2, 1] = 8.0 * x[0] * x[1] - return y - - -# Run tests - - -def test_finite_diff_jacobian_inputs(): - """Test invalid inputs to FiniteDiffJacobian.""" - # Test `f` argument - npt.assert_raises(TypeError, ForwardDiffJacobian, "string", 2) - # Test `m` argument - npt.assert_raises(TypeError, ForwardDiffJacobian, f1, 3.14159) - # Test `n` argument - npt.assert_raises(TypeError, ForwardDiffJacobian, f1, 2, 3.14159) - # Test `eps` argument - npt.assert_raises(TypeError, ForwardDiffJacobian, f1, 2, eps="string") - # Test negative `n` - npt.assert_raises(ValueError, ForwardDiffJacobian, f1, -2) - # Test negative `m` - npt.assert_raises(ValueError, ForwardDiffJacobian, f1, 2, -2) - # Test `eps` array size - npt.assert_raises(ValueError, ForwardDiffJacobian, f1, 2, eps=np.full((5,), 0.0078125)) - # Test negative `eps` scalar - npt.assert_raises(ValueError, ForwardDiffJacobian, f1, 2, eps=-0.24681357) - # Test negative `eps` array - npt.assert_raises(ValueError, ForwardDiffJacobian, f1, 2, eps=np.full((2,), -0.0078125)) - - -def test_forward_diff_jacobian_square(): - """Test ForwardDiffJacobian against square system `f1` and `J1`.""" - j = ForwardDiffJacobian(f1, 2) - x = np.random.rand(2) - jx = j(x) - j1x = j1(x) - diff = np.abs(jx - j1x) - assert np.all(diff < 1.0e-3) - - j = ForwardDiffJacobian(f1, 2) - x = np.random.rand(2) - f1x = f1(x) - jx = j(x, f1x) - j1x = j1(x) - diff = np.abs(jx - j1x) - assert np.all(diff < 1.0e-3) - - j = ForwardDiffJacobian(f1, 2, eps=np.full((2,), 1.0e-3)) - x = np.random.rand(2) - jx = j(x) - j1x = j1(x) - diff = np.abs(jx - j1x) - assert np.all(diff < 5.0e-3) - - -def test_central_diff_jacobian_square(): - """Test CentralDiffJacobian against square system `f1` and `J1`.""" - j = CentralDiffJacobian(f1, 2) - x = np.random.rand(2) - jx = j(x) - j1x = j1(x) - diff = np.abs(jx - j1x) - assert np.all(diff < 1.0e-3) - - j = CentralDiffJacobian(f1, 2, eps=np.full((2,), 1.0e-3)) - x = np.random.rand(2) - jx = j(x) - j1x = j1(x) - diff = np.abs(jx - j1x) - assert np.all(diff < 5.0e-3) - - -def test_forward_diff_jacobian_rectangular(): - """Test ForwardDiffJacobian against rectangular system `f2` and `J2`.""" - j = ForwardDiffJacobian(f2, 3, 2) - x = np.random.rand(2) - jx = j(x) - j2x = j2(x) - diff = np.abs(jx - j2x) - assert np.all(diff < 1.0e-3) - - j = ForwardDiffJacobian(f2, 3, 2) - x = np.random.rand(2) - f2x = f2(x) - jx = j(x, f2x) - j2x = j2(x) - diff = np.abs(jx - j2x) - assert np.all(diff < 1.0e-3) - - j = ForwardDiffJacobian(f2, 3, 2, np.full((2,), 1.0e-4)) - x = np.random.rand(2) - jx = j(x) - j2x = j2(x) - diff = np.abs(jx - j2x) - assert np.all(diff < 5.0e-3) - - -def test_central_diff_jacobian_rectangular(): - """Test CentralDiffJacobian against rectangular system `f2` and `J2`.""" - j = CentralDiffJacobian(f2, 3, 2) - x = np.random.rand(2) - jx = j(x) - j2x = j2(x) - diff = np.abs(jx - j2x) - assert np.all(diff < 1.0e-3) - - j = CentralDiffJacobian(f2, 3, 2, np.full((2,), 1.0e-3)) - x = np.random.rand(2) - jx = j(x) - j2x = j2(x) - diff = np.abs(jx - j2x) - assert np.all(diff < 5.0e-3) diff --git a/flik/test/test_gauss_newton.py b/flik/test/test_gauss_newton.py deleted file mode 100644 index d90b2f2..0000000 --- a/flik/test/test_gauss_newton.py +++ /dev/null @@ -1,196 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -"""Test file for `flik.nonlinear.nonlinear_solve`.""" - - -import numpy as np - -from nose.tools import assert_raises - -from flik import nonlinear_solve - - -__all__ = [ - "test_gauss_newton_linear_solve", - "test_gauss_newton_nonlinear_solve", - "test_gauss_newton_nonlinear_overdetermined_solve", - "test_gauss_newton_nonlinear_overdetermined_solve", - "test_gauss_newton_approximation_overdetermined", - ] - - -def f1(x): - """Test function.""" - return np.array([136 - x[0] - x[1] ** 2. - x[2] ** 3., - 1038 - x[0] - 4. * x[1] ** 2. - 8. * x[2] ** 3., - 3458 - x[0] - 9. * x[1] ** 2. - 27. * x[2] ** 3.]) - - -def j1(x): - """Test function.""" - return np.array([[-1., -2. * x[1], -3. * x[2] ** 2.], - [-1., -8. * x[1], -24. * x[2]**2.], - [-1., -18. * x[1], -81. * x[2]**2.]]) - - -def f_lin(x): - """Test function.""" - return np.array([644. * x[0] + 52. * x[1] - 227, 52. * x[0] + 5. * x[1] - 17]) - - -def j_lin(_): - """Test function.""" - return np.array([[644., 52.], [52., 5.]]) - - -def f3(x): - """Test function.""" - return np.array([136 - x[0] - x[1]**2. - x[2]**3., - 1038 - x[0] - 4. * x[1]**2. - 8. * x[2]**3., - 3458 - x[0] - 9. * x[1]**2. - 27. * x[2]**3., - 8146 - x[0] - 16. * x[1]**2. - 64 * x[2]**3.]) - - -def j3(x): - """Test function.""" - return np.array([[-1., -2. * x[1], -3. * x[2]**2.], - [-1., -8. * x[1], -24. * x[2]**2.], - [-1., -18. * x[1], -81. * x[2]**2.], - [-1., -32. * x[1], -192 * x[2]**2.]]) - - -def f4(x): - """Test function.""" - return np.array([4232. - x[0] - x[1]**2. - x[2]**3. - x[3]**4., - 66574. - x[0] - 4. * x[1]**2. - 8. * x[2]**3. - 16. * x[3]**4., - 335234. - x[0] - 9. * x[1]**2. - 27. * x[2]**3. - 81. * x[3]**4.]) - - -def j4(x): - """Test function.""" - return np.array([[-1, -2. * x[1], -3. * x[2]**2., -4. * x[3]**3.], - [-1., -8. * x[1], -24. * x[2]**2., -64 * x[3]**3.], - [-1., -18. * x[1], -81. * x[2]**2., -324 * x[3]**4.]]) - - -def test_gauss_newton_linear_solve(): - """Test that gauss_newton solves linear systems in 1 step.""" - # Obtained from pg 489 of Numerical Mathematics and Computing Sixth Edition. - x0 = np.array([100., -200.]) - result = nonlinear_solve(f_lin, x0, j_lin, eps=1e-20, method="gaussnewton") - x_expt = np.array([0.4864, -1.6589]) - f_expt = np.array([0., 0.]) - jac_expt = j_lin(1.) - message_expt = "Convergence obtained." - assert np.allclose(f_expt, result['f'], rtol=1e-5, atol=1e-5) - assert np.allclose(jac_expt, result['J'], rtol=1e-5, atol=1e-5) - assert np.allclose(x_expt, result['x'], rtol=1e-4, atol=1e-4) - assert result['success'] - assert result['message'] == message_expt - - -def test_gauss_newton_nonlinear_solve(): - """Test that gauss_newton solves nonlinear systems.""" - # The function being optimized is c0 + c1^2 * x^2 + c2^3 * x^3 - # The points are evaluated on [1., 2., 3.] - x0 = np.array([20., 2.5, 500.]) - result = nonlinear_solve(f1, x0, j1, maxiter=5000, method="gaussnewton") - f_expt = np.array([0.] * 3) - x_expt = np.array([2., 3., 5.]) - message_expt = "Convergence obtained." - jac_expt = j1(x_expt) - assert np.allclose(f_expt, result['f'], rtol=1e-5, atol=1e-5) - assert np.allclose(jac_expt, result['J'], rtol=1e-5, atol=1e-5) - assert np.allclose(x_expt, result['x']) - assert result['success'] - assert result['message'] == message_expt - # Because of symmetry the second coefficient of -3 should work - x0 = np.array([200., -200., 200.]) - result = nonlinear_solve(f1, x0, j1, maxiter=5000, method="gaussnewton") - x_expt = np.array([2., -3., 5.]) - jac_expt = j1(x_expt) - assert np.allclose(f_expt, result['f'], rtol=1e-5, atol=1e-5) - assert np.allclose(jac_expt, result['J'], rtol=1e-5, atol=1e-5) - assert np.allclose(x_expt, result['x']) - assert result['success'] - assert result['message'] == message_expt - - -def test_gauss_newton_nonlinear_overdetermined_solve(): - """Test that gauss_newton solves nonlinear, overdetermined systems.""" - # The function being optimized is c0 + c1^2 * x^2 + c2^3 * x^3 - x0 = np.array([20., 2.5, 500.]) - result = nonlinear_solve(f3, x0, j3, method="gaussnewton") - f_expt = np.array([0.] * 4) - x_expt = np.array([2., 3., 5.]) - message_expt = "Convergence obtained." - jac_expt = j3(x_expt) - assert np.allclose(f_expt, result['f'], rtol=1e-5, atol=1e-5) - assert np.allclose(jac_expt, result['J'], rtol=1e-5, atol=1e-5) - assert np.allclose(x_expt, result['x']) - assert result['success'] - assert result['message'] == message_expt - # Because of symmetry the second coefficient of -3 should work - x0 = np.array([1000., -200., 500.]) - result = nonlinear_solve(f3, x0, j3, method="gaussnewton") - x_expt = np.array([2., -3., 5.]) - message_expt = "Convergence obtained." - jac_expt = j3(x_expt) - assert np.allclose(f_expt, result['f'], rtol=1e-5, atol=1e-5) - assert np.allclose(jac_expt, result['J'], rtol=1e-5, atol=1e-5) - assert np.allclose(x_expt, result['x']) - assert result['success'] - assert result['message'] == message_expt - - -def test_gauss_newton_nonlinear_underdetermined_solve(): - """Test that gauss_newton solves nonlinear, underdetermined systems.""" - # Non-Exact Initial Guess. Here a really Good Initial guess is needed - x_expt = np.array([2., 3., 5., 8.]) - f_expt = np.array([0.] * 3) - x0 = np.array([2.001, 3.001, 5.001, 10.01]) - result = nonlinear_solve(f4, x0, j4, maxiter=12, method="gaussnewton") - message_expt = "Maximum number of iterations reached." - assert result['message'] == message_expt - # Repeat with a higher number of iterations and a really good initial guess - # Still it only converges to the first decimal place. - x0 = np.array([2.00001, 3.00001, 5.00001, 8.00001]) - result = nonlinear_solve(f4, x0, j4, maxiter=10000, method="gaussnewton") - jac_expt = j4(result["x"]) - assert np.allclose(f_expt, result['f'], rtol=1e-5, atol=1e-5) - assert np.allclose(jac_expt, result['J'], rtol=1e-5, atol=1e-5) - assert np.allclose(x_expt, result['x'], rtol=1e-1, atol=1e-1) - assert result['message'] == "Convergence obtained." - assert result['success'] - - -def test_gauss_newton_approximation_overdetermined(): - """Test gauss newton using an approximation finite diff Jacobian.""" - x0 = np.array([20., 2.5, 500.]) - result = nonlinear_solve(f3, x0, method="gaussnewton") - f_expt = np.array([0.] * 4) - x_expt = np.array([2., 3., 5.]) - message_expt = "Convergence obtained." - jac_expt = j3(x_expt) - assert np.allclose(f_expt, result['f'], rtol=1e-5, atol=1e-5) - assert np.allclose(jac_expt, result['J'], rtol=1e-5, atol=1e-5) - assert np.allclose(x_expt, result['x']) - assert result['success'] - assert result['message'] == message_expt diff --git a/flik/test/test_jacobian.py b/flik/test/test_jacobian.py deleted file mode 100644 index 883db84..0000000 --- a/flik/test/test_jacobian.py +++ /dev/null @@ -1,245 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -"""Test file for `flik.jacobian`.""" - - -import numpy as np -import numpy.testing as npt - -from flik import Jacobian - - -__all__ = [ - "test_jacobian_inputs", - "test_update_badbroyden_j1", - "test_update_goodbroyden_j1", - "test_update_dfp", - "test_secant_condition_dfp", - "test_positive_definiteness_dfp", - "test_update_inv_bfgs_j1", - "test_update_inv_bfgs_j3", - "test_bfgs_secant_condition", - ] - - -# Seed the numpy rng for consistency -np.random.seed(101010101) - - -# Define some analytical test functions and Jacobians - - -def f1(x): - """Test function.""" - y = np.copy(x) - y **= 2 - y[0] += x[1] - y[1] -= x[0] - return y - - -def j1(x): - """Test function.""" - y = np.empty((2, 2), dtype=x.dtype) - y[0, 0] = 2.0 * x[0] - y[0, 1] = 1.0 - y[1, 0] = -1.0 - y[1, 1] = 2.0 * x[1] - return y - - -def f2(x): - """Test function.""" - return np.array([x[0]**2 + x[1]**2 - 1, 3. * x[0]**2. + x[0] * x[1]**2. + x[2], - x[2]**2. + x[1]]) - - -def j2(x): - """Test function.""" - return np.array([[2. * x[0], 2. * x[1], 0.], - [6. * x[0] + x[1]**2., 2. * x[0] * x[1], 1.], - [0., 1., 2. * x[2]]]) - - -def test_jacobian_inputs(): - """Test invalid inputs to Jacobian class.""" - # Test adding a non-callable jacobian to the class. - npt.assert_raises(TypeError, Jacobian, jac=5.) - # Test Jacobian calls the correct function. - jac_obj = Jacobian(jac=j1) - x_0 = np.random.rand(2) - assert np.allclose(jac_obj(x_0), j1(x_0)) - jac_obj_obj = Jacobian(jac_obj) - assert np.allclose(jac_obj_obj(x_0), j1(x_0)) - - -def test_update_goodbroyden_j1(): - """Test Jacobian.update_goodbroyden using analytical Jacobian and function.""" - # Initial point x_k - x_0 = np.random.rand(2) - # Function at initial x_k - f1_0 = f1(x_0) - # Jacobian at initial x_k - j1_0 = j1(x_0) - # Define arbitrary step - dx = np.asarray([0.01, 0.01]) - # Take step dx - x_1 = x_0 + dx - # Function at x_{k+1} - f1_1 = f1(x_1) - # Compute analytical Jacobian at x_{k+1}` - j1_1 = j1(x_1) - delta_f = f1_1 - f1_0 - expected_ans = j1_0 + \ - np.outer((delta_f - j1_0.dot(dx)) / np.dot(dx, dx), dx) - # Approximate Jacobian at (x_0 + dx) with Good Broyden method - # from Jacobian and vector function at initial x_0 - Jacobian.update_goodbroyden(j1_0, x_1, x_1 - x_0, f1_1 - f1_0) - assert np.allclose(j1_0, expected_ans) - assert np.allclose(j1_1, expected_ans, atol=1e-1) - - -def test_update_badbroyden_j1(): - """Test Jacobian.update_badbroyden against analytical Jacobian and function.""" - # Initial point x_0 - x_0 = np.random.rand(2) - # Function at initial x_k - f1_0 = f1(x_0) - # Inverse Jacobian at initial x_0 - j1_inv_0 = np.linalg.inv(j1(x_0)) - # Define arbitrary step and take step - dx = np.asarray([0.01, 0.01]) - x_1 = x_0 + dx - # Function at x_{k+1} - f1_1 = f1(x_1) - delta_f = f1_1 - f1_0 - # Compute analytical inverse Jacobian at `x_0 + dx` - expected_ans = j1_inv_0 + np.outer((dx - j1_inv_0.dot(delta_f)) / - np.dot(delta_f, delta_f), delta_f) - # Approximate inverse Jacobian at (x_0 + dx) with Bad Broyden method - # from inverse Jacobian and vector function at initial x_0 - Jacobian.update_badbroyden(j1_inv_0, x_1, dx, delta_f) - assert np.allclose(j1_inv_0, expected_ans, atol=1e-1) - # assert that the approximation jacobian matches the analytic jacobian. - j1_1 = np.linalg.inv(j1(x_1)) - assert np.allclose(j1_1, expected_ans, atol=1e-1) - - -def test_update_dfp(): - """Test the definition of DFP.""" - # Test wikipedia definition/format. - random_vecs = np.random.rand(10, 3) * 50. + 0.5 - x = random_vecs[0] - dx = np.array([0.0001, 0.0001, 0.0001]) - df = f2(x + dx) - f2(x) - b0 = j2(x) - gamma = 1. / np.dot(df, dx) - expected_answer = (np.eye(3) - - gamma * np.outer(df, dx)).dot(b0).dot(np.eye(3) - - gamma * np.outer(dx, df)) - expected_answer += gamma * np.outer(df, df) - Jacobian.update_dfp(b0, x + dx, dx, df) - assert np.allclose(b0, expected_answer) - - -def test_positive_definiteness_dfp(): - """Test that DFP is positive definite.""" - # It's crucial that the x values are greater than 0.5 - # or else positive definiteness is not required. - random_vecs = np.random.rand(10, 3) * 50. + 0.5 - dx = np.array([0.5, 0.5, 0.5]) - x = random_vecs[0] - dfp = np.eye(3) - Jacobian.update_dfp(dfp, x + dx, dx, f2(x + dx) - f2(x)) - assert np.all(np.asarray([np.dot(x, dfp.dot(x)) for x in random_vecs]) > 0.) - - -def test_secant_condition_dfp(): - """Test that dfp update satisfies secant condition.""" - x = (np.random.rand(1, 3) * 50. - 10.)[0] - dx = np.array([0.5, 0.5, 0.5]) - df = f2(x + dx) - f2(x) - dfp = np.eye(3) - Jacobian.update_dfp(dfp, x + dx, dx, f2(x + dx) - f2(x)) - assert np.allclose(df, dfp.dot(dx)) - - -def test_update_sr1_j1(): - """Test Jacobian.update_sr1 using analytical Jacobian and function.""" - # Initial point x_k - x_0 = np.random.rand(2) - # Function at initial x_k - f1_0 = f1(x_0) - # Jacobian at initial x_k - j1_0 = j1(x_0) - # Define arbitrary step - dx = np.asarray([0.01, 0.01]) - # Take step dx - x_1 = x_0 + dx - # Function at x_{k+1} - f1_1 = f1(x_1) - # Compute analytical Jacobian at x_{k+1}` - j1_1 = j1(x_1) - delta_f = f1_1 - f1_0 - tmp = delta_f - j1_0.dot(dx) - expected_ans = j1_0 + (np.outer(tmp, tmp.T)) / np.dot(tmp.T, dx) - # Approximate Jacobian at (x_0 + dx) with Good Broyden method - # from Jacobian and vector function at initial x_0 - Jacobian.update_sr1(j1_0, x_1, x_1 - x_0, f1_1 - f1_0) - assert np.allclose(j1_0, expected_ans) - assert np.allclose(j1_1, expected_ans, atol=1e-1) - - -def test_update_sr1_inv_j1(): - """Test Jacobian.update_sr1_inv against analytical Jacobian and function.""" - # Initial point x_0 - x_0 = np.random.rand(2) - # Function at initial x_k - f1_0 = f1(x_0) - # Inverse Jacobian at initial x_0 - j1_inv_0 = np.linalg.inv(j1(x_0)) - # Define arbitrary step and take step - dx = np.asarray([0.01, 0.01]) - x_1 = x_0 + dx - # Function at x_{k+1} - f1_1 = f1(x_1) - delta_f = f1_1 - f1_0 - # Compute analytical inverse Jacobian at `x_0 + dx` - tmp = dx - j1_inv_0.dot(delta_f) - expected_ans = j1_inv_0 + (np.outer(tmp, tmp.T)) / np.dot(tmp.T, dx) - # Approximate inverse Jacobian at (x_0 + dx) with Bad Broyden method - # from inverse Jacobian and vector function at initial x_0 - Jacobian.update_sr1inv(j1_inv_0, x_1, x_1 - x_0, f1_1 - f1_0) - assert np.allclose(j1_inv_0, expected_ans, atol=1e-1) - # assert that the approximation jacobian matches the analytic jacobian. - j1_1 = np.linalg.inv(j1(x_1)) - assert np.allclose(j1_1, expected_ans, atol=1e-1) - - -def test_bfgs_secant_condition(): - """Test the secant condition for BFGS.""" - x = (np.random.rand(1, 3) * 50. - 10.)[0] - dx = np.array([0.5, 0.5, 0.5]) - df = f2(x + dx) - f2(x) - # Compute the updated bfgs Jacobian - bfgs = np.eye(3) - Jacobian.update_bfgs(bfgs, x + dx, dx, f2(x + dx) - f2(x)) - bfgs = np.linalg.inv(bfgs) - assert np.allclose(df, bfgs.dot(dx), atol=1e-1) diff --git a/flik/test/test_newton.py b/flik/test/test_newton.py deleted file mode 100644 index 50e9a7d..0000000 --- a/flik/test/test_newton.py +++ /dev/null @@ -1,234 +0,0 @@ -# An experimental local optimization package -# Copyright (C) 2018 Ayers Lab . -# -# This file is part of Flik. -# -# Flik 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. -# -# Flik 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 - - -"""Test file for `flik.nonlinear.nonlinear_solve`.""" - - -import numpy as np - -from numpy.testing import assert_raises - -from flik import nonlinear_solve - - -__all__ = [ - "test_nonlinear_inputs", - "test_nonlinear_linear_solve", - "test_nonlinear_no_solution1", - "test_nonlinear_nonlinear_solve", - "test_nonlinear_singular_matrix_error", - "test_nonlinear_goodbroyden", - "test_nonlinear_badbroyden", - "test_nonlinear_sr1", - "test_nonlinear_sr1inv", - "test_nonlinear_dfp", - ] - - -def f1(_): - """Test function.""" - return np.zeros((2, 1)) - - -def j1(_): - """Test function.""" - return np.array([[1, 1], [-1, 1]]) - - -def f4(x): - """Test function.""" - return np.array([x[0] + x[1] - 3, x[1] - x[0] + 1]) - - -def f5(x): - """Test function.""" - return np.array([np.power(x[0], 3) + x[1] - 1, np.power(x[1], 3) - x[0] + 1]) - - -def j5(x): - """Test function.""" - return np.array([[3 * x[0] ** 2, 1], [-1, 3 * x[1] ** 2]]) - - -def f6(x): - """Test function.""" - return np.array([x[0] ** 2 + x[1] ** 2 - 1, -x[0] ** 2 + x[1] + 10]) - - -def j6(x): - """Test function.""" - return np.array([[2 * x[0], 2 * x[1]], [-2 * x[0], 1]]) - - -def f7(x): - """Test function.""" - return np.array([-2. * x[0] ** 2. - (4. / 3.) * x[1] ** 3., - -2. * x[0] ** 2. - (4./3.) * x[1] ** 3.]) - - -def j7(x): - """Test function.""" - return np.array([[4. * x[0], 4. * x[1] ** 2.], - [4. * x[0], 4 * x[1]**2.]]) - - -def test_nonlinear_inputs(): - """Test invalid inputs to nonlinear_solve.""" - # Initial guess - x_0 = np.array([1.0, 1.0]) - # Test function coverage - f1(x_0) - # Check validity of inputs - assert_raises(TypeError, nonlinear_solve, "string", x_0, J=j1) - assert_raises(TypeError, nonlinear_solve, f1, "string", J=j1) - assert_raises(TypeError, nonlinear_solve, f1, x_0, J="string") - assert_raises(TypeError, nonlinear_solve, f1, x_0, J=j1, stepsize="string") - assert_raises(TypeError, nonlinear_solve, f1, x_0, J=j1, stepsize=np.ones((2, 2))) - assert_raises(TypeError, nonlinear_solve, f1, x_0, J=j1, eps="string") - assert_raises(ValueError, nonlinear_solve, f1, x_0, J=j1, eps=-1.0) - assert_raises(ValueError, nonlinear_solve, f1, x_0, J=j1, eps=-1) - assert_raises(ValueError, nonlinear_solve, f1, x_0, J=j1, method="tomato") - assert_raises(ValueError, nonlinear_solve, f1, x_0, J=j1, stepsize=-0.1) - assert_raises(ValueError, nonlinear_solve, f1, x_0, J=j1, stepsize=-np.abs(x_0)) - assert_raises(TypeError, nonlinear_solve, f1, x_0, J=j1, maxiter=0.) - assert_raises(ValueError, nonlinear_solve, f1, x_0, J=j1, maxiter=-1) - assert_raises(TypeError, nonlinear_solve, f1, x_0, J=j1, method=0) - - -def test_nonlinear_linear_solve(): - """Test that newton solves linear systems in 1 step.""" - x_0 = np.array([1.0, 1.0]) - result = nonlinear_solve(f4, x_0, j1, stepsize=np.array([1., 1.]), eps=1.0e-9, maxiter=1) - assert result["success"] - assert result["message"] == "Convergence obtained." - assert result["niter"] == 1 - assert np.allclose(result["f"], [0., 0.], atol=1.0e-6) - assert np.allclose(result["x"], [2., 1.], atol=1.0e-9) - assert np.allclose(result["f"], [0., 0.], atol=1.0e-6) - assert np.allclose(result["J"], [[1., 1.], [-1., 1]], atol=1.0e-6) - - -def test_nonlinear_nonlinear_solve(): - """Test that newton solves nonlinear systems.""" - x_0 = np.array([0.5, 0.5]) - result = nonlinear_solve(f5, x_0, j5, eps=1.0e-9, maxiter=100) - assert result["success"] - assert result["niter"] < 101 - assert np.allclose(result["f"], [0., 0.], atol=1.0e-6) - assert result["message"] == "Convergence obtained." - assert result["niter"] < 100 - assert np.allclose(result["x"], [1., 0.], atol=1.0e-6) - assert np.allclose(result["f"], [0., 0.], atol=1.0e-6) - assert np.allclose(result["J"], [[3., 1.], [-1., 0]], atol=1.0e-6) - - -def test_nonlinear_no_solution1(): - """Test that no solution in nonlinear_solve raises error.""" - x_0 = np.array([1., 1.]) - result = nonlinear_solve(f6, x_0, j6, eps=1.0e-9, maxiter=100) - assert not result["success"] - assert not np.allclose(result["f"], [0., 0.], atol=1.0e-3) - assert result["message"] == "Maximum number of iterations reached." - assert result["niter"] == 100 - - -def test_nonlinear_singular_matrix_error(): - """Test that newton raises error when Jacobian is singular.""" - x0 = np.array([5., 5.]) - result = nonlinear_solve(f7, x0, j7) - message_expt = "Singular Jacobian; no solution found." - assert not result['success'] - assert result["message"] == message_expt - assert result["eps"] == 1e-6 - - -def test_nonlinear_goodbroyden(): - """Test that newton solves system with good broyden update.""" - x_0 = np.array([0.5, 0.5]) - result = nonlinear_solve(f5, x_0, stepsize=1, eps=1.0e-6, maxiter=100, method="goodbroyden") - assert result["success"] - assert result["niter"] < 101 - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-5) - assert result["message"] == "Convergence obtained." - assert result["niter"] < 100 - assert np.allclose(result["x"], [1., 0.], rtol=1.0e-5, atol=1.0e-5) - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-5) - assert np.allclose(result["J"], [[3., 1.], [-1., 0]], rtol=1.0e-3, atol=1.0e-3) - - -def test_nonlinear_badbroyden(): - """Test that newton solves system with bad broyden update.""" - x_0 = np.array([1.1, 0.5]) - result = nonlinear_solve(f5, x_0, stepsize=1, eps=1.0e-6, maxiter=100, method="badbroyden") - assert result["success"] - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-5) - assert result["message"] == "Convergence obtained." - assert result["niter"] < 100 - assert np.allclose(result["x"], [1., 0.], rtol=1.0e-5, atol=1.0e-5) - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-5) - # assert np.allclose(result["J"], [[3., 1.], [-1., 0]], rtol=1.0e-4, atol=1.0e-2) - - -def test_nonlinear_dfp(): - """Test that newton solves system with dfp update.""" - x_0 = np.array([3., 1.]) - result = nonlinear_solve(f5, x_0, j5, stepsize=1, eps=1.0e-5, maxiter=1000, method="dfp") - assert result["success"] - assert result["message"] == "Convergence obtained." - assert np.allclose(result["f"], [0., 0.], atol=1.0e-3) - assert np.allclose(result["x"], [1., 0.], atol=1.0e-3) - assert result["niter"] < 1000 - - -def test_nonlinear_sr1(): - """Test that newton solves system with sr1 update.""" - x_0 = np.array([0.5, 0.5]) - result = nonlinear_solve(f5, x_0, j5, stepsize=0.5, eps=1.0e-6, maxiter=100, method="sr1") - assert result["success"] - assert result["niter"] < 101 - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-5) - assert result["message"] == "Convergence obtained." - assert result["niter"] < 100 - assert np.allclose(result["x"], [1., 0.], rtol=1.0e-5, atol=1.0e-5) - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-5) - # assert np.allclose(result["J"], [[3., 1.], [-1., 0]], rtol=1.0e-3, atol=1.0e-2) - - -def test_nonlinear_sr1inv(): - """Test that newton solves system with sr1 inverse update.""" - x_0 = np.array([0.5, 0.5]) - result = nonlinear_solve(f5, x_0, j5, stepsize=0.5, eps=1.0e-6, maxiter=1000, method="sr1inv") - assert result["success"] - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-5) - assert result["message"] == "Convergence obtained." - assert result["niter"] < 100 - assert np.allclose(result["x"], [1., 0.], rtol=1.0e-5, atol=1.0e-5) - # assert np.allclose(result["J"], [[3., 1.], [-1., 0]], rtol=1.0e-5, atol=1.0e-1) - - -def test_nonlinear_bfgs(): - """Test that newton solves system with bfgs inverse update.""" - x_0 = np.array([0.5, 0.5]) - result = nonlinear_solve(f5, x_0, j5, stepsize=0.5, eps=1.0e-5, maxiter=10000, method="bfgs") - assert result["success"] - assert result["niter"] < 10000 - assert np.allclose(result["f"], [0., 0.], rtol=1.0e-5, atol=1.0e-3) - assert result["message"] == "Convergence obtained." - assert np.allclose(result["x"], [1., 0.], rtol=1.0e-5, atol=1.0e-2) - # assert np.allclose(result["J"], j5([1., 0.]), atol=1e-1) diff --git a/flik/trustregion.py b/flik/trustregion.py new file mode 100644 index 0000000..fb0ac14 --- /dev/null +++ b/flik/trustregion.py @@ -0,0 +1,13 @@ +class TrustRegion: + def __init__(self, radius): + self.radius = radius + # or whatever other structure for storing data + self.cache = {} + + def update(self, x, step, model): + # change self.radius + pass + + def check_step(self, step): + # check if step should be accepted + return True