diff --git a/doc/src/hubs.rst b/doc/src/hubs.rst index 5df6f2251..45bcfbe8b 100644 --- a/doc/src/hubs.rst +++ b/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 -------------- diff --git a/examples/netdes/netdes_cylinders.py b/examples/netdes/netdes_cylinders.py index 232a775aa..cd2ba9b2b 100644 --- a/examples/netdes/netdes_cylinders.py +++ b/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.subgradient_hub: + 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"]\ diff --git a/examples/run_all.py b/examples/run_all.py index 1f93be2e4..42b91b724 100644 --- a/examples/run_all.py +++ b/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 --subgradient-hub --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", diff --git a/examples/sslp/sslp_cylinders.py b/examples/sslp/sslp_cylinders.py index 139053106..bf7623886 100644 --- a/examples/sslp/sslp_cylinders.py +++ b/examples/sslp/sslp_cylinders.py @@ -34,7 +34,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() diff --git a/mpisppy/cylinders/hub.py b/mpisppy/cylinders/hub.py index ebbbb0a22..85fdee486 100644 --- a/mpisppy/cylinders/hub.py +++ b/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") diff --git a/mpisppy/extensions/reduced_costs_fixer.py b/mpisppy/extensions/reduced_costs_fixer.py index a76b7aef4..4226b4595 100644 --- a/mpisppy/extensions/reduced_costs_fixer.py +++ b/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: diff --git a/mpisppy/generic_cylinders.py b/mpisppy/generic_cylinders.py index 89e876d6c..8c6358b60 100644 --- a/mpisppy/generic_cylinders.py +++ b/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.subgradient_hub: + # 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, diff --git a/mpisppy/opt/aph.py b/mpisppy/opt/aph.py index 9539c9c5d..bcfb833e3 100644 --- a/mpisppy/opt/aph.py +++ b/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() diff --git a/mpisppy/opt/ph.py b/mpisppy/opt/ph.py index d08373105..0f6d39f6c 100644 --- a/mpisppy/opt/ph.py +++ b/mpisppy/opt/ph.py @@ -28,13 +28,16 @@ class PH(mpisppy.phbase.PHBase): # uncomment the line below to get per-rank profile outputs, which can # be examined with snakeviz (or your favorite profile output analyzer) #@profile(filename="profile_out") - def ph_main(self, finalize=True): + def ph_main(self, finalize=True, attach_prox=True): """ Execute the PH algorithm. Args: finalize (bool, optional, default=True): If True, call `PH.post_loops()`, if False, do not, and return None for Eobj + attach_prox (bool, optional, default=True): + If True, use a proximal term, if False, no proximal + term is added to the objective function Returns: tuple: @@ -56,13 +59,15 @@ def ph_main(self, finalize=True): """ verbose = self.options['verbose'] smoothed = self.options['smoothed'] - self.PH_Prep(attach_smooth = smoothed) + self.PH_Prep(attach_prox=attach_prox, attach_smooth = smoothed) if (verbose): - print('Calling PH Iter0 on global rank {}'.format(global_rank)) + print(f'Calling {self.__class__.__name__} Iter0 on global rank {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)) + print(f'Completed {self.__class__.__name__} Iter0 on global rank {global_rank}') if ('asynchronousPH' in self.options) and (self.options['asynchronousPH']): raise RuntimeError("asynchronousPH is deprecated; use APH") diff --git a/mpisppy/opt/subgradient.py b/mpisppy/opt/subgradient.py new file mode 100644 index 000000000..ecfdd0e44 --- /dev/null +++ b/mpisppy/opt/subgradient.py @@ -0,0 +1,111 @@ +############################################################################### +# 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.opt.ph +import mpisppy.MPI as _mpi + +_global_rank = _mpi.COMM_WORLD.Get_rank() + +class Subgradient(mpisppy.opt.ph.PH): + """ 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). + """ + return self.ph_main(finalize=finalize) + + def ph_main(self, finalize=True): + # for use with the PH hub + if self.options["smoothed"] != 0: + raise RuntimeError("Cannnot use smoothing with Subgradient algorithm") + return super().ph_main(finalize=finalize, attach_prox=False) + + 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) + + diff --git a/mpisppy/phbase.py b/mpisppy/phbase.py index bdd8a5a04..795ce9d29 100644 --- a/mpisppy/phbase.py +++ b/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") diff --git a/mpisppy/spbase.py b/mpisppy/spbase.py index 92f25a4e7..11479e824 100644 --- a/mpisppy/spbase.py +++ b/mpisppy/spbase.py @@ -98,6 +98,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") diff --git a/mpisppy/spopt.py b/mpisppy/spopt.py index e458183c2..6da2815e9 100644 --- a/mpisppy/spopt.py +++ b/mpisppy/spopt.py @@ -19,6 +19,7 @@ import pyomo.environ as pyo from pyomo.opt import SolverFactory +from pyomo.common.collections import ComponentSet from mpisppy.spbase import SPBase import mpisppy.utils.sputils as sputils @@ -65,6 +66,7 @@ def __init__( # object to get garbage collected to # free the memory the C++ model uses. SPPresolve(self).presolve() + self._create_fixed_nonant_cache() self.current_solver_options = None self.extensions = extensions self.extension_kwargs = extension_kwargs @@ -944,6 +946,20 @@ def _create_solvers(self, presolve=True): print("Set instance times: \tmin=%4.2f mean=%4.2f max=%4.2f" % (np.min(asit), np.mean(asit), np.max(asit))) + def _create_fixed_nonant_cache(self): + self._initial_fixed_varibles = ComponentSet() + for s in self.local_scenarios.values(): + for v in s._mpisppy_data.nonant_indices.values(): + if v.fixed: + self._initial_fixed_varibles.add(v) + + def _can_update_best_bound(self): + for s in self.local_scenarios.values(): + for v in s._mpisppy_data.nonant_indices.values(): + if v.fixed: + if v not in self._initial_fixed_varibles: + return False + return True def subproblem_scenario_generator(self): """ diff --git a/mpisppy/utils/cfg_vanilla.py b/mpisppy/utils/cfg_vanilla.py index 393b2a015..9f412e786 100644 --- a/mpisppy/utils/cfg_vanilla.py +++ b/mpisppy/utils/cfg_vanilla.py @@ -18,6 +18,7 @@ from mpisppy.opt.ph import PH from mpisppy.opt.aph import APH from mpisppy.opt.lshaped import LShapedMethod +from mpisppy.opt.subgradient import Subgradient from mpisppy.fwph.fwph import FWPH from mpisppy.utils.xhat_eval import Xhat_Eval import mpisppy.utils.sputils as sputils @@ -175,6 +176,47 @@ def aph_hub(cfg, return hub_dict +def subgradient_hub(cfg, + scenario_creator, + scenario_denouement, + all_scenario_names, + scenario_creator_kwargs=None, + ph_extensions=None, + extension_kwargs=None, + ph_converger=None, + rho_setter=None, + variable_probability=None, + all_nodenames=None, +): + shoptions = shared_options(cfg) + options = copy.deepcopy(shoptions) + options["convthresh"] = cfg.intra_hub_conv_thresh + options["bundles_per_rank"] = cfg.bundles_per_rank + options["smoothed"] = 0 + + hub_dict = { + "hub_class": PHHub, + "hub_kwargs": {"options": {"rel_gap": cfg.rel_gap, + "abs_gap": cfg.abs_gap, + "max_stalled_iters": cfg.max_stalled_iters}}, + "opt_class": Subgradient, + "opt_kwargs": { + "options": options, + "all_scenario_names": all_scenario_names, + "scenario_creator": scenario_creator, + "scenario_creator_kwargs": scenario_creator_kwargs, + "scenario_denouement": scenario_denouement, + "rho_setter": rho_setter, + "variable_probability": variable_probability, + "extensions": ph_extensions, + "extension_kwargs": extension_kwargs, + "ph_converger": ph_converger, + "all_nodenames": all_nodenames + } + } + add_wxbar_read_write(hub_dict, cfg) + add_ph_tracking(hub_dict, cfg) + return hub_dict def extension_adder(hub_dict,ext_class): # TBD March 2023: this is not really good enough diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index f68bde8f8..73cd117dd 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -420,6 +420,12 @@ def aph_args(self): domain=float, default=0.01) + def subgradient_args(self): + + self.add_to_config(name="subgradient_hub", + description="Use subgradient hub instead of PH (default False)", + domain=bool, + default=False) def fixer_args(self): @@ -630,7 +636,7 @@ def lagranger_args(self): default=None) - def subgradient_args(self): + def subgradient_bounder_args(self): self.add_to_config('subgradient', description="have a subgradient spoke",