From f4cc7e62271cd79466c4a4f83cfb2f7ef9f74d44 Mon Sep 17 00:00:00 2001 From: David L Woodruff Date: Tue, 1 Oct 2024 10:16:42 -0700 Subject: [PATCH] Dyn rho (#435) * dynamic code is factored; gradient works; sensi is next * sensi rho and grad rho share dynamics * ruff * change import in test * add a little test for dynamic gradient-based rho --- examples/generic_cylinders.bash | 18 +++- mpisppy/extensions/dyn_rho_base.py | 113 +++++++++++++++++++++++ mpisppy/extensions/gradient_extension.py | 54 +---------- mpisppy/extensions/sensi_rho.py | 28 ++++-- mpisppy/generic_cylinders.py | 2 +- mpisppy/tests/test_gradient_rho.py | 15 ++- mpisppy/utils/cfg_vanilla.py | 2 +- mpisppy/utils/config.py | 18 ++-- 8 files changed, 179 insertions(+), 71 deletions(-) create mode 100644 mpisppy/extensions/dyn_rho_base.py diff --git a/examples/generic_cylinders.bash b/examples/generic_cylinders.bash index d999ce3bc..61492b8d4 100644 --- a/examples/generic_cylinders.bash +++ b/examples/generic_cylinders.bash @@ -3,6 +3,21 @@ SOLVER="cplex" +cd farmer +mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01 --solution-base-name farmer_nonants + +mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --bundles-per-rank=0 --max-iterations=100 --default-rho=1 --solver-name=${SOLVER} --xhatpath=./farmer_nonants.npy --grad-order-stat 0.0 --xhatxbar --ph-ob --max-stalled-iters 5000 --grad-rho-setter --rel-gap 0.001 + +# now do it again, but this time using dynamic rho +echo "^^^ grad rho dynamic ^^^" +mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --bundles-per-rank=0 --max-iterations=100 --default-rho=1 --solver-name=${SOLVER} --xhatpath=./farmer_nonants.npy --grad-order-stat 0.0 --xhatxbar --ph-ob --max-stalled-iters 5000 --grad-rho-setter --rel-gap 0.001 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1 + +echo "^^^ sensi rho dynamic ^^^" +mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --bundles-per-rank=0 --max-iterations=100 --default-rho=1 --solver-name=${SOLVER} --xhatpath=./farmer_nonants.npy --grad-order-stat 0.0 --xhatxbar --ph-ob --max-stalled-iters 5000 --sensi-rho --rel-gap 0.001 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1 + +cd .. +exit + echo "^^^ netdes sensi-rho ^^^" cd netdes mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name netdes --netdes-data-path ./data --instance-name network-10-20-L-01 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01 --sensi-rho @@ -11,7 +26,6 @@ cd .. echo "^^^ farmer sensi-rho ^^^" mpiexec -np 3 python -m mpi4py ../mpisppy/generic_cylinders.py --module-name farmer/farmer --num-scens 3 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01 --sensi-rho -exit # sslp EF echo "^^^ sslp ef ^^^" @@ -94,7 +108,7 @@ mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name # now do it again, but this time using dynamic rho -mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --bundles-per-rank=0 --max-iterations=100 --default-rho=1 --solver-name=${SOLVER} --xhatpath=./farmer_nonants.npy --grad-order-stat 0.0 --xhatxbar --ph-ob --max-stalled-iters 5000 --grad-rho-setter --rel-gap 0.001 --grad-dynamic-dual-crit --grad-dynamic-dual-thresh 0.1 +mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --bundles-per-rank=0 --max-iterations=100 --default-rho=1 --solver-name=${SOLVER} --xhatpath=./farmer_nonants.npy --grad-order-stat 0.0 --xhatxbar --ph-ob --max-stalled-iters 5000 --grad-rho-setter --rel-gap 0.001 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1 cd .. # end gradient based rho demo diff --git a/mpisppy/extensions/dyn_rho_base.py b/mpisppy/extensions/dyn_rho_base.py new file mode 100644 index 000000000..8a0ab842e --- /dev/null +++ b/mpisppy/extensions/dyn_rho_base.py @@ -0,0 +1,113 @@ +############################################################################### +# 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. +############################################################################### + +# A dynamic rho base class that assumes a member function compute_and_update_rho +# As of Oct 2024, it only mainly provides services and needs a better mid_iter +# Children should call the parent for post_iter0 + + +import numpy as np + +import mpisppy.extensions.extension +from mpisppy.utils.wtracker import WTracker +from mpisppy import global_toc + +# for trapping numpy warnings +import warnings + +class Dyn_Rho_extension_base(mpisppy.extensions.extension.Extension): + """ + This extension enables dynamic rho, e.g. gradient or sensitivity based. + + Args: + opt (PHBase object): gives the problem + cfg (Config object): config object + + """ + def __init__(self, opt, comm=None): + super().__init__(opt) + self.cylinder_rank = self.opt.cylinder_rank + + self.primal_conv_cache = [] + self.dual_conv_cache = [] + self.wt = WTracker(self.opt) + + def _display_rho_values(self): + for sname, scenario in self.opt.local_scenarios.items(): + rho_list = [scenario._mpisppy_model.rho[ndn_i]._value + for ndn_i, _ in scenario._mpisppy_data.nonant_indices.items()] + print(sname, 'rho values: ', rho_list[:5]) + break + + def _display_W_values(self): + for (sname, scenario) in self.opt.local_scenarios.items(): + W_list = [w._value for w in scenario._mpisppy_model.W.values()] + print(sname, 'W values: ', W_list) + break + + + def _update_rho_primal_based(self): + with warnings.catch_warnings(): + warnings.filterwarnings('error') + curr_conv, last_conv = self.primal_conv_cache[-1], self.primal_conv_cache[-2] + try: + primal_diff = np.abs((last_conv - curr_conv) / last_conv) + except Warning: + if self.cylinder_rank == 0: + print(f"dyn_rho extension reports {last_conv=} {curr_conv=} - no rho updates recommended") + return False + return (primal_diff <= self.cfg.dynamic_rho_primal_thresh) + + def _update_rho_dual_based(self): + curr_conv, last_conv = self.dual_conv_cache[-1], self.dual_conv_cache[-2] + dual_diff = np.abs((last_conv - curr_conv) / last_conv) if last_conv != 0 else 0 + #print(f'{dual_diff =}') + return (dual_diff <= self.cfg.dynamic_rho_dual_thresh) + + def _update_recommended(self): + return (self.cfg.dynamic_rho_primal_crit and self._update_rho_primal_based()) or \ + (self.cfg.dynamic_rho_dual_crit and self._update_rho_dual_based()) + + def pre_iter0(self): + pass + + def post_iter0(self): + self.primal_conv_cache.append(self.opt.convergence_diff()) + self.dual_conv_cache.append(self.wt.W_diff()) + + def miditer(self): + self.primal_conv_cache.append(self.opt.convergence_diff()) + self.dual_conv_cache.append(self.wt.W_diff()) + if self.opt._PHIter == 1: + self.grad_object.write_grad_cost() + if self.opt._PHIter == 1 or self._update_recommended(): + self.grad_object.write_grad_rho() + rho_setter_kwargs = self.opt.options['rho_setter_kwargs'] \ + if 'rho_setter_kwargs' in self.opt.options \ + else dict() + + sum_rho = 0.0 + num_rhos = 0 + for sname, scenario in self.opt.local_scenarios.items(): + rholist = self.rho_setter(scenario, **rho_setter_kwargs) + for (vid, rho) in rholist: + (ndn, i) = scenario._mpisppy_data.varid_to_nonant_index[vid] + scenario._mpisppy_model.rho[(ndn, i)] = rho + sum_rho += rho + num_rhos += 1 + + rho_avg = sum_rho / num_rhos + + global_toc(f"Rho values recomputed - average rank 0 rho={rho_avg}") + + def enditer(self): + pass + + def post_everything(self): + pass diff --git a/mpisppy/extensions/gradient_extension.py b/mpisppy/extensions/gradient_extension.py index 236e9b016..5a45628d4 100644 --- a/mpisppy/extensions/gradient_extension.py +++ b/mpisppy/extensions/gradient_extension.py @@ -8,18 +8,14 @@ ############################################################################### import os -import numpy as np -import mpisppy.extensions.extension +import mpisppy.extensions.dyn_rho_base import mpisppy.utils.gradient as grad import mpisppy.utils.find_rho as find_rho -from mpisppy.utils.wtracker import WTracker from mpisppy import global_toc -# for trapping numpy warnings -import warnings -class Gradient_extension(mpisppy.extensions.extension.Extension): +class Gradient_extension(mpisppy.extensions.dyn_rho_base.Dyn_Rho_extension_base): """ This extension makes PH use gradient-rho and the corresponding rho setter. @@ -32,8 +28,8 @@ class Gradient_extension(mpisppy.extensions.extension.Extension): """ def __init__(self, opt, comm=None): - super().__init__(opt) - self.cylinder_rank = self.opt.cylinder_rank + super().__init__(opt, comm=comm) + self.cfg = opt.options["gradient_extension_options"]["cfg"] # This is messy because we want to be able to use or give rhos as requested. # (e.g., if the user gave us an input file, use that) @@ -60,53 +56,13 @@ def __init__(self, opt, comm=None): self.grad_object = grad.Find_Grad(opt, self.cfg) self.rho_setter = find_rho.Set_Rho(self.cfg).rho_setter - self.primal_conv_cache = [] - self.dual_conv_cache = [] - self.wt = WTracker(self.opt) - - def _display_rho_values(self): - for sname, scenario in self.opt.local_scenarios.items(): - rho_list = [scenario._mpisppy_model.rho[ndn_i]._value - for ndn_i, _ in scenario._mpisppy_data.nonant_indices.items()] - print(sname, 'rho values: ', rho_list[:5]) - break - - def _display_W_values(self): - for (sname, scenario) in self.opt.local_scenarios.items(): - W_list = [w._value for w in scenario._mpisppy_model.W.values()] - print(sname, 'W values: ', W_list) - break - - - def _update_rho_primal_based(self): - with warnings.catch_warnings(): - warnings.filterwarnings('error') - curr_conv, last_conv = self.primal_conv_cache[-1], self.primal_conv_cache[-2] - try: - primal_diff = np.abs((last_conv - curr_conv) / last_conv) - except Warning: - if self.cylinder_rank == 0: - print(f"Gradient extension reports {last_conv=} {curr_conv=} - no rho updates recommended") - return False - return (primal_diff <= self.cfg.grad_dynamic_primal_thresh) - - def _update_rho_dual_based(self): - curr_conv, last_conv = self.dual_conv_cache[-1], self.dual_conv_cache[-2] - dual_diff = np.abs((last_conv - curr_conv) / last_conv) if last_conv != 0 else 0 - #print(f'{dual_diff =}') - return (dual_diff <= self.cfg.grad_dynamic_dual_thresh) - - def _update_recommended(self): - return (self.cfg.grad_dynamic_primal_crit and self._update_rho_primal_based()) or \ - (self.cfg.grad_dynamic_dual_crit and self._update_rho_dual_based()) def pre_iter0(self): pass def post_iter0(self): global_toc("Using gradient-based rho setter") - self.primal_conv_cache.append(self.opt.convergence_diff()) - self.dual_conv_cache.append(self.wt.W_diff()) + super().post_iter0() def miditer(self): self.primal_conv_cache.append(self.opt.convergence_diff()) diff --git a/mpisppy/extensions/sensi_rho.py b/mpisppy/extensions/sensi_rho.py index 9a63fec1c..57366109e 100644 --- a/mpisppy/extensions/sensi_rho.py +++ b/mpisppy/extensions/sensi_rho.py @@ -8,19 +8,19 @@ ############################################################################### import numpy as np - - -import mpisppy.extensions.extension +from mpisppy import global_toc +import mpisppy.extensions.dyn_rho_base import mpisppy.MPI as MPI from mpisppy.utils.nonant_sensitivities import nonant_sensitivies -class SensiRho(mpisppy.extensions.extension.Extension): +class SensiRho(mpisppy.extensions.dyn_rho_base.Dyn_Rho_extension_base): """ Rho determination algorithm using nonant sensitivities """ - def __init__(self, ph): + def __init__(self, ph, comm=None): + super().__init__(ph, comm=comm) self.ph = ph self.multiplier = 1.0 @@ -30,6 +30,7 @@ def __init__(self, ph): and "multiplier" in ph.options["sensi_rho_options"] ): self.multiplier = ph.options["sensi_rho_options"]["multiplier"] + self.cfg = ph.options["sensi_rho_options"]["cfg"] @staticmethod def _compute_rho_min_max(ph, npop, mpiop, start): @@ -134,7 +135,6 @@ def _compute_rho_max(ph): def _compute_rho_min(ph): return SensiRho._compute_rho_min_max(ph, np.minimum, MPI.MIN, np.inf) - # Implementation of a virtual function def compute_and_update_rho(self): ph = self.ph @@ -167,10 +167,24 @@ def pre_iter0(self): pass def post_iter0(self): + global_toc("Using sensi-rho rho setter") + super().post_iter0() self.compute_and_update_rho() def miditer(self): - pass + self.primal_conv_cache.append(self.opt.convergence_diff()) + self.dual_conv_cache.append(self.wt.W_diff()) + + if self._update_recommended(): + self.compute_and_update_rho() + sum_rho = 0.0 + num_rhos = 0 # could be computed... + for sname, s in self.opt.local_scenarios.items(): + for ndn_i, nonant in s._mpisppy_data.nonant_indices.items(): + sum_rho += s._mpisppy_model.rho[ndn_i]._value + num_rhos += 1 + rho_avg = sum_rho / num_rhos + global_toc(f"Rho values recomputed - average rank 0 rho={rho_avg}") def enditer(self): pass diff --git a/mpisppy/generic_cylinders.py b/mpisppy/generic_cylinders.py index 70ab324cd..ad18286a7 100644 --- a/mpisppy/generic_cylinders.py +++ b/mpisppy/generic_cylinders.py @@ -68,7 +68,7 @@ def _parse_args(m): cfg.wxbar_read_write_args() cfg.tracking_args() cfg.gradient_args() - cfg.dynamic_gradient_args() + cfg.dynamic_rho_args() cfg.reduced_costs_args() cfg.sep_rho_args() cfg.coeff_rho_args() diff --git a/mpisppy/tests/test_gradient_rho.py b/mpisppy/tests/test_gradient_rho.py index 19af19d16..f775b36fe 100644 --- a/mpisppy/tests/test_gradient_rho.py +++ b/mpisppy/tests/test_gradient_rho.py @@ -39,7 +39,7 @@ def _create_cfg(): cfg.popular_args() cfg.two_sided_args() cfg.ph_args() - cfg.dynamic_gradient_args() + cfg.dynamic_rho_args() cfg.solver_name = solver_name cfg.default_rho = 1 return cfg @@ -124,7 +124,6 @@ def test_grad_cost_and_rho(self): self.assertAlmostEqual(float(rows[3][1]), 2.0616161616161617) os.remove(self.cfg.grad_cost_file_out) os.remove(self.cfg.grad_rho_file_out) - ############################################### @@ -268,6 +267,18 @@ def test_grad_extensions(self): self.ph_object = self._run_ph_farmer() self.assertAlmostEqual(self.ph_object.conv, 2.128968965420187, places=1) + def test_dyn_grad_extensions(self): + print("** test dynamic grad extensions **") + self.cfg.grad_rho_file_out = './_test_rho.csv' + #self.cfg.grad_cost_file_in = './examples/rho_test_data/grad_cost.csv' + self.cfg.xhatpath = './examples/rho_test_data/farmer_cyl_nonants.npy' + self.cfg.max_iterations = 4 + self.cfg.dynamic_rho_dual_crit = True + self.cfg.dynamic_rho_dual_thresh = 0.1 + self.ph_object = self._run_ph_farmer() + print(f"{self.ph_object.conv=}") + self.assertAlmostEqual(self.ph_object.conv, 0.039598, places=1) + if __name__ == '__main__': unittest.main() diff --git a/mpisppy/utils/cfg_vanilla.py b/mpisppy/utils/cfg_vanilla.py index 59ba3d868..ffdd9d1f1 100644 --- a/mpisppy/utils/cfg_vanilla.py +++ b/mpisppy/utils/cfg_vanilla.py @@ -212,7 +212,7 @@ def add_coeff_rho(hub_dict, cfg): def add_sensi_rho(hub_dict, cfg): hub_dict = extension_adder(hub_dict,SensiRho) - hub_dict["opt_kwargs"]["options"]["sensi_rho_options"] = {"multiplier" : cfg.sensi_rho_multiplier} + hub_dict["opt_kwargs"]["options"]["sensi_rho_options"] = {"multiplier" : cfg.sensi_rho_multiplier, "cfg": cfg} def add_cross_scenario_cuts(hub_dict, cfg, diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index 0e72b1f19..a8b334a2a 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -826,27 +826,27 @@ def gradient_args(self): domain=float, default=1e3) - def dynamic_gradient_args(self): # AKA adaptive + def dynamic_rho_args(self): # AKA adaptive self.gradient_args() - self.add_to_config('grad_dynamic_primal_crit', - description="Use dynamic gradient-based primal criterion for update", + self.add_to_config('dynamic_rho_primal_crit', + description="Use dynamic primal criterion for some rho updates", domain=bool, default=False) - self.add_to_config('grad_dynamic_dual_crit', - description="Use dynamic gradient-based dual criterion for update", + self.add_to_config('dynamic_rho_dual_crit', + description="Use dynamic dual criterion for some rho updates", domain=bool, default=False) - self.add_to_config("grad_dynamic_primal_thresh", - description="primal threshold for diff during gradient calcs", + self.add_to_config("dynamic_rho_primal_thresh", + description="primal threshold for diff during dynamic rho calcs", domain=float, default=0.1) - self.add_to_config("grad_dynamic_dual_thresh", - description="dual threshold for abs norm during gradient calcs", + self.add_to_config("dynamic_rho_dual_thresh", + description="dual threshold for dirr during dynamic rho calcs", domain=float, default=0.1)