diff --git a/pygradflow/newton.py b/pygradflow/newton.py index 5c6217c..1a15528 100644 --- a/pygradflow/newton.py +++ b/pygradflow/newton.py @@ -8,7 +8,7 @@ from pygradflow.params import Params, NewtonType from pygradflow.problem import Problem from pygradflow.step import step_solver -from pygradflow.step.step_solver import StepSolver +from pygradflow.step.step_solver import StepSolver, StepResult logger = lgg.getChild("newton") @@ -29,7 +29,7 @@ def params(self) -> Params: return self.orig_iterate.params @abc.abstractmethod - def step(self, iterate): + def step(self, iterate: Iterate) -> StepResult: raise NotImplementedError() @@ -56,14 +56,8 @@ def __init__( self.step_solver.update_active_set(active_set) self.step_solver.update_derivs(orig_iterate) - def step(self, iterate: Iterate) -> Iterate: - (xn, yn) = self.step_solver.solve(iterate) - - return Iterate(self.problem, - self.params, - xn, - yn, - self.orig_iterate.eval) + def step(self, iterate): + return self.step_solver.solve(iterate) class FullNewtonMethod(NewtonMethod): @@ -84,20 +78,14 @@ def __init__( super().__init__(problem, orig_iterate, dt, rho) self.step_solver = step_solver - def step(self, iterate: Iterate) -> Iterate: + def step(self, iterate: Iterate) -> StepResult: p = self.func.projection_initial(iterate, self.rho) active_set = self.func.compute_active_set(p) self.step_solver.update_active_set(active_set) self.step_solver.update_derivs(iterate) - (xn, yn) = self.step_solver.solve(iterate) - - return Iterate(self.problem, - self.params, - xn, - yn, - self.orig_iterate.eval) + return self.step_solver.solve(iterate) class FixedActiveSetNewtonMethod(NewtonMethod): @@ -121,23 +109,22 @@ def __init__(self, problem, active_set, orig_iterate, dt, rho): np.sum(active_set), ) - def split_sol(self, s): + def split_sol(self, s: np.ndarray): n = self.problem.num_vars m = self.problem.num_cons assert s.shape == (n + m,) return s[:n], s[n:] - def create_iterate(self, iterate, s): + def create_step(self, iterate: Iterate, s: np.ndarray) -> StepResult: xn, yn = self.split_sol(s) x = iterate.x y = iterate.y - return Iterate(self.problem, - self.params, - x - xn, - y - yn, - self.orig_iterate.eval) + dx = x - xn + dy = y - yn + + return StepResult(iterate, dx, dy) @staticmethod def active_set_from_iterate(problem, iterate): @@ -177,7 +164,8 @@ def step(self, iterate): solver = sp.sparse.linalg.splu(mat) s = solver.solve(rhs) - next_iterate = self.create_iterate(iterate, s) + next_step = self.create_step(iterate, s) + next_iterate = next_step.iterate logger.info( "Initial rhs norm: {0}, final: {1}".format( @@ -185,7 +173,7 @@ def step(self, iterate): ) ) - return next_iterate + return next_step class ActiveSetNewtonMethod(NewtonMethod): @@ -209,19 +197,13 @@ def __init__( self.step_solver = step_solver self.step_solver.update_derivs(orig_iterate) - def step(self, iterate: Iterate) -> Iterate: + def step(self, iterate): p = self.func.projection_initial(iterate, self.rho) active_set = self.func.compute_active_set(p) self.step_solver.update_active_set(active_set) - (xn, yn) = self.step_solver.solve(iterate) - - return Iterate(self.problem, - self.params, - xn, - yn, - self.orig_iterate.eval) + return self.step_solver.solve(iterate) def newton_method( diff --git a/pygradflow/solver.py b/pygradflow/solver.py index 7c5e77a..062094e 100644 --- a/pygradflow/solver.py +++ b/pygradflow/solver.py @@ -82,14 +82,14 @@ def compute_step( method = newton_method(problem, params, iterate, dt, self.rho) - def next_iterates(): + def next_steps(): curr_iterate = iterate while True: - next_iterate = method.step(curr_iterate) - yield next_iterate - curr_iterate = next_iterate + next_step = method.step(curr_iterate) + yield next_step + curr_iterate = next_step.iterate - return controller.step(iterate, self.rho, dt, next_iterates()) + return controller.step(iterate, self.rho, dt, next_steps()) def _deriv_check(self, x: np.ndarray, y: np.ndarray) -> None: from pygradflow.deriv_check import deriv_check diff --git a/pygradflow/step/scaled_step_solver.py b/pygradflow/step/scaled_step_solver.py index 95a1610..fae0368 100644 --- a/pygradflow/step/scaled_step_solver.py +++ b/pygradflow/step/scaled_step_solver.py @@ -4,7 +4,7 @@ import numpy as np from pygradflow.implicit_func import ScaledImplicitFunc -from pygradflow.step.step_solver import StepSolver +from pygradflow.step.step_solver import StepSolver, StepResult from pygradflow.iterate import Iterate from pygradflow.params import Params from pygradflow.problem import Problem @@ -68,7 +68,7 @@ def update_active_set(self, active_set: np.ndarray) -> None: self.active_set = copy.copy(active_set) self.reset_deriv() - def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]: + def solve(self, iterate: Iterate) -> StepResult: (b0, b1, b2) = self.initial_rhs(iterate) n = self.n @@ -90,7 +90,4 @@ def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]: dx = sx dy = fact * (sy - rho * b2) - x = iterate.x - y = iterate.y - - return (x - dx, y - dy) + return StepResult(iterate, dx, dy) diff --git a/pygradflow/step/standard_step_solver.py b/pygradflow/step/standard_step_solver.py index 1d806d8..9e0d27e 100644 --- a/pygradflow/step/standard_step_solver.py +++ b/pygradflow/step/standard_step_solver.py @@ -4,7 +4,7 @@ import numpy as np from pygradflow.implicit_func import ImplicitFunc -from pygradflow.step.step_solver import StepSolver +from pygradflow.step.step_solver import StepSolver, StepResult from pygradflow.iterate import Iterate from pygradflow.params import Params from pygradflow.problem import Problem @@ -54,7 +54,7 @@ def update_active_set(self, active_set: np.ndarray) -> None: self.active_set = copy.copy(active_set) self._reset_deriv() - def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]: + def solve(self, iterate: Iterate) -> StepResult: if self.deriv is None: self._compute_deriv() @@ -73,7 +73,4 @@ def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]: dx = s[:n] dy = s[n:] - x = iterate.x - y = iterate.y - - return (x - dx, y - dy) + return StepResult(iterate, dx, dy) diff --git a/pygradflow/step/step_control.py b/pygradflow/step/step_control.py index a82243d..948bac0 100644 --- a/pygradflow/step/step_control.py +++ b/pygradflow/step/step_control.py @@ -5,6 +5,7 @@ from pygradflow.log import logger +from pygradflow.step.step_solver import StepResult from pygradflow.controller import ControllerSettings, LogController from pygradflow.implicit_func import ImplicitFunc from pygradflow.iterate import Iterate @@ -12,7 +13,7 @@ from pygradflow.problem import Problem -class StepResult: +class StepControlResult: def __init__(self, iterate: Iterate, lamb: float, accepted: bool) -> None: self.iterate = iterate self.lamb = lamb @@ -26,7 +27,11 @@ def __init__(self, problem: Problem, params: Params) -> None: self.lamb = params.lamb_init @abc.abstractmethod - def step(self, iterate, rho, dt, next_iterates): + def step(self, + iterate: Iterate, + rho: float, + dt: float, + next_steps: Iterator[StepResult]) -> StepControlResult: raise NotImplementedError() @@ -34,7 +39,7 @@ class ExactController(StepController): def __init__(self, problem, params): super().__init__(problem, params) - def step(self, iterate, rho, dt, next_iterates): + def step(self, iterate, rho, dt, next_steps): assert dt > 0.0 lamb = 1.0 / dt @@ -46,20 +51,20 @@ def func_val(iterate): cur_func_val = func_val(iterate) for i in range(10): - next_iterate = next(next_iterates) + next_iterate = next(next_steps).iterate next_func_val = func_val(next_iterate) logger.info(f"Func val: {next_func_val}") if next_func_val <= self.params.newton_tol: logger.debug("Newton method converged in %d iterations", i + 1) - return StepResult(next_iterate, 0.5 * lamb, True) + return StepControlResult(next_iterate, 0.5 * lamb, True) elif next_func_val > cur_func_val: break logger.debug("Newton method did not converge") - return StepResult(next_iterate, 2.0 * lamb, False) + return StepControlResult(next_iterate, 2.0 * lamb, False) class DistanceRatioController(StepController): @@ -68,9 +73,7 @@ def __init__(self, problem: Problem, params: Params) -> None: settings = ControllerSettings.from_params(params) self.controller = LogController(settings, params.theta_ref) - def step( - self, iterate: Iterate, rho: float, dt: float, next_iterates: Iterator[Iterate] - ) -> StepResult: + def step(self, iterate, rho, dt, next_steps): assert dt > 0.0 lamb = 1.0 / dt @@ -79,24 +82,30 @@ def step( func = ImplicitFunc(problem, iterate, dt) - mid_iterate = next(next_iterates) + # TODO: Fix + mid_step = next(next_steps) + mid_iterate = mid_step.iterate - if np.linalg.norm(func.value_at(mid_iterate, rho)) <= params.newton_tol: + mid_func_norm = np.linalg.norm(func.value_at(mid_iterate, rho)) + + if mid_func_norm <= params.newton_tol: lamb_n = max(lamb * params.lamb_red, params.lamb_min) - logger.info("Newton converged during first iteration, lamb_n = %f", lamb_n) - return StepResult(mid_iterate, lamb_n, True) + logger.info("Newton converged during first iteration, lamb_n = %f", + lamb_n) + return StepControlResult(mid_iterate, lamb_n, True) - first_diff = mid_iterate.dist(iterate) + first_diff = mid_step.diff if first_diff == 0.: - return StepResult(mid_iterate, lamb, True) + return StepControlResult(mid_iterate, lamb, True) - final_iterate = next(next_iterates) + final_step = next(next_steps) + final_iterate = final_step.iterate - second_diff = final_iterate.dist(mid_iterate) + second_diff = final_step.diff if second_diff == 0.: - return StepResult(final_iterate, lamb, True) + return StepControlResult(final_iterate, lamb, True) theta = second_diff / first_diff @@ -120,4 +129,4 @@ def step( self.lamb = lamb_n - return StepResult(final_iterate, lamb_n, accepted) + return StepControlResult(final_iterate, lamb_n, accepted) diff --git a/pygradflow/step/step_solver.py b/pygradflow/step/step_solver.py index 8153e9b..864eced 100644 --- a/pygradflow/step/step_solver.py +++ b/pygradflow/step/step_solver.py @@ -1,4 +1,5 @@ import abc +import functools import numpy as np import scipy as sp @@ -7,6 +8,33 @@ from pygradflow.params import Params from pygradflow.problem import Problem from pygradflow.step.linear_solver import LinearSolver +from pygradflow.util import norm_mult + + +class StepResult: + def __init__(self, orig_iterate, dx, dy): + self.orig_iterate = orig_iterate + self.dx = dx + self.dy = dy + + @functools.cached_property + def iterate(self): + iterate = self.orig_iterate + + return Iterate(iterate.problem, + iterate.params, + iterate.x - self.dx, + iterate.y - self.dy, + iterate.eval) + + @functools.cached_property + def diff(self): + diff = norm_mult(self.dx, self.dy) + + assert np.allclose(self.iterate.dist(self.orig_iterate), + diff) + + return diff class StepSolver(abc.ABC): @@ -24,10 +52,10 @@ def linear_solver(self, mat: sp.sparse.spmatrix) -> LinearSolver: @abc.abstractmethod def update_active_set(self, active_set: np.ndarray): - raise NotImplementedError + raise NotImplementedError() def update_derivs(self, iterate: Iterate): - raise NotImplementedError + raise NotImplementedError() - def solve(self, iterate: Iterate): - raise NotImplementedError + def solve(self, iterate: Iterate) -> StepResult: + raise NotImplementedError() diff --git a/pyproject.toml b/pyproject.toml index ea5d7d6..05823d7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "pygradflow" -version = "0.2.1" +version = "0.2.2" description = "PyGradFlow is a simple implementation of the sequential homotopy method to be used to solve general nonlinear programs." authors = ["Christoph Hansknecht "] readme = "README.md" diff --git a/tests/pygradflow/test_solver.py b/tests/pygradflow/test_solver.py index 5386c7d..6a7784e 100644 --- a/tests/pygradflow/test_solver.py +++ b/tests/pygradflow/test_solver.py @@ -180,7 +180,7 @@ def test_newton_step_unconstrained( assert np.allclose(deriv.toarray(), np.eye(s)) - next_iterate = newton.step(iterate) + next_iterate = newton.step(iterate).iterate assert np.allclose(next_iterate.x, iterate.x) assert np.allclose(next_iterate.y, iterate.y) @@ -219,7 +219,7 @@ def test_newton_step_constrained( assert np.allclose(deriv.toarray(), np.eye(s)) - next_iterate = newton.step(iterate) + next_iterate = newton.step(iterate).iterate xclip = np.clip(x_0, problem.var_lb, problem.var_ub) @@ -253,7 +253,7 @@ def test_custom_step_solver(rosenbrock_instance): assert np.allclose(deriv.toarray(), np.eye(s)) - next_iterate = newton.step(iterate) + next_iterate = newton.step(iterate).iterate xclip = np.clip(x_0, problem.var_lb, problem.var_ub) @@ -338,7 +338,7 @@ def test_one_step_convergence(newton_type, step_solver_type, linear_solver_type) iterate = Iterate(problem, params, x_0, y_0) method = newton_method(problem, params, iterate, dt, rho) - next_iterate = method.step(iterate) + next_iterate = method.step(iterate).iterate func = ImplicitFunc(problem, iterate, dt) assert np.allclose(func.value_at(next_iterate, rho=rho), 0.0)