trying to pass LB/UB from RC Spoke

Author: bknueven

mpisppy/cylinders/hub.py

@@ -12,22 +12,21 @@
 import mpisppy.log
 
 from mpisppy.cylinders.spcommunicator import RecvArray, SPCommunicator
-from math import inf
 
 from mpisppy import global_toc
 
 from mpisppy.cylinders.spwindow import Field
 
 # Could also pass, e.g., sys.stdout instead of a filename
-mpisppy.log.setup_logger("mpisppy.cylinders.Hub",
+mpisppy.log.setup_logger(__name__,
                          "hub.log",
                          level=logging.CRITICAL)
-logger = logging.getLogger("mpisppy.cylinders.Hub")
+logger = logging.getLogger(__name__)
 
 class Hub(SPCommunicator):
 
     send_fields = (*SPCommunicator.send_fields, Field.SHUTDOWN, Field.BEST_OBJECTIVE_BOUNDS,)
-    receive_fields = (*SPCommunicator.receive_fields, Field.OBJECTIVE_INNER_BOUND, Field.OBJECTIVE_OUTER_BOUND, )
+    receive_fields = (*SPCommunicator.receive_fields,)
 
     _hub_algo_best_bound_provider = False
 
@@ -37,16 +36,10 @@ def __init__(self, spbase_object, fullcomm, strata_comm, cylinder_comm, communic
         logger.debug(f"Built the hub object on global rank {fullcomm.Get_rank()}")
         # for logging
         self.print_init = True
-        self.latest_ib_char = None
-        self.latest_ob_char = None
-        self.last_ib_idx = None
-        self.last_ob_idx = None
         # for termination based on stalling out
         self.stalled_iter_cnt = 0
         self.last_gap = float('inf')  # abs_gap tracker
 
-        self.initialize_bound_values()
-
         return
 
     @abc.abstractmethod
@@ -169,66 +162,6 @@ def hub_finalize(self):
             global_toc("Statistics at termination", True)
             self.screen_trace()
 
-    def update_innerbounds(self):
-        """ Update the inner bounds after receiving them from the spokes
-        """
-        logging.debug("Hub is trying to update from InnerBounds")
-        for idx, cls, recv_buf in self.receive_field_spcomms[Field.OBJECTIVE_INNER_BOUND]:
-            if recv_buf.is_new():
-                bound = recv_buf[0]
-                logging.debug("!! new InnerBound to opt {}".format(bound))
-                self.BestInnerBound = self.InnerBoundUpdate(bound, cls, idx)
-        logging.debug("ph back from InnerBounds")
-
-    def update_outerbounds(self):
-        """ Update the outer bounds after receiving them from the spokes
-        """
-        logging.debug("Hub is trying to update from OuterBounds")
-        for idx, cls, recv_buf in self.receive_field_spcomms[Field.OBJECTIVE_OUTER_BOUND]:
-            if recv_buf.is_new():
-                bound = recv_buf[0]
-                logging.debug("!! new OuterBound to opt {}".format(bound))
-                self.BestOuterBound = self.OuterBoundUpdate(bound, cls, idx)
-        logging.debug("ph back from OuterBounds")
-
-    def OuterBoundUpdate(self, new_bound, cls=None, idx=None, char='*'):
-        current_bound = self.BestOuterBound
-        if self._outer_bound_update(new_bound, current_bound):
-            if cls is None:
-                self.latest_ob_char = char
-                self.last_ob_idx = 0
-            else:
-                self.latest_ob_char = cls.converger_spoke_char
-                self.last_ob_idx = idx
-            return new_bound
-        else:
-            return current_bound
-
-    def InnerBoundUpdate(self, new_bound, cls=None, idx=None, char='*'):
-        current_bound = self.BestInnerBound
-        if self._inner_bound_update(new_bound, current_bound):
-            if cls is None:
-                self.latest_ib_char = char
-                self.last_ib_idx = 0
-            else:
-                self.latest_ib_char = cls.converger_spoke_char
-                self.last_ib_idx = idx
-            return new_bound
-        else:
-            return current_bound
-
-    def initialize_bound_values(self):
-        if self.opt.is_minimizing:
-            self.BestInnerBound = inf
-            self.BestOuterBound = -inf
-            self._inner_bound_update = lambda new, old : (new < old)
-            self._outer_bound_update = lambda new, old : (new > old)
-        else:
-            self.BestInnerBound = -inf
-            self.BestOuterBound = inf
-            self._inner_bound_update = lambda new, old : (new > old)
-            self._outer_bound_update = lambda new, old : (new < old)
-
     def _populate_boundsout_cache(self, buf):
         """ Populate a given buffer with the current bounds
         """

mpisppy/cylinders/reduced_costs_spoke.py

@@ -6,6 +6,7 @@
 # All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
 # full copyright and license information.
 ###############################################################################
+import math
 import pyomo.environ as pyo
 import numpy as np
 from mpisppy.cylinders.lagrangian_bounder import LagrangianOuterBound
@@ -15,7 +16,8 @@
 
 class ReducedCostsSpoke(LagrangianOuterBound):
 
-    send_fields = (*LagrangianOuterBound.send_fields, Field.EXPECTED_REDUCED_COST, Field.SCENARIO_REDUCED_COST,)
+    send_fields = (*LagrangianOuterBound.send_fields, Field.EXPECTED_REDUCED_COST, Field.SCENARIO_REDUCED_COST,
+                   Field.NONANT_LOWER_BOUNDS, Field.NONANT_UPPER_BOUNDS,)
     receive_fields = (*LagrangianOuterBound.receive_fields,)
 
     converger_spoke_char = 'R'
@@ -57,8 +59,40 @@ def register_send_fields(self) -> None:
             scenario_buffer_len += len(s._mpisppy_data.nonant_indices)
         self._scenario_rc_buffer = np.zeros(scenario_buffer_len)
 
+        self.initialize_bound_fields()
+        self.create_integer_variable_where()
+
         return
 
+    def initialize_bound_fields(self):
+        self._nonant_lower_bounds = self.send_buffers[Field.NONANT_LOWER_BOUNDS].value_array()
+        self._nonant_upper_bounds = self.send_buffers[Field.NONANT_UPPER_BOUNDS].value_array()
+
+        self._nonant_lower_bounds[:] = -np.inf
+        self._nonant_upper_bounds[:] = np.inf
+
+        for s in self.opt.local_scenarios.values():
+            scenario_lower_bounds = np.fromiter(
+                    _lb_generator(s._mpisppy_data.nonant_indices.values()),
+                    dtype=float,
+                    count=len(s._mpisppy_data.nonant_indices),
+                )
+            self._nonant_lower_bounds = np.maximum(self._nonant_lower_bounds, scenario_lower_bounds)
+
+            scenario_upper_bounds = np.fromiter(
+                    _ub_generator(s._mpisppy_data.nonant_indices.values()),
+                    dtype=float,
+                    count=len(s._mpisppy_data.nonant_indices),
+                )
+            self._nonant_upper_bounds = np.minimum(self._nonant_upper_bounds, scenario_upper_bounds)
+
+    def create_integer_variable_where(self):
+        self._integer_variable_where = np.full(len(self._nonant_lower_bounds), False)
+        for s in self.opt.local_scenarios.values():
+            for idx, xvar in enumerate(s._mpisppy_data.nonant_indices.values()):
+                if xvar.is_integer():
+                    self._integer_variable_where[idx] = True
+
     @property
     def rc_global(self):
         return self.send_buffers[Field.EXPECTED_REDUCED_COST].value_array()
@@ -84,24 +118,19 @@ def lagrangian_prep(self):
         same as base class, but relax the integer variables and
         attach the reduced cost suffix
         """
-        # Split up PH_Prep? Prox option is important for APH.
-        # Seems like we shouldn't need the Lagrangian stuff, so attach_prox=False
-        # Scenarios are created here
-        self.opt.PH_Prep(attach_prox=False)
-        self.opt._reenable_W()
-
         relax_integer_vars = pyo.TransformationFactory("core.relax_integer_vars")
         for s in self.opt.local_subproblems.values():
             relax_integer_vars.apply_to(s)
             s.rc = pyo.Suffix(direction=pyo.Suffix.IMPORT)
-        self.opt._create_solvers(presolve=False)
+        super().lagrangian_prep()
 
     def lagrangian(self, need_solution=True):
         if not need_solution:
             raise RuntimeError("ReducedCostsSpoke always needs a solution to work")
         bound = super().lagrangian(need_solution=need_solution)
         if bound is not None:
             self.extract_and_store_reduced_costs()
+            self.extract_and_store_updated_nonant_bounds(bound)
         return bound
 
     def extract_and_store_reduced_costs(self):
@@ -176,6 +205,56 @@ def extract_and_store_reduced_costs(self):
             Field.SCENARIO_REDUCED_COST,
         )
 
+    def extract_and_store_updated_nonant_bounds(self, lr_outer_bound):
+        # update the best bound from the hub
+        # as a side effect, calls update_receive_buffers
+        # TODO: fix this side effect 
+        if self.got_kill_signal():
+            return
+
+        self.update_innerbounds()        
+
+        if math.isinf(self.BestInnerBound):
+            # can do anything with no bound
+            return
+        
+        nonzero_rc = np.where(self.rc_global==0, np.nan, self.rc_global)
+        bound_tightening = np.divide(self.BestInnerBound - lr_outer_bound, nonzero_rc)
+
+        tighten_upper = np.where(bound_tightening>0, np.nan, bound_tightening)
+        tighten_lower = np.where(bound_tightening<0, np.nan, bound_tightening)
+
+        tighten_upper += self._nonant_lower_bounds 
+        tighten_lower += self._nonant_upper_bounds
+
+        # max of existing lower and new lower
+        np.fmax(tighten_lower, self._nonant_lower_bounds, out=self._nonant_lower_bounds)
+
+        # min of existing upper and new upper
+        np.fmin(tighten_upper, self._nonant_upper_bounds, out=self._nonant_upper_bounds)
+
+        # ceiling of lower bounds for integer variables
+        np.ceil(self._nonant_lower_bounds, out=self._nonant_lower_bounds, where=self._integer_variable_where)
+        # floor of upper bounds for integer variables
+        np.floor(self._nonant_upper_bounds, out=self._nonant_upper_bounds, where=self._integer_variable_where)
+         
+
     def main(self):
         # need the solution for ReducedCostsSpoke
         super().main(need_solution=True)
+
+
+def _lb_generator(var_iterable):
+    for v in var_iterable:
+        lb = v.lb
+        if lb is None:
+            yield -np.inf
+        yield lb
+
+
+def _ub_generator(var_iterable):
+    for v in var_iterable:
+        ub = v.ub
+        if ub is None:
+            yield np.inf
+        yield ub

mpisppy/cylinders/spcommunicator.py

@@ -19,13 +19,17 @@
     Separate hub and spoke classes for memory/window management?
 """
 
-import numpy as np
 import abc
 import time
+import logging
+import numpy as np
+from math import inf
 
-from mpisppy import MPI
+from mpisppy import MPI, global_toc
 from mpisppy.cylinders.spwindow import Field, FieldLengths, SPWindow
 
+logger = logging.getLogger(__name__)
+
 def communicator_array(size):
     arr = np.empty(size+1, dtype='d')
     arr[:] = np.nan
@@ -116,7 +120,8 @@ class SPCommunicator:
         or expects to receive from another SPCommunicator object.
     """
     send_fields = ()
-    receive_fields = (Field.NONANT_LOWER_BOUNDS, Field.NONANT_UPPER_BOUNDS,)
+    receive_fields = (Field.OBJECTIVE_INNER_BOUND, Field.OBJECTIVE_OUTER_BOUND,
+                      Field.NONANT_LOWER_BOUNDS, Field.NONANT_UPPER_BOUNDS,)
 
     def __init__(self, spbase_object, fullcomm, strata_comm, cylinder_comm, communicators, options=None):
         self.fullcomm = fullcomm
@@ -152,6 +157,13 @@ def __init__(self, spbase_object, fullcomm, strata_comm, cylinder_comm, communic
         # the SPBase object
         self.opt.spcomm = self
 
+        # for communicating with bounders
+        self.latest_ib_char = None
+        self.latest_ob_char = None
+        self.last_ib_idx = None
+        self.last_ob_idx = None
+        self.initialize_bound_values()
+
         return
 
     def _make_key(self, field: Field, origin: int):
@@ -385,20 +397,88 @@ def update_nonant_bounds(self):
         upper bound buffers should be up-to-date, which can be done by calling
         `SPCommunicator.update_receive_buffers`.
         """
-        _INF = float("inf")
+        bounds_modified = 0
         for _, _, recv_buf in self.receive_field_spcomms[Field.NONANT_LOWER_BOUNDS]:
-            for s in self.opt.local_scenarios.items():
+            for s in self.opt.local_scenarios.values():
                 for ci, (ndn_i, xvar) in enumerate(s._mpisppy_data.nonant_indices.items()):
                     xvarlb = xvar.lb
                     if xvarlb is None:
-                        xvarlb = -_INF
+                        xvarlb = -inf
                     if recv_buf[ci] > xvarlb:
                         xvar.lb = recv_buf[ci]
+                        bounds_modified += 1
         for _, _, recv_buf in self.receive_field_spcomms[Field.NONANT_UPPER_BOUNDS]:
-            for s in self.opt.local_scenarios.items():
+            for s in self.opt.local_scenarios.values():
                 for ci, (ndn_i, xvar) in enumerate(s._mpisppy_data.nonant_indices.items()):
                     xvarub = xvar.ub
                     if xvarub is None:
-                        xvarub = _INF
+                        xvarub = inf 
                     if recv_buf[ci] < xvarub:
                         xvar.ub = recv_buf[ci]
+                        bounds_modified += 1
+
+        bounds_modified /= len(self.opt.local_scenarios)
+
+        if bounds_modified > 0:
+            global_toc(f"{self.__class__.__name__}: tightened {int(bounds_modified)} variable bounds", self.cylinder_rank == 0)
+
+    def update_innerbounds(self):
+        """ Update the inner bounds after receiving them from the spokes
+        """
+        logger.debug(f"{self.__class__.__name__} is trying to update from InnerBounds")
+        for idx, cls, recv_buf in self.receive_field_spcomms[Field.OBJECTIVE_INNER_BOUND]:
+            if recv_buf.is_new():
+                bound = recv_buf[0]
+                logger.debug("!! new InnerBound to opt {}".format(bound))
+                self.BestInnerBound = self.InnerBoundUpdate(bound, cls, idx)
+        logger.debug(f"{self.__class__.__name__} back from InnerBounds")
+
+    def update_outerbounds(self):
+        """ Update the outer bounds after receiving them from the spokes
+        """
+        logger.debug(f"{self.__class__.__name__} is trying to update from OuterBounds")
+        for idx, cls, recv_buf in self.receive_field_spcomms[Field.OBJECTIVE_OUTER_BOUND]:
+            if recv_buf.is_new():
+                bound = recv_buf[0]
+                logger.debug("!! new OuterBound to opt {}".format(bound))
+                self.BestOuterBound = self.OuterBoundUpdate(bound, cls, idx)
+        logger.debug(f"{self.__class__.__name__} back from OuterBounds")
+
+    def OuterBoundUpdate(self, new_bound, cls=None, idx=None, char='*'):
+        current_bound = self.BestOuterBound
+        if self._outer_bound_update(new_bound, current_bound):
+            if cls is None:
+                self.latest_ob_char = char
+                self.last_ob_idx = 0
+            else:
+                self.latest_ob_char = cls.converger_spoke_char
+                self.last_ob_idx = idx
+            return new_bound
+        else:
+            return current_bound
+
+    def InnerBoundUpdate(self, new_bound, cls=None, idx=None, char='*'):
+        current_bound = self.BestInnerBound
+        if self._inner_bound_update(new_bound, current_bound):
+            if cls is None:
+                self.latest_ib_char = char
+                self.last_ib_idx = 0
+            else:
+                self.latest_ib_char = cls.converger_spoke_char
+                self.last_ib_idx = idx
+            return new_bound
+        else:
+            return current_bound
+
+    def initialize_bound_values(self):
+        if self.opt.is_minimizing:
+            self.BestInnerBound = inf
+            self.BestOuterBound = -inf
+            self._inner_bound_update = lambda new, old : (new < old)
+            self._outer_bound_update = lambda new, old : (new > old)
+        else:
+            self.BestInnerBound = -inf
+            self.BestOuterBound = inf
+            self._inner_bound_update = lambda new, old : (new > old)
+            self._outer_bound_update = lambda new, old : (new < old)
+