PHTracking extension

examples/farmer/farmer_cylinders.py

@@ -15,6 +15,7 @@
 from mpisppy.extensions.norm_rho_updater import NormRhoUpdater
 from mpisppy.convergers.norm_rho_converger import NormRhoConverger
 from mpisppy.convergers.primal_dual_converger import PrimalDualConverger
+from mpisppy.utils.cfg_vanilla import extension_adder
 
 write_solution = True
 
@@ -33,6 +34,8 @@ def _parse_args():
     cfg.lagranger_args()
     cfg.xhatshuffle_args()
     cfg.converger_args()
+    cfg.wxbar_read_write_args()
+    cfg.tracking_args()
     cfg.add_to_config("crops_mult",
                          description="There will be 3x this many crops (default 1)",
                          domain=int,
@@ -101,11 +104,12 @@ def main():
         hub_dict['opt_kwargs']['options']\
             ['primal_dual_converger_options'] = {
                 'verbose': True,
-                'tol': cfg.primal_dual_converger_tol}
+                'tol': cfg.primal_dual_converger_tol,
+                'tracking': True}
 
     ## hack in adaptive rho
     if cfg.use_norm_rho_updater:
-        hub_dict['opt_kwargs']['extensions'] = NormRhoUpdater
+        extension_adder(hub_dict, NormRhoUpdater)
         hub_dict['opt_kwargs']['options']['norm_rho_options'] = {'verbose': True}
 
     # FWPH spoke

examples/run_all.py

@@ -158,10 +158,16 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring):
        "--num-scens 3 --bundles-per-rank=0 --max-iterations=50 --default-rho=1 --solver-name={} "
        "--primal-dual-converger --primal-dual-converger-tol=0.5 --lagrangian --xhatshuffle "
        "--intra-hub-conv-thresh -0.1 --rel-gap=1e-6".format(solver_name))
-do_one("farmer", "farmer_cylinders.py",  4,
+do_one("farmer", "farmer_cylinders.py",  5,
        "--num-scens 3 --bundles-per-rank=0 --max-iterations=50 --default-rho=1 --solver-name={} "
-       "--use-norm-rho-converger --use-norm-rho-updater --lagrangian --xhatshuffle --fwph "
-       "--display-convergence-detail".format(solver_name))
+       "--use-norm-rho-converger --use-norm-rho-updater --rel-gap=1e-6 --lagrangian --lagranger "
+       "--xhatshuffle --fwph --W-fname=out_ws.txt --Xbar-fname=out_xbars.txt "
+       "--ph-track-progress --track-convergence=4 --track-xbar=4 --track-nonants=4 "
+       "--track-duals=4".format(solver_name))
+do_one("farmer", "farmer_cylinders.py",  5,
+       "--num-scens 3 --bundles-per-rank=0 --max-iterations=50 --default-rho=1 --solver-name={} "
+       "--use-norm-rho-converger --use-norm-rho-updater --lagrangian --lagranger --xhatshuffle --fwph "
+       "--init-W-fname=out_ws.txt --init-Xbar-fname=out_xbars.txt --ph-track-progress --track-convergence=4 "  "--track-xbar=4 --track-nonants=4 --track-duals=4 ".format(solver_name))
 do_one("farmer", "farmer_lshapedhub.py", 2,
        "--num-scens 3 --bundles-per-rank=0 --max-iterations=50 "
        "--solver-name={} --rel-gap=0.0 "
@@ -361,14 +367,14 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring):
            "--default-rho=1 --num-scens=3 --max-solver-threads=2 "
            "--lagrangian-iter0-mipgap=1e-7 --lagrangian --xhatshuffle "
            "--ph-mipgaps-json=phmipgaps.json "
-           "--solver-name={}".format(solver_name))    
+           "--solver-name={}".format(solver_name))
     # as of May 2022, this one works well, but outputs some crazy messages
     do_one("uc", "uc_ama.py", 3,
            "--bundles-per-rank=0 --max-iterations=2 "
            "--default-rho=1 --num-scens=3 "
            "--fixer-tol=1e-2 --lagranger --xhatshuffle "
            "--solver-name={}".format(solver_name))
-    
+
     # 10-scenario UC
     time_one("UC_cylinder10scen", "uc", "uc_cylinders.py", 3,
              "--bundles-per-rank=5 --max-iterations=2 "

mpisppy/convergers/converger.py

@@ -27,9 +27,15 @@ def __init__(self, opt):
     @abc.abstractmethod
     def is_converged(self):
         ''' Indicated whether the algorithm has converged.
-            
+
             Must return a boolean. If True, the algorithm will terminate at the
             current iteration--no more solves will be performed by SPBase.
             Otherwise, the iterations will continue.
         '''
         pass
+
+    def post_loops(self):
+        '''Method called after the termination of the algorithm.
+            This method is called after the post_loops of any extensions
+        '''
+        pass

mpisppy/convergers/primal_dual_converger.py

@@ -1,35 +1,51 @@
-import math
 import numpy as np
+import os
+import pandas as pd
 import mpisppy.convergers.converger
+import matplotlib.pyplot as plt
 from mpisppy import MPI
+from mpisppy.extensions.phtracker import TrackedData
 
 class PrimalDualConverger(mpisppy.convergers.converger.Converger):
     """ Convergence checker for the primal-dual metrics.
         Primal convergence is measured as weighted sum over all scenarios s
         p_{s} * ||x_{s} - \bar{x}||_1.
-        Dual convergence is measured as 
+        Dual convergence is measured as
         rho * ||\bar{x}_{t} - \bar{x}_{t-1}||_1
     """
     def __init__(self, ph):
         """ Initialization method for the PrimalDualConverger class."""
         super().__init__(ph)
-        if 'primal_dual_converger_options' in ph.options and \
-                'verbose' in ph.options['primal_dual_converger_options'] and \
-                ph.options['primal_dual_converger_options']['verbose']:
-            self._verbose = True
-        else:
-            self._verbose = False
-        self._ph = ph
-        self.convergence_threshold = ph.options['primal_dual_converger_options']\
-            ['tol']
 
+        self.options = ph.options.get('primal_dual_converger_options', {})
+        self._verbose = self.options.get('verbose', False)
+        self._ph = ph
+        self.convergence_threshold = self.options.get('tol', 1)
+        self.tracking = self.options.get('tracking', False)
         self.prev_xbars = self._get_xbars()
+        self._rank = self._ph.cylinder_rank
+
+        if self.tracking and self._rank == 0:
+            # if phtracker is set up, save the results in the phtracker/hub folder
+            if 'phtracker_options' in self._ph.options:
+                tracker_options = self._ph.options["phtracker_options"]
+                cylinder_name = tracker_options.get(
+                    "cylinder_name", type(self._ph.spcomm).__name__)
+                results_folder = tracker_options.get(
+                    "results_folder", "results")
+                results_folder = os.path.join(results_folder, cylinder_name)
+            else:
+                results_folder = self.options.get('results_folder', 'results')
+            self.tracker = TrackedData('pd', results_folder, plot=True, verbose=self._verbose)
+            os.makedirs(results_folder, exist_ok=True)
+            self.tracker.initialize_fnames(name=self.options.get('pd_fname', None))
+            self.tracker.initialize_df(['iteration', 'primal_gap', 'dual_gap'])
 
     def _get_xbars(self):
         """
         Get the current xbar values from the local scenarios
         Returns:
-            xbars (dict): dictionary of xbar values indexed by 
+            xbars (dict): dictionary of xbar values indexed by
                           (decision node name, index)
         """
         xbars = {}
@@ -38,7 +54,7 @@ def _get_xbars(self):
                 xbars[ndn_i] = xbar.value
             break
         return xbars
-    
+
     def _compute_primal_convergence(self):
         """
         Compute the primal convergence metric
@@ -55,12 +71,12 @@ def _compute_primal_convergence(self):
                 nlen = s._mpisppy_data.nlens[ndn]
                 x_bars = np.fromiter((s._mpisppy_model.xbars[ndn,i]._value
                                       for i in range(nlen)), dtype='d')
-                
+
                 nonants_array = np.fromiter(
                     (v._value for v in node.nonant_vardata_list),
                     dtype='d', count=nlen)
                 _l1 = np.abs(x_bars - nonants_array)
-                
+
                 # invariant to prob_coeff being a scalar or array
                 prob = s._mpisppy_data.prob_coeff[ndn] * np.ones(nlen)
                 local_sum_diff[0] += np.dot(prob, _l1)
@@ -107,13 +123,40 @@ def is_converged(self):
         self.prev_xbars = self._get_xbars()
         ret_val = max(primal_gap, dual_gap) <= self.convergence_threshold
 
-        if self._verbose and self._ph.cylinder_rank == 0:
+        if self._verbose and self._rank == 0:
             print(f"primal gap = {round(primal_gap, 5)}, dual gap = {round(dual_gap, 5)}")
 
             if ret_val:
                 print("Dual convergence check passed")
             else:
                 print("Dual convergence check failed "
                       f"(requires primal + dual gaps) <= {self.convergence_threshold}")
-
+        if self.tracking and self._rank == 0:
+            self.tracker.add_row([self._ph._PHIter, primal_gap, dual_gap])
+            self.tracker.write_out_data()
         return ret_val
+
+    def plot_results(self):
+        """
+        Plot the results of the convergence checks
+        by reading in csv file and plotting
+        """
+        plot_fname = self.tracker.plot_fname
+        conv_data = pd.read_csv(self.tracker.fname)
+
+        # Create a log-scale plot
+        plt.semilogy(conv_data['iteration'], conv_data['primal_gap'], label='Primal Gap')
+        plt.semilogy(conv_data['iteration'], conv_data['dual_gap'], label='Dual Gap')
+
+        plt.xlabel('Iteration')
+        plt.ylabel('Convergence Metric')
+        plt.legend()
+        plt.savefig(plot_fname)
+        plt.close()
+
+    def post_everything(self):
+        '''
+        Reading the convergence data and plotting the results
+        '''
+        if self.tracking and self._rank == 0:
+            self.plot_results()

mpisppy/cylinders/lagranger_bounder.py

@@ -93,13 +93,13 @@ def main(self):
 
         self.lagrangian_prep()
 
+        self.A_iter = 1
         self.trivial_bound = self._lagrangian(0)
 
         self.bound = self.trivial_bound
 
         self.opt.current_solver_options = self.opt.iterk_solver_options
 
-        self.A_iter = 1
         while not self.got_kill_signal():
             # because of aph, do not check for new data, just go for it
             self.bound = self._update_weights_and_solve(self.A_iter)
@@ -113,4 +113,7 @@ def finalize(self):
         '''
         self.final_bound = self._update_weights_and_solve(self.A_iter)
         self.bound = self.final_bound
+        if self.opt.extensions is not None and \
+            hasattr(self.opt.extobject, 'post_everything'):
+            self.opt.extobject.post_everything()
         return self.final_bound

mpisppy/cylinders/lagrangian_bounder.py

@@ -67,19 +67,19 @@ def main(self):
 
         self.lagrangian_prep()
 
+        self.dk_iter = 1
         self.trivial_bound = self.lagrangian()
 
         self.opt.current_solver_options = self.opt.iterk_solver_options
 
         self.bound = self.trivial_bound
 
-        dk_iter = 1
         while not self.got_kill_signal():
             if self.new_Ws:
                 bound = self._set_weights_and_solve()
                 if bound is not None:
                     self.bound = bound
-            dk_iter += 1
+                self.dk_iter += 1
 
     def finalize(self):
         '''
@@ -89,4 +89,7 @@ def finalize(self):
         '''
         self.final_bound = self._set_weights_and_solve()
         self.bound = self.final_bound
+        if self.opt.extensions is not None and \
+            hasattr(self.opt.extobject, 'post_everything'):
+            self.opt.extobject.post_everything()
         return self.final_bound

mpisppy/extensions/extension.py

@@ -44,6 +44,18 @@ def post_solve(self, subproblem, results):
         '''
         return results
 
+    def pre_solve_loop(self):
+        ''' Method called before every solve loop within
+            mpisppy.spot.SPOpt.solve_loop()
+        '''
+        pass
+
+    def post_solve_loop(self):
+        ''' Method called after every solve loop within
+            mpisppy.spot.SPOpt.solve_loop()
+        '''
+        pass
+
     def pre_iter0(self):
         ''' Method called at the end of PH_Prep().
             When this method is called, all scenarios have been created, and
@@ -62,7 +74,7 @@ def post_iter0(self):
 
     def miditer(self):
         ''' Method called after x-bar has been computed and the dual weights
-            have been updated, but before solve_loop(). 
+            have been updated, but before solve_loop().
             If a converger is present, this method is called between the
             convergence_value() method and the is_converged() method.
         '''
@@ -97,10 +109,27 @@ def __init__(self, ph, ext_classes):
             name = constr.__name__
             self.extdict[name] = constr(ph)
 
+    def pre_solve(self, subproblem):
+        for lobject in self.extdict.values():
+            lobject.pre_solve(subproblem)
+
+    def post_solve(self, subproblem, results):
+        for lobject in self.extdict.values():
+            results = lobject.post_solve(subproblem, results)
+        return results
+
+    def pre_solve_loop(self):
+        for lobject in self.extdict.values():
+            lobject.pre_solve_loop()
+
+    def post_solve_loop(self):
+        for lobject in self.extdict.values():
+            lobject.post_solve_loop()
+
     def pre_iter0(self):
         for lobject in self.extdict.values():
             lobject.pre_iter0()
-                                        
+
     def post_iter0(self):
         for lobject in self.extdict.values():
             lobject.post_iter0()