Merge pull request #428 from bknueven/sep_rho

SepRho and CoeffRho

doc/src/extensions.rst

@@ -116,6 +116,23 @@ constructor or in the hub dictionary under ``opt_kwargs`` as the
 
 There is an example of the function in the sizes example (``_rho_setter``).
 
+SepRho
+^^^^^^
+
+Set per variable rho values using the "SEP" algorithm from
+
+Progressive hedging innovations for a class of stochastic mixed-integer resource allocation problems
+Jean-Paul Watson, David L. Woodruff, Compu Management Science, 2011
+DOI 10.1007/s10287-010-0125-4
+
+One can additional specify a multiplier on the computed value (default = 1.0).
+If the cost coefficient on a non-anticipative variable is 0, the default rho value is used instead.
+
+CoeffRho
+^^^^^^^^
+
+Set per variable rho values proportional to the cost coefficient on each non-anticipative variable,
+with an optional multiplier (default = 1.0). If the coefficient is 0, the default rho value is used instead.
 
 wtracker_extension
 ^^^^^^^^^^^^^^^^^^

examples/afew.py

@@ -47,7 +47,7 @@ def do_one(dirname, progname, np, argstring):
 # for farmer, the first arg is num_scens and is required
 do_one("farmer", "farmer_cylinders.py", 3,
        "--num-scens=3 --bundles-per-rank=0 --max-iterations=50 "
-       "--default-rho=1 --display-convergence-detail "
+       "--default-rho=1 --sep-rho --display-convergence-detail "
        "--solver-name={} --xhatshuffle --lagrangian --use-norm-rho-updater".format(solver_name))
 do_one("farmer", "farmer_lshapedhub.py", 2,
        "--num-scens=3 --bundles-per-rank=0 --max-iterations=50 "
@@ -58,7 +58,7 @@ def do_one(dirname, progname, np, argstring):
        4,
        "--num-scens=3 --bundles-per-rank=0 --max-iterations=5 "
        "--iter0-mipgap=0.01 --iterk-mipgap=0.001 --linearize-proximal-terms "
-       " --xhatshuffle --lagrangian --fwph "
+       " --xhatshuffle --lagrangian --fwph --smoothing "
        "--default-rho=1 --solver-name={} --display-progress".format(solver_name))
 
 do_one("hydro", "hydro_cylinders_pysp.py", 3,

examples/farmer/farmer_cylinders.py

@@ -45,6 +45,8 @@ def _parse_args():
     cfg.wxbar_read_write_args()
     cfg.tracking_args()
     cfg.reduced_costs_args()
+    cfg.sep_rho_args()
+    cfg.coeff_rho_args()
     cfg.add_to_config("crops_mult",
                          description="There will be 3x this many crops (default 1)",
                          domain=int,
@@ -125,6 +127,11 @@ def main():
     if cfg.reduced_costs:
         vanilla.add_reduced_costs_fixer(hub_dict, cfg)
 
+    if cfg.sep_rho:
+        vanilla.add_sep_rho(hub_dict, cfg)
+    if cfg.coeff_rho:
+        vanilla.add_coeff_rho(hub_dict, cfg)
+
     # FWPH spoke
     if cfg.fwph:
         fw_spoke = vanilla.fwph_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs)
@@ -145,6 +152,10 @@ def main():
         ph_ob_spoke = vanilla.ph_ob_spoke(*beans,
                                           scenario_creator_kwargs=scenario_creator_kwargs,
                                           rho_setter = rho_setter)
+        if cfg.sep_rho:
+            vanilla.add_sep_rho(ph_ob_spoke, cfg)
+        if cfg.coeff_rho:
+            vanilla.add_coeff_rho(ph_ob_spoke, cfg)
 
     # xhat looper bound spoke
     if cfg.xhatlooper:

examples/run_all.py

@@ -286,11 +286,10 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring):
        4,
        "--instance-name=sslp_15_45_10 --bundles-per-rank=2 "
        "--max-iterations=5 --default-rho=1 "
-       "--subgradient --xhatshuffle --fwph "
+       "--subgradient --xhatshuffle --fwph  --coeff-rho "
        "--linearize-proximal-terms "
        "--rel-gap=0.0 "
        "--solver-name={} --fwph-stop-check-tol 0.01".format(solver_name))
-
 do_one("sslp",
        "sslp_cylinders.py",
        3,

examples/sslp/sslp_cylinders.py

@@ -32,6 +32,7 @@ def _parse_args():
     cfg.xhatshuffle_args()
     cfg.subgradient_args()
     cfg.reduced_costs_args()
+    cfg.coeff_rho_args()
     cfg.parse_command_line("sslp_cylinders")
     return cfg
 
@@ -87,9 +88,15 @@ def main():
     if reduced_costs:
         vanilla.add_reduced_costs_fixer(hub_dict, cfg)
 
+    if cfg.coeff_rho:
+        vanilla.add_coeff_rho(hub_dict, cfg)
+
     # FWPH spoke
     if fwph:
         fw_spoke = vanilla.fwph_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs)
+        # Need to fix FWPH to support extensions
+        # if cfg.coeff_rho:
+        #     vanilla.add_coeff_rho(fw_spoke, cfg)
 
     # Standard Lagrangian bound spoke
     if lagrangian:
@@ -100,6 +107,8 @@ def main():
         subgradient_spoke = vanilla.subgradient_spoke(*beans,
                                               scenario_creator_kwargs=scenario_creator_kwargs,
                                               rho_setter = None)
+        if cfg.coeff_rho:
+            vanilla.add_coeff_rho(subgradient_spoke, cfg)
 
     # xhat looper bound spoke
     if xhatlooper:

mpisppy/extensions/coeff_rho.py

@@ -0,0 +1,40 @@
+###############################################################################
+# 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.extensions.extension
+
+from mpisppy.utils.sputils import nonant_cost_coeffs
+
+
+class CoeffRho(mpisppy.extensions.extension.Extension):
+    """
+    Determine rho as a linear function of the objective coefficient
+    """
+
+    def __init__(self, ph):
+        self.ph = ph
+        self.multiplier = 1.0
+        if (
+            "coeff_rho_options" in ph.options
+            and "multiplier" in ph.options["coeff_rho_options"]
+        ):
+            self.multiplier = ph.options["coeff_rho_options"]["multiplier"]
+
+    def post_iter0(self):
+        for s in self.ph.local_scenarios.values():
+            cc = nonant_cost_coeffs(s)
+            for ndn_i, rho in s._mpisppy_model.rho.items():
+                if cc[ndn_i] != 0:
+                    rho._value = abs(cc[ndn_i]) * self.multiplier
+                # if self.ph.cylinder_rank==0:
+                #     nv = s._mpisppy_data.nonant_indices[ndn_i] # var_data object
+                #     print(ndn_i,nv.getname(),cc[ndn_i],rho._value)
+
+        if self.ph.cylinder_rank == 0:
+            print("Rho values updated by CoeffRho Extension")

mpisppy/extensions/reduced_costs_fixer.py

@@ -196,7 +196,7 @@ def reduced_costs_fixing(self, reduced_costs):
             if self.opt.cylinder_rank == 0 and self.verbose:
                 print("All reduced costs are nan, heuristic fixing will not be applied")
             return
-        
+
         # compute the quantile target
         abs_reduced_costs = np.abs(reduced_costs)
 

mpisppy/extensions/sep_rho.py

@@ -0,0 +1,170 @@
+###############################################################################
+# 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.extensions.extension
+import numpy as np
+import mpisppy.MPI as MPI
+
+from mpisppy.utils.sputils import nonant_cost_coeffs
+
+
+class SepRho(mpisppy.extensions.extension.Extension):
+    """
+    Rho determination algorithm "SEP" from
+    Progressive hedging innovations for a class of stochastic mixed-integer
+        resource allocation problems
+    Jean-Paul Watson, David L. Woodruff, Compu Management Science, 2011
+    DOI 10.1007/s10287-010-0125-4
+    """
+
+    def __init__(self, ph):
+        self.ph = ph
+
+        self.multiplier = 1.0
+
+        if (
+            "sep_rho_options" in ph.options
+            and "multiplier" in ph.options["sep_rho_options"]
+        ):
+            self.multiplier = ph.options["sep_rho_options"]["multiplier"]
+
+    def _compute_primal_residual_norm(self, ph):
+        local_nodenames = []
+        local_primal_residuals = {}
+        global_primal_residuals = {}
+
+        for k, s in ph.local_scenarios.items():
+            nlens = s._mpisppy_data.nlens
+            for node in s._mpisppy_node_list:
+                if node.name not in local_nodenames:
+                    ndn = node.name
+                    local_nodenames.append(ndn)
+                    nlen = nlens[ndn]
+
+                    local_primal_residuals[ndn] = np.zeros(nlen, dtype="d")
+                    global_primal_residuals[ndn] = np.zeros(nlen, dtype="d")
+
+        for k, s in ph.local_scenarios.items():
+            nlens = s._mpisppy_data.nlens
+            xbars = s._mpisppy_model.xbars
+            for node in s._mpisppy_node_list:
+                ndn = node.name
+                primal_residuals = local_primal_residuals[ndn]
+
+                unweighted_primal_residuals = np.fromiter(
+                    (
+                        abs(v._value - xbars[ndn, i]._value)
+                        for i, v in enumerate(node.nonant_vardata_list)
+                    ),
+                    dtype="d",
+                    count=nlens[ndn],
+                )
+                primal_residuals += s._mpisppy_probability * unweighted_primal_residuals
+
+        for nodename in local_nodenames:
+            ph.comms[nodename].Allreduce(
+                [local_primal_residuals[nodename], MPI.DOUBLE],
+                [global_primal_residuals[nodename], MPI.DOUBLE],
+                op=MPI.SUM,
+            )
+
+        primal_resid = {}
+        for ndn, global_primal_resid in global_primal_residuals.items():
+            for i, v in enumerate(global_primal_resid):
+                primal_resid[ndn, i] = v
+
+        return primal_resid
+
+    @staticmethod
+    def _compute_min_max(ph, npop, mpiop, start):
+        local_nodenames = []
+        local_xmaxmin = {}
+        global_xmaxmin = {}
+
+        for k, s in ph.local_scenarios.items():
+            nlens = s._mpisppy_data.nlens
+            for node in s._mpisppy_node_list:
+                if node.name not in local_nodenames:
+                    ndn = node.name
+                    local_nodenames.append(ndn)
+                    nlen = nlens[ndn]
+
+                    local_xmaxmin[ndn] = start * np.ones(nlen, dtype="d")
+                    global_xmaxmin[ndn] = np.zeros(nlen, dtype="d")
+
+        for k, s in ph.local_scenarios.items():
+            nlens = s._mpisppy_data.nlens
+            for node in s._mpisppy_node_list:
+                ndn = node.name
+                xmaxmin = local_xmaxmin[ndn]
+
+                xmaxmin_partial = np.fromiter(
+                    (v._value for v in node.nonant_vardata_list),
+                    dtype="d",
+                    count=nlens[ndn],
+                )
+                xmaxmin = npop(xmaxmin, xmaxmin_partial)
+                local_xmaxmin[ndn] = xmaxmin
+
+        for nodename in local_nodenames:
+            ph.comms[nodename].Allreduce(
+                [local_xmaxmin[nodename], MPI.DOUBLE],
+                [global_xmaxmin[nodename], MPI.DOUBLE],
+                op=mpiop,
+            )
+
+        xmaxmin_dict = {}
+        for ndn, global_xmaxmin_dict in global_xmaxmin.items():
+            for i, v in enumerate(global_xmaxmin_dict):
+                xmaxmin_dict[ndn, i] = v
+
+        return xmaxmin_dict
+
+    @staticmethod
+    def _compute_xmax(ph):
+        return SepRho._compute_min_max(ph, np.maximum, MPI.MAX, -np.inf)
+
+    @staticmethod
+    def _compute_xmin(ph):
+        return SepRho._compute_min_max(ph, np.minimum, MPI.MIN, np.inf)
+
+    def pre_iter0(self):
+        pass
+
+    def post_iter0(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)
+            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 miditer(self):
+        pass
+
+    def enditer(self):
+        pass
+
+    def post_everything(self):
+        pass