Dynsep

examples/farmer/farmer_cylinders.py

@@ -47,6 +47,7 @@ def _parse_args():
     cfg.reduced_costs_args()
     cfg.sep_rho_args()
     cfg.coeff_rho_args()
+    cfg.dynamic_rho_args()
     cfg.add_to_config("crops_mult",
                          description="There will be 3x this many crops (default 1)",
                          domain=int,

examples/generic_cylinders.bash

@@ -4,8 +4,18 @@
 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
 
+echo "^^^ sep rho not 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 --sep-rho --rel-gap 0.001
+
+
+echo "^^^ sep 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 --sep-rho --rel-gap 0.001 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1
+
+exit
+
 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
@@ -16,7 +26,6 @@ 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

mpisppy/extensions/dyn_rho_base.py

@@ -55,7 +55,7 @@ def _display_W_values(self):
     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]
+            curr_conv, last_conv = self.primal_conv_cache[0], self.primal_conv_cache[-1]
             try:
                 primal_diff =  np.abs((last_conv - curr_conv) / last_conv)
             except Warning:
@@ -65,14 +65,16 @@ def _update_rho_primal_based(self):
             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]
+        if len(self.dual_conv_cache) < 4:
+            # first two entries will be 0 by construction, so ignore
+            return False
+        curr_conv, last_conv = self.dual_conv_cache[0], self.dual_conv_cache[-1]
         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())
+        return (hasattr(self.cfg, "dynamic_rho_primal_crit") and self.cfg.dynamic_rho_primal_crit and self._update_rho_primal_based()) or \
+               (hasattr(self.cfg, "dynamic_rho_dual_crit") and self.cfg.dynamic_rho_dual_crit and self._update_rho_dual_based())
 
     def pre_iter0(self):
         pass

mpisppy/extensions/sep_rho.py

@@ -7,14 +7,14 @@
 # full copyright and license information.
 ###############################################################################
 
-import mpisppy.extensions.extension
+import mpisppy.extensions.dyn_rho_base
 import numpy as np
 import mpisppy.MPI as MPI
-
+from mpisppy import global_toc
 from mpisppy.utils.sputils import nonant_cost_coeffs
 
 
-class SepRho(mpisppy.extensions.extension.Extension):
+class SepRho(mpisppy.extensions.dyn_rho_base.Dyn_Rho_extension_base):
     """
     Rho determination algorithm "SEP" from
     Progressive hedging innovations for a class of stochastic mixed-integer
@@ -23,7 +23,8 @@ class SepRho(mpisppy.extensions.extension.Extension):
     DOI 10.1007/s10287-010-0125-4
     """
 
-    def __init__(self, ph):
+    def __init__(self, ph, comm=None):
+        super().__init__(ph, comm=comm)
         self.ph = ph
 
         self.multiplier = 1.0
@@ -33,6 +34,9 @@ def __init__(self, ph):
             and "multiplier" in ph.options["sep_rho_options"]
         ):
             self.multiplier = ph.options["sep_rho_options"]["multiplier"]
+        self.cfg = ph.options["sep_rho_options"]["cfg"]
+
+        self._nonant_cost_coeffs = None
 
     def _compute_primal_residual_norm(self, ph):
         local_nodenames = []
@@ -134,34 +138,55 @@ def _compute_xmax(ph):
     def _compute_xmin(ph):
         return SepRho._compute_min_max(ph, np.minimum, MPI.MIN, np.inf)
 
-    def pre_iter0(self):
-        pass
+    def nonant_cost_coeffs(self, s):
+        if self._nonant_cost_coeffs is None:
+            # Should be the same in every scenario...
+            self._nonant_cost_coeffs = nonant_cost_coeffs(s)
+        return self._nonant_cost_coeffs
 
-    def post_iter0(self):
+    def _compute_and_update_rho(self):
         ph = self.ph
         primal_resid = self._compute_primal_residual_norm(ph)
         xmax = self._compute_xmax(ph)
         xmin = self._compute_xmin(ph)
 
         for s in ph.local_scenarios.values():
-            cc = nonant_cost_coeffs(s)
+            cc = self.nonant_cost_coeffs(s)
             for ndn_i, rho in s._mpisppy_model.rho.items():
                 if cc[ndn_i] != 0:
                     nv = s._mpisppy_data.nonant_indices[ndn_i]  # var_data object
                     if nv.is_integer():
                         rho._value = abs(cc[ndn_i]) / (xmax[ndn_i] - xmin[ndn_i] + 1)
                     else:
                         rho._value = abs(cc[ndn_i]) / max(1, primal_resid[ndn_i])
-
                     rho._value *= self.multiplier
 
-                # if ph.cylinder_rank==0:
-                #     print(ndn_i,nv.getname(),xmax[ndn_i],xmin[ndn_i],primal_resid[ndn_i],cc[ndn_i],rho._value)
-        if ph.cylinder_rank == 0:
-            print("Rho values updated by SepRho Extension")
+    def compute_and_update_rho(self):
+        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 pre_iter0(self):
+        pass
+
+    def post_iter0(self):
+        global_toc("Using sep-rho rho setter")
+        super().post_iter0()        
+        self.compute_and_update_rho()
 
     def miditer(self):
-        pass
+        # this code could be factored; see sens_rho.py
+        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()
 
     def enditer(self):
         pass

mpisppy/utils/cfg_vanilla.py

@@ -204,7 +204,7 @@ def add_fixer(hub_dict,
 
 def add_sep_rho(hub_dict, cfg):
     hub_dict = extension_adder(hub_dict,SepRho)
-    hub_dict["opt_kwargs"]["options"]["sep_rho_options"] = {"multiplier" : cfg.sep_rho_multiplier}
+    hub_dict["opt_kwargs"]["options"]["sep_rho_options"] = {"multiplier" : cfg.sep_rho_multiplier, "cfg": cfg}
 
 def add_coeff_rho(hub_dict, cfg):
     hub_dict = extension_adder(hub_dict,CoeffRho)

mpisppy/utils/config.py

@@ -150,9 +150,9 @@ def _bad_rho_setters(msg):
 
         if self.grad_rho_setter and self.sensi_rho:
             _bad_rho_setters("Only one rho setter can be active.")
-        if not self.grad_rho_setter and not self.sensi_rho:
+        if not self.grad_rho_setter and not self.sensi_rho and not self.sep_rho:
             if self.dynamic_rho_primal_crit or self.dynamic_rho_dual_crit:
-                _bad_rho_setters("dynamic rho only works with grad- and sensi-")
+                _bad_rho_setters("dynamic rho only works with grad-, sensi-, and sep-rho")
 
     def add_solver_specs(self, prefix=""):
         sstr = f"{prefix}_solver" if prefix != "" else "solver"
@@ -845,12 +845,12 @@ def dynamic_rho_args(self): # AKA adaptive
         self.gradient_args()
 
         self.add_to_config('dynamic_rho_primal_crit',
-                           description="Use dynamic primal criterion for some rho updates",
+                           description="Use dynamic primal criterion for some types of rho updates",
                            domain=bool,
                            default=False)
 
         self.add_to_config('dynamic_rho_dual_crit',
-                           description="Use dynamic dual criterion for some rho updates",
+                           description="Use dynamic dual criterion for some stypes of rho updates",
                            domain=bool,
                            default=False)