diff --git a/pygradflow/implicit_func.py b/pygradflow/implicit_func.py index 912c3a4..71d9617 100644 --- a/pygradflow/implicit_func.py +++ b/pygradflow/implicit_func.py @@ -116,6 +116,25 @@ def active_set_at_point(self, p): ub = problem.var_ub return super().compute_active_set_box(p, lb, ub) + def flow_rhs(self, iterate: Iterate, rho: float) -> np.ndarray: + + lb = self.problem.var_lb + ub = self.problem.var_ub + x = iterate.x + + xflow = -iterate.aug_lag_deriv_x(rho) + + # Project flow onto tnagent cone of box constraints at given x + at_lb = np.where(x < lb - 1e-8)[0] + at_ub = np.where(x > ub + 1e-8)[0] + + xflow[at_lb] = np.maximum(0.0, xflow[at_lb]) + xflow[at_ub] = np.minimum(0.0, xflow[at_ub]) + + yflow = iterate.aug_lag_deriv_y() + + return np.concatenate([xflow, yflow]) + def projection_initial(self, iterate: Iterate, rho: float, tau=None): x_0 = self.orig_iterate.x dt = self.dt diff --git a/pygradflow/step/residuum_ratio_control.py b/pygradflow/step/residuum_ratio_control.py index 33737c3..7253afd 100644 --- a/pygradflow/step/residuum_ratio_control.py +++ b/pygradflow/step/residuum_ratio_control.py @@ -6,6 +6,7 @@ from pygradflow.params import Params from pygradflow.problem import Problem from pygradflow.step.newton_control import NewtonController +from pygradflow.step.solver.step_solver import StepResult from pygradflow.step.step_control import StepControlResult @@ -29,7 +30,8 @@ def step(self, iterate, rho, dt, display, timer): mid_step = next(next_steps) mid_iterate = mid_step.iterate - mid_norm = np.linalg.norm(func.value_at(mid_iterate, rho)) + mid_value = func.value_at(mid_iterate, rho) + mid_norm = np.linalg.norm(mid_value) self.display_step(0, mid_step) @@ -43,21 +45,81 @@ def step(self, iterate, rho, dt, display, timer): theta = mid_norm / orig_norm accepted = theta <= params.theta_max - if accepted: - lamb_mod = self.controller.update(theta) - lamb_n = max(params.lamb_min, lamb / lamb_mod) - else: - lamb_n = lamb * params.lamb_inc - if self.controller.error_sum > 0.0: - self.controller.reset() - logger.debug( - "StepController: theta: %e, accepted: %s, lamb: %e, lamb_n: %e", + "StepController: theta: %e, accepted: %s, lamb: %e", theta, accepted, lamb, - lamb_n, ) - self.lamb = lamb_n - return StepControlResult.from_step_result(mid_step, lamb_n, accepted) + if accepted: + lamb_mod = self.controller.update(theta) + lamb_n = max(params.lamb_min, lamb / lamb_mod) + self.lamb = lamb_n + return StepControlResult.from_step_result(mid_step, lamb_n, accepted) + + # Step would be rejected + lamb_n = lamb * params.lamb_inc + + if self.controller.error_sum > 0.0: + self.controller.reset() + + # Recovery starts here + + active_set = func.compute_active_set(iterate, rho) + mid_active_set = func.compute_active_set(mid_iterate, rho) + + if (active_set == mid_active_set).all(): + return StepControlResult.from_step_result(mid_step, lamb_n, accepted) + + mid_primal_value = mid_value[: self.problem.num_vars] + + mid_primal_norm = np.linalg.norm(mid_primal_value) + + mid_unchanged_primal_value = mid_primal_value[mid_active_set == active_set] + mid_unchanged_primal_norm = np.linalg.norm(mid_unchanged_primal_value) + + norm_factor = mid_unchanged_primal_norm / mid_primal_norm + + if norm_factor >= 1e-6: + self.lamb = lamb_n + return StepControlResult.from_step_result(mid_step, lamb_n, accepted) + + dir = func.flow_rhs(mid_iterate, rho) + + # Track back until first active step change... + + primal_dir = dir[: self.problem.num_vars] + dual_dir = dir[self.problem.num_vars :] + + pos_dir = primal_dir > 0.0 + neg_dir = primal_dir < 0.0 + + lb = problem.var_lb + ub = problem.var_ub + + pos_ratio = (ub[pos_dir] - mid_primal_value[pos_dir]) / primal_dir[pos_dir] + neg_ratio = (lb[neg_dir] - mid_primal_value[neg_dir]) / primal_dir[neg_dir] + + recovery_dt = np.minimum(np.min(pos_ratio), np.min(neg_ratio)) + recovery_dt = np.minimum(dt, recovery_dt) + + recovery_dx = recovery_dt * primal_dir + recovery_dy = recovery_dt * dual_dir + + recovery_step = StepResult(iterate, recovery_dx, recovery_dy, active_set) + + recovery_func = ImplicitFunc(problem, iterate, recovery_dt) + + recovery_value = recovery_func.value_at(recovery_step.iterate, rho) + recovery_norm = np.linalg.norm(recovery_value) + + recovery_theta = recovery_norm / orig_norm + accepted = recovery_theta <= params.theta_max + + if accepted: + lamb_n = lamb + else: + lamb_n = 2.0 * lamb + + return StepControlResult.from_step_result(recovery_step, lamb_n, accepted=accepted) diff --git a/tests/pygradflow/linprog.py b/tests/pygradflow/linprog.py new file mode 100644 index 0000000..cb56edf --- /dev/null +++ b/tests/pygradflow/linprog.py @@ -0,0 +1,48 @@ +import numpy as np +import scipy as sp + +from pygradflow.problem import Problem + + +class LinearProgram(Problem): + """ + Linear program: minimize c^T x subject to Ax <= b + """ + + def __init__(self, A, b, c): + (num_cons, num_vars) = A.shape + + var_lb = np.full(shape=(num_vars,), fill_value=-np.inf) + var_ub = np.full(shape=(num_vars,), fill_value=np.inf) + + cons_lb = np.full(shape=(num_cons,), fill_value=-np.inf) + cons_ub = b + + self.c = c + self.b = b + self.A = sp.sparse.coo_matrix(A) + + super().__init__(var_lb, var_ub, cons_lb=cons_lb, cons_ub=cons_ub) + + def obj(self, x): + return np.dot(self.c, x) + + def obj_grad(self, x): + return self.c + + def cons(self, x): + return self.A.dot(x) + + def cons_jac(self, x): + return self.A + + def lag_hess(self, x, y): + return sp.sparse.coo_matrix( + np.zeros( + ( + self.num_vars, + self.num_vars, + ), + dtype=float, + ) + ) diff --git a/tests/pygradflow/test_solve_linear.py b/tests/pygradflow/test_solve_linear.py new file mode 100644 index 0000000..2458aeb --- /dev/null +++ b/tests/pygradflow/test_solve_linear.py @@ -0,0 +1,154 @@ +import numpy as np +import pytest + +from pygradflow.integration.integration_solver import IntegrationSolver +from pygradflow.params import Params, PenaltyUpdate, StepControlType +from pygradflow.solver import Solver + +from .linprog import LinearProgram + + +@pytest.fixture +def linprog(): + A = np.array( + [ + [9.51056516e-01, 6.90983006e-01], + [-3.63271264e-01, 1.11803399e00], + [-1.17557050e00, 2.22044605e-16], + [-3.63271264e-01, -1.11803399e00], + [9.51056516e-01, -6.90983006e-01], + [-2.17234844e-01, 1.07058388e00], + [-7.11984239e-02, 1.02313377e00], + [7.48379961e-02, 9.75683661e-01], + [2.20874416e-01, 9.28233552e-01], + [3.66910836e-01, 8.80783443e-01], + [5.12947256e-01, 8.33333333e-01], + [6.58983676e-01, 7.85883224e-01], + [8.05020096e-01, 7.38433115e-01], + [-1.08531503e00, 1.24225999e-01], + [-9.95059562e-01, 2.48451997e-01], + [-9.04804091e-01, 3.72677996e-01], + [-8.14548620e-01, 4.96903995e-01], + [-7.24293149e-01, 6.21129994e-01], + [-6.34037678e-01, 7.45355992e-01], + [-5.43782206e-01, 8.69581991e-01], + [-4.53526735e-01, 9.93807990e-01], + [-4.53526735e-01, -9.93807990e-01], + [-5.43782206e-01, -8.69581991e-01], + [-6.34037678e-01, -7.45355992e-01], + [-7.24293149e-01, -6.21129994e-01], + [-8.14548620e-01, -4.96903995e-01], + [-9.04804091e-01, -3.72677996e-01], + [-9.95059562e-01, -2.48451997e-01], + [-1.08531503e00, -1.24225999e-01], + [8.05020096e-01, -7.38433115e-01], + [6.58983676e-01, -7.85883224e-01], + [5.12947256e-01, -8.33333333e-01], + [3.66910836e-01, -8.80783443e-01], + [2.20874416e-01, -9.28233552e-01], + [7.48379961e-02, -9.75683661e-01], + [-7.11984239e-02, -1.02313377e00], + [-2.17234844e-01, -1.07058388e00], + [9.51056516e-01, 5.37431227e-01], + [9.51056516e-01, 3.83879448e-01], + [9.51056516e-01, 2.30327669e-01], + [9.51056516e-01, 7.67758895e-02], + [9.51056516e-01, -7.67758895e-02], + [9.51056516e-01, -2.30327669e-01], + [9.51056516e-01, -3.83879448e-01], + [9.51056516e-01, -5.37431227e-01], + ] + ) + + b = np.array( + [ + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + 0.95105652, + ] + ) + + c = np.array([0.0, 0.1]) + + return LinearProgram(A, b, c) + + +@pytest.fixture +def x0(): + return np.array([0.0, 0.5]) + + +@pytest.fixture +def y0(): + return 0.0 + + +def test_solve_with_recovery(linprog, x0, y0): + params = Params( + collect_path=True, + rho=1e-1, + step_control_type=StepControlType.ResiduumRatio, + penalty_update=PenaltyUpdate.Constant, + ) + + solver = Solver(linprog, params) + + result = solver.solve(x0, y0) + + import matplotlib.pyplot as plt + + primal_path = result.primal_path + + path_x = primal_path[0, :] + path_y = primal_path[1, :] + + # plt.plot(path_x, path_y, "o-") + # plt.show() + + # import pdb + + # pdb.set_trace() + + assert result.success