Subgradient hub

doc/src/hubs.rst

@@ -64,6 +64,14 @@ APH
 The implementation of Asynchronous Projective Hedging is described in a
 forthcoming paper.
 
+Subgradient
+-----------
+
+The Subgradient implemenation can be used with most spokes becuase it
+also supplies x and/or W values at every iteration, and is largely based
+on the PH implementation. It utilizes a constant step size rule based on
+`rho` unless modified by an extension.
+
 Hub Convergers
 --------------
 

examples/netdes/netdes_cylinders.py

@@ -22,9 +22,10 @@ def _parse_args():
     cfg.two_sided_args()
     cfg.xhatlooper_args()
     cfg.ph_args()
+    cfg.subgradient_args()
     cfg.fwph_args()
     cfg.lagrangian_args()
-    cfg.subgradient_args()
+    cfg.subgradient_bounder_args()
     cfg.xhatshuffle_args()
     cfg.slammax_args()
     cfg.cross_scenario_cuts_args()
@@ -71,11 +72,17 @@ def main():
     else:
         ph_ext = None
 
-    # Vanilla PH hub
-    hub_dict = vanilla.ph_hub(*beans,
-                              scenario_creator_kwargs=scenario_creator_kwargs,
-                              ph_extensions=ph_ext,
-                              rho_setter = None)
+    if cfg.SUBGRAD:
+        hub_dict = vanilla.subgradient_hub(*beans,
+                                  scenario_creator_kwargs=scenario_creator_kwargs,
+                                  ph_extensions=ph_ext,
+                                  rho_setter = None)
+    else:
+        # Vanilla PH hub
+        hub_dict = vanilla.ph_hub(*beans,
+                                  scenario_creator_kwargs=scenario_creator_kwargs,
+                                  ph_extensions=ph_ext,
+                                  rho_setter = None)
 
     if cross_scenario_cuts:
         hub_dict["opt_kwargs"]["options"]["cross_scen_options"]\

examples/run_all.py

@@ -248,10 +248,10 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring):
        f"--num-scens 3 --crops-multiplier=1  --EF-solver-name={solver_name} "
        "--BPL-c0 25 --BPL-eps 100 --confidence-level 0.95 --BM-vs-BPL BPL")
 
-do_one("netdes", "netdes_cylinders.py", 5,
+do_one("netdes", "netdes_cylinders.py", 4,
        "--max-iterations=3 --instance-name=network-10-20-L-01 "
-       "--solver-name={} --rel-gap=0.0 --default-rho=1 --presolve "
-       "--slammax --lagrangian --xhatshuffle --cross-scenario-cuts --max-solver-threads=2".format(solver_name))
+       "--solver-name={} --rel-gap=0.0 --default-rho=10000 --presolve "
+       "--slammax --SUBGRAD --xhatshuffle --cross-scenario-cuts --max-solver-threads=2".format(solver_name))
 
 # sizes is slow for xpress so try linearizing the proximal term.
 do_one("sizes",

examples/sslp/sslp_cylinders.py

@@ -30,7 +30,7 @@ def _parse_args():
     cfg.fwph_args()
     cfg.lagrangian_args()
     cfg.xhatshuffle_args()
-    cfg.subgradient_args()
+    cfg.subgradient_bounder_args()
     cfg.reduced_costs_args()
     cfg.coeff_rho_args()
     cfg.integer_relax_then_enforce_args()

mpisppy/cylinders/hub.py

@@ -10,6 +10,7 @@
 import abc
 import logging
 import mpisppy.log
+from mpisppy.opt.subgradient import Subgradient
 from mpisppy.opt.aph import APH
 
 from mpisppy import MPI
@@ -469,11 +470,6 @@ def setup_hub(self):
                 "Cannot call setup_hub before memory windows are constructed"
             )
 
-        # attribute to set False if some extension
-        # modified the iteration 0 subproblems such
-        # that the trivial bound is no longer valid
-        self.use_trivial_bound = True
-
         self.initialize_spoke_indices()
         self.initialize_bound_values()
 
@@ -535,9 +531,8 @@ def sync_with_spokes(self):
         self.sync()
 
     def is_converged(self):
-        ## might as well get a bound, in this case
-        if self.opt._PHIter == 1 and self.use_trivial_bound:
-            self.BestOuterBound = self.OuterBoundUpdate(self.opt.trivial_bound)
+        if self.opt.best_bound_obj_val is not None:
+            self.BestOuterBound = self.OuterBoundUpdate(self.opt.best_bound_obj_val)
 
         if not self.has_innerbound_spokes:
             if self.opt._PHIter == 1:
@@ -553,7 +548,7 @@ def is_converged(self):
             return False
 
         if not self.has_outerbound_spokes:
-            if self.opt._PHIter == 1:
+            if self.opt._PHIter == 1 and not isinstance(self.opt, Subgradient):
                 global_toc(
                     "Without outer bound spokes, no progress "
                     "will be made on the Best Bound")

mpisppy/extensions/reduced_costs_fixer.py

@@ -90,7 +90,6 @@ def iter0_post_solver_creation(self):
                 while not self.opt.spcomm.hub_from_spoke(self.opt.spcomm.outerbound_receive_buffers[self.reduced_costs_spoke_index], self.reduced_costs_spoke_index):
                     continue
             self.sync_with_spokes(pre_iter0 = True)
-            self.opt.spcomm.use_trivial_bound = False
         self.fix_fraction_target = self._fix_fraction_target_iter0
 
     def post_iter0_after_sync(self):
@@ -109,11 +108,14 @@ def sync_with_spokes(self, pre_iter0 = False):
             self._last_serial_number = serial_number
             reduced_costs = spcomm.outerbound_receive_buffers[idx][1:1+self.nonant_length]
             this_outer_bound = spcomm.outerbound_receive_buffers[idx][0]
-            new_outer_bound = self._update_best_outer_bound(this_outer_bound)
+            is_new_outer_bound = self._update_best_outer_bound(this_outer_bound)
+            if pre_iter0:
+                # make sure we set the bound we compute prior to iteration 0
+                self.opt.spcomm.BestOuterBound = self.opt.spcomm.OuterBoundUpdate(self._best_outer_bound, idx=idx)
             if not pre_iter0 and self._use_rc_bt:
                 self.reduced_costs_bounds_tightening(reduced_costs, this_outer_bound)
             if self._use_rc_fixer and self.fix_fraction_target > 0.0:
-                if new_outer_bound or not self._rc_fixer_require_improving_lagrangian:
+                if is_new_outer_bound or not self._rc_fixer_require_improving_lagrangian:
                     self.reduced_costs_fixing(reduced_costs)
         else:
             if self.opt.cylinder_rank == 0 and self.verbose:

mpisppy/generic_cylinders.py

@@ -62,13 +62,14 @@ def _parse_args(m):
     cfg.two_sided_args()
     cfg.ph_args()
     cfg.aph_args()
+    cfg.subgradient_args()
     cfg.fixer_args()    
     cfg.integer_relax_then_enforce_args()
     cfg.gapper_args()    
     cfg.fwph_args()
     cfg.lagrangian_args()
     cfg.ph_ob_args()
-    cfg.subgradient_args()
+    cfg.subgradient_bounder_args()
     cfg.xhatshuffle_args()
     cfg.xhatxbar_args()
     cfg.converger_args()
@@ -143,6 +144,16 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_
                                    rho_setter = rho_setter,
                                    all_nodenames = all_nodenames,
                                    )
+    elif cfg.SUBGRAD:
+        # Vanilla Subgradient hub
+        hub_dict = vanilla.subgradient_hub(
+                       *beans,
+                       scenario_creator_kwargs=scenario_creator_kwargs,
+                       ph_extensions=None,
+                       ph_converger=ph_converger,
+                       rho_setter = rho_setter,
+                       all_nodenames = all_nodenames,
+                   )
     else:
         # Vanilla PH hub
         hub_dict = vanilla.ph_hub(*beans,

mpisppy/opt/aph.py

@@ -1067,6 +1067,8 @@ def APH_main(self, spcomm=None, finalize=True):
         # End APH-specific Prep
 
         trivial_bound = self.Iter0()
+        if self._can_update_best_bound():
+            self.best_bound_obj_val = trivial_bound
 
         self.setup_Lens()
         self.setup_dispatchrecord()

mpisppy/opt/ph.py

@@ -61,6 +61,8 @@ def ph_main(self, finalize=True):
         if (verbose):
             print('Calling PH Iter0 on global rank {}'.format(global_rank))
         trivial_bound = self.Iter0()
+        if self._can_update_best_bound():
+            self.best_bound_obj_val = trivial_bound
         if (verbose):
             print ('Completed PH Iter0 on global rank {}'.format(global_rank))
         if ('asynchronousPH' in self.options) and (self.options['asynchronousPH']):

mpisppy/opt/subgradient.py

@@ -0,0 +1,131 @@
+###############################################################################
+# mpi-sppy: MPI-based Stochastic Programming in PYthon
+#
+# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
+# Sustainable Energy, LLC, The Regents of the University of California, et al.
+# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
+# full copyright and license information.
+###############################################################################
+
+import mpisppy.phbase
+import mpisppy.MPI as _mpi
+
+_global_rank = _mpi.COMM_WORLD.Get_rank()
+
+class Subgradient(mpisppy.phbase.PHBase):
+    """ Subgradient Algorithm """
+
+    def subgradient_main(self, finalize=True):
+        """ Execute the subgradient algorithm.
+
+        Args:
+            finalize (bool, optional, default=True):
+                If True, call `PH.post_loops()`, if False, do not,
+                and return None for Eobj
+
+        Returns:
+            tuple:
+                Tuple containing
+
+                conv (float):
+                    The convergence value (not easily interpretable).
+                Eobj (float or `None`):
+                    If `finalize=True`, this is the expected, weighted
+                    objective value. This value is not directly useful.
+                    If `finalize=False`, this value is `None`.
+                trivial_bound (float):
+                    The "trivial bound", computed by solving the model with no
+                    nonanticipativity constraints (immediately after iter 0).
+        """
+        verbose = self.options['verbose']
+        smoothed = self.options['smoothed']
+        if smoothed != 0:
+            raise RuntimeError("Cannnot use smoothing with Subgradient algorithm")
+        self.PH_Prep(attach_prox=False, attach_smooth=smoothed)
+
+        if (verbose):
+            print('Calling Subgradient Iter0 on global rank {}'.format(_global_rank))
+        trivial_bound = self.Iter0()
+        # set self.best_bound_obj_val if we don't have any additional fixed variables
+        if self._can_update_best_bound():
+            self.best_bound_obj_val = trivial_bound
+        if (verbose):
+            print('Completed Subgradient Iter0 on global rank {}'.format(_global_rank))
+
+        self.iterk_loop()
+
+        if finalize:
+            Eobj = self.post_loops(self.extensions)
+        else:
+            Eobj = None
+
+        return self.conv, Eobj, trivial_bound
+
+    def ph_main(self, finalize=True):
+        # for working with a PHHub
+        return self.subgradient_main(finalize=finalize)
+
+    def solve_loop(self,
+                   solver_options=None,
+                   use_scenarios_not_subproblems=False,
+                   dtiming=False,
+                   dis_W=False,
+                   dis_prox=False,
+                   gripe=False,
+                   disable_pyomo_signal_handling=False,
+                   tee=False,
+                   verbose=False,
+                   need_solution=True,
+                   ):
+        """ Loop over `local_subproblems` and solve them in a manner
+        dicated by the arguments.
+
+        In addition to changing the Var values in the scenarios, this function
+        also updates the `_PySP_feas_indictor` to indicate which scenarios were
+        feasible/infeasible.
+
+        Args:
+            solver_options (dict, optional):
+                The scenario solver options.
+            use_scenarios_not_subproblems (boolean, optional):
+                If True, solves individual scenario problems, not subproblems.
+                This distinction matters when using bundling. Default is False.
+            dtiming (boolean, optional):
+                If True, reports solve timing information. Default is False.
+            dis_W (boolean, optional):
+                If True, duals weights (Ws) are disabled before solve, then
+                re-enabled after solve. Default is False.
+            dis_prox (boolean, optional):
+                If True, prox terms are disabled before solve, then
+                re-enabled after solve. Default is False.
+            gripe (boolean, optional):
+                If True, output a message when a solve fails. Default is False.
+            disable_pyomo_signal_handling (boolean, optional):
+                True for asynchronous PH; ignored for persistent solvers.
+                Default False.
+            tee (boolean, optional):
+                If True, displays solver output. Default False.
+            verbose (boolean, optional):
+                If True, displays verbose output. Default False.
+            need_solution (boolean, optional):
+                If True, raises an exception if a solution is not available.
+                Default True
+        """
+        super().solve_loop(
+            solver_options=solver_options,
+            use_scenarios_not_subproblems=use_scenarios_not_subproblems,
+            dtiming=dtiming,
+            dis_W=dis_W,
+            dis_prox=dis_prox,
+            gripe=gripe,
+            disable_pyomo_signal_handling=disable_pyomo_signal_handling,
+            tee=tee,
+            verbose=verbose,
+            need_solution=need_solution,
+        )
+
+        # set self.best_bound_obj_val if we don't have any additional fixed variables
+        if self._can_update_best_bound():
+            self.best_bound_obj_val = self.Ebound(verbose)
+
+

mpisppy/phbase.py

@@ -931,6 +931,10 @@ def _vb(msg):
         if have_extensions:
             self.extobject.post_iter0()
 
+        self.trivial_bound = self.Ebound(verbose)
+        if self._can_update_best_bound():
+            self.best_bound_obj_val = self.trivial_bound
+
         if self.spcomm is not None:
             self.spcomm.sync()
 
@@ -955,8 +959,6 @@ def _vb(msg):
         #global_toc('Rank: {} - Before iter loop'.format(self.cylinder_rank), True)
         self.conv = None
 
-        self.trivial_bound = self.Ebound(verbose)
-
         if dprogress and self.cylinder_rank == 0:
             print("")
             print("After PH Iteration",self._PHIter)
@@ -1000,6 +1002,10 @@ def iterk_loop(self):
         self.conv = None
 
         max_iterations = int(self.options["PHIterLimit"])
+        if self.spcomm is not None:
+            # print a screen trace for iteration 0
+            if self.spcomm.is_converged():
+                global_toc("Cylinder convergence", self.cylinder_rank == 0)
 
         for self._PHIter in range(1, max_iterations+1):
             iteration_start_time = time.time()
@@ -1159,6 +1165,5 @@ def attach_xbars(self):
                 scenario._mpisppy_data.nonant_indices.keys(), initialize=0.0, mutable=True
             )
 
-
 if __name__ == "__main__":
     print ("No main for PHBase")

mpisppy/spbase.py

@@ -97,6 +97,9 @@ def __init__(
         self.first_stage_solution_available = False
         self.best_solution_obj_val = None
 
+        # sometimes we know a best bound
+        self.best_bound_obj_val = None
+
         if options.get("toc", True):
             global_toc("Initializing SPBase")