Lp writer (#470)

* start a special extension to write lp files, etc. for each scenario

* fix a potential bug in the use of proper bundles

* the lp and nonant writer extension seems to be working

* deal with some proper bundle issues raised by, and for, the sizes example

* add a test and a little doc for the scenario lp writer extension

* the scenario_lpfiles extension

* remove some notes

* use a Pyomo function to convert variable names per Ben's suggestion

doc/src/extensions.rst

@@ -15,6 +15,13 @@ extension can be included in some applications, but not others.
 There are a number of extensions, particularly for PH, that are provided
 with ``mpi-sppy`` and they provide examples that can be used for the
 creation of more. Extensions can be found in ``mpisppy.extensions``.
+Note that some things (e.g. some xhatters) can be used as a cylinder
+or as an extension. A few other things (e.g., cross scenario cuts) need
+both an extension and a cylinder.
+
+Many extensions are supported in :ref:`generic_cylinders`. The rest of
+this help file describes extensions released with mpisppy along with
+some hints for including them in your own cylinders driver program.
 
 Multiple Extensions
 -------------------
@@ -266,3 +273,12 @@ If some variables have zero probability in all scenarios, then you will need to
 ``do_not_check_variable_probabilities`` to True in the options for ``spbase``. This will result in skipping the checks for
 all variable probabilities! So you might want to set this to False to verify that the probabilities sum to one
 only for the Vars you expect before setting it to True.
+
+Scenario_lpwriter
+-----------------
+
+This extension writes an lp file with the model and json file with (a) list(s) of
+scenario tree node names and nonanticaptive variables for each scenario before
+the iteration zero solve of PH or APH. Note that for two-stage problems, all
+json files will be the same. See ``mpisppy.generic_cylinders.py``
+for an example of use.

examples/generic_cylinders.bash

@@ -4,6 +4,25 @@
 SOLVER="cplex"
 SPB=1
 
+echo "^^^ write scenario lp and nonant json files ^^^"
+cd sizes
+python ../../mpisppy/generic_cylinders.py --module-name sizes --num-scens 3 --default-rho 1 --solver-name ${SOLVER} --max-iterations 0 --scenario-lpfiles
+cd ..
+
+echo "^^^ pickle sizes bundles ^^^"
+cd sizes
+python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name sizes --num-scens 10 --pickle-bundles-dir sizes_pickles --scenarios-per-bundle 5 --default-rho 1
+cd ..
+
+echo "^^^ unpickle the sizes bundles and write the lp and nonant files ^^^"
+# note that numscens need to match the number before pickling...
+# so does scenarios per bundle
+cd sizes
+python ../../mpisppy/generic_cylinders.py --module-name sizes --num-scens 10 --default-rho 1 --solver-name ${SOLVER} --max-iterations 0 --scenario-lpfiles --unpickle-bundles-dir sizes_pickles --scenarios-per-bundle 5
+cd ..
+echo "xxxx Early exit. xxxx"
+exit
+
 echo "^^^ pickle the scenarios ^^^"
 cd farmer
 python ../../mpisppy/generic_cylinders.py --module-name farmer --pickle-scenarios-dir farmer_pickles --crops-mult 2 --num-scens 10
@@ -53,8 +72,6 @@ mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name
 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

examples/sizes/sizes.py

@@ -39,11 +39,11 @@ def scenario_denouement(rank, scenario_name, scenario):
 
 #=========
 def scenario_names_creator(num_scens,start=None):
-    # (only for Amalgamator): return the full list of num_scens scenario names
     # if start!=None, the list starts with the 'start' labeled scenario
+    # note that the scenarios for the sizes problem are one-based
     if (start is None) :
-        start=0
-    return [f"Scenario{i+1}" for i in range(start,start+num_scens)]
+        start=1
+    return [f"Scenario{i}" for i in range(start, start+num_scens)]
 
 
 #=========
@@ -97,8 +97,13 @@ def _rho_setter(scen, **kwargs):
         scen (pyo.ConcreteModel): the scenario
     Returns:
         a list of (id(vardata), rho)
+    Note:
+        This rho_setter will not work with proper bundles.
     """
     retlist = []
+    if not hasattr(scen, "UnitReductionCost"):
+        print("WARNING: _rho_setter not used (probably because of proper bundles)")
+        return retlist        
     RF = 0.001  # a factor for rho, if you like
 
     if "RF" in kwargs and isinstance(kwargs["RF"], float):

mpisppy/extensions/scenario_lpfiles.py

@@ -0,0 +1,42 @@
+###############################################################################
+# 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.
+###############################################################################
+""" Extension to write an lp file for each scenario and a json file for
+the nonant structure for each scenario (yes, for two stage problems this
+json file will be the same for all scenarios.)
+"""
+
+import json
+import mpisppy.extensions.extension
+import pyomo.core.base.label as pyomo_label
+
+
+def lpize(varname):
+    # convert varname to the string that will appear in the lp file
+    # return varname.replace("[", "(").replace("]", ")").replace(",", "_").replace(".","_")
+    return pyomo_label.cpxlp_label_from_name(varname)
+
+
+class Scenario_lpfiles(mpisppy.extensions.extension.Extension):
+
+    def __init__(self, ph):
+        self.ph = ph
+
+    def pre_iter0(self):
+        for k, s in self.ph.local_subproblems.items():
+            s.write(f"{k}.lp", io_options={'symbolic_solver_labels': True})
+            nonants_by_node = {nd.name: [lpize(var.name) for var in nd.nonant_vardata_list] for nd in s._mpisppy_node_list}
+            with open(f"{k}_nonants.json", "w") as jfile:
+                json.dump(nonants_by_node, jfile)
+
+    def post_iter0(self):
+        return
+
+
+
+

mpisppy/generic_cylinders.py

@@ -24,6 +24,7 @@
 from mpisppy.extensions.fixer import Fixer
 from mpisppy.extensions.mipgapper import Gapper
 from mpisppy.extensions.gradient_extension import Gradient_extension
+from mpisppy.extensions.scenario_lpfiles import Scenario_lpfiles
 import mpisppy.utils.solver_spec as solver_spec
 from mpisppy import global_toc
 from mpisppy import MPI
@@ -45,6 +46,11 @@ def _parse_args(m):
                       description="The string used for a directory of ouput along with a csv and an npv file (default None, which means no soltion output)",
                       domain=str,
                       default=None)
+    cfg.add_to_config(name="scenario_lpfiles",
+                      description="Invokes an extension that writes an model lp file and a nonants json file for each scenario before iteration 0",
+                      domain=bool,
+                      default=False)
+
     m.inparser_adder(cfg)
     # many models, e.g., farmer, need num_scens_required
     #  in which case, it should go in the inparser_adder function
@@ -80,7 +86,7 @@ def _parse_args(m):
     cfg.checker()  # looks for inconsistencies 
     return cfg
 
-def _name_lists(module, cfg):
+def _name_lists(module, cfg, bundle_wrapper=None):
 
     # Note: high level code like this assumes there are branching factors for
     # multi-stage problems. For other trees, you will need lower-level code
@@ -99,15 +105,16 @@ def _name_lists(module, cfg):
     # proper bundles should be almost magic
     if cfg.unpickle_bundles_dir or cfg.scenarios_per_bundle is not None:
         num_buns = cfg.num_scens // cfg.scenarios_per_bundle
-        all_scenario_names = wrapper.bundle_names_creator(num_buns, cfg=cfg)
+        all_scenario_names = bundle_wrapper.bundle_names_creator(num_buns, cfg=cfg)
+        all_nodenames = None  # This is seldom used; also, proper bundles result in two stages
     else:
         all_scenario_names = module.scenario_names_creator(num_scens)
 
     return all_scenario_names, all_nodenames
 
 
 #==========
-def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement):
+def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement, bundle_wrapper=None):
     rho_setter = module._rho_setter if hasattr(module, '_rho_setter') else None
     if cfg.default_rho is None and rho_setter is None:
         if cfg.sep_rho or cfg.coeff_rho or cfg.sensi_rho:
@@ -125,7 +132,7 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_
     else:
         ph_converger = None
 
-    all_scenario_names, all_nodenames = _name_lists(module, cfg)    
+    all_scenario_names, all_nodenames = _name_lists(module, cfg, bundle_wrapper=bundle_wrapper)    
 
     # Things needed for vanilla cylinders
     beans = (cfg, scenario_creator, scenario_denouement, all_scenario_names)
@@ -179,6 +186,9 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_
         ext_classes.append(Gradient_extension)
         hub_dict['opt_kwargs']['options']['gradient_extension_options'] = {'cfg': cfg}        
 
+    if cfg.scenario_lpfiles:
+        ext_classes.append(Scenario_lpfiles)
+
     if cfg.sep_rho:
         vanilla.add_sep_rho(hub_dict, cfg)
 
@@ -397,9 +407,9 @@ def _write_bundles(module,
     global_toc(f"Bundles written to {cfg.pickle_bundles_dir}")
 
 #==========
-def _do_EF(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement):
+def _do_EF(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement, bundle_wrapper=None):
 
-    all_scenario_names, _ = _name_lists(module, cfg)
+    all_scenario_names, _ = _name_lists(module, cfg, bundle_wrapper=bundle_wrapper)
     ef = sputils.create_EF(
         all_scenario_names,
         module.scenario_creator,
@@ -481,16 +491,17 @@ def _proper_bundles(cfg):
 
     cfg = _parse_args(module)
 
+    bundle_wrapper = None  # the default
     if _proper_bundles(cfg):
         # TBD: remove the need for dill if you are not reading or writing
         import mpisppy.utils.pickle_bundle as pickle_bundle
         import mpisppy.utils.proper_bundler as proper_bundler
 
-        wrapper = proper_bundler.ProperBundler(module)
-        scenario_creator = wrapper.scenario_creator
+        bundle_wrapper = proper_bundler.ProperBundler(module)
+        scenario_creator = bundle_wrapper.scenario_creator
         # The scenario creator is wrapped, so these kw_args will not go the original
         # creator (the kw_creator will keep the original args)
-        scenario_creator_kwargs = wrapper.kw_creator(cfg)
+        scenario_creator_kwargs = bundle_wrapper.kw_creator(cfg)
     elif cfg.unpickle_scenarios_dir is not None:
         # So reading pickled scenarios cannot be composed with proper bundles
         import mpisppy.utils.pickle_bundle as pickle_bundle
@@ -520,6 +531,6 @@ def _proper_bundles(cfg):
                          scenario_denouement,
                          global_comm)
     elif cfg.EF:
-        _do_EF(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement)
+        _do_EF(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement, bundle_wrapper=bundle_wrapper)
     else:
-        _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement)
+        _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement, bundle_wrapper=bundle_wrapper)

mpisppy/tests/test_ef_ph.py

@@ -15,6 +15,7 @@
 
 import os
 import glob
+import json
 import unittest
 import pandas as pd
 import pyomo.environ as pyo
@@ -379,6 +380,45 @@ def test_xhat_extension(self):
         xhatobj1 = round_pos_sig(ph.extobject._xhat_looper_obj_final, 1)
         self.assertEqual(xhatobj1, 200000)
 
+
+    @unittest.skipIf(not pyo.SolverFactory('glpk').available(),
+                     "glpk is not available")
+    def test_scenario_lpwriter_extension(self):
+        print("test scenarip_lpwriter")
+        from mpisppy.extensions.scenario_lpfiles import Scenario_lpfiles
+        options = self._copy_of_base_options()
+        options["iter0_solver_options"] = {"mipgap": 0.1}    
+        options["PHIterLimit"] = 0
+        options["solver_name"] = "glpk"
+        options["tee_rank0_solves"] = True
+
+        ph = mpisppy.opt.ph.PH(
+            options,
+            self.all3_scenario_names,
+            scenario_creator,
+            scenario_denouement,
+            scenario_creator_kwargs={"scenario_count": 3},
+            extensions=Scenario_lpfiles,
+        )
+        conv, basic_obj, tbound = ph.ph_main()
+        # The idea is to detect a change in Pyomo's writing of lp files
+        with open("Scenario1_nonants.json", "r") as jfile:
+            nonants_by_node = json.load(jfile)
+        vname = nonants_by_node["ROOT"][0]  # first name in the file
+        gotit = False
+        with open("Scenario1.lp", 'r') as lpfile:
+            for line in lpfile:
+                if vname in line:
+                    gotit = True
+                    break
+        assert gotit, f"The first nonant in Scenario1_nonants.json ({vname}) not found in Scenario1.lp"
+        print("   deleting Scenario*.p and Scenario*_nonants.json")        
+        for fn in glob.glob("Scenario*.lp"):
+            os.remove(fn)
+        for fn in glob.glob("Scenario*_nonants.json"):
+            os.remove(fn)
+
+
     @unittest.skipIf(not solver_available,
                      "no solver is available")
     def test_wtracker_extension(self):

mpisppy/utils/proper_bundler.py

@@ -94,7 +94,6 @@ def scenario_creator(self, sname, **kwargs):
             # snames are scenario names
             snames = self.module.scenario_names_creator(lastnum-firstnum+1,
                                                         firstnum)
-
             # We are assuming seeds are managed by the *scenario* creator.
             bundle = sputils.create_EF(snames, self.module.scenario_creator,
                                        scenario_creator_kwargs=self.original_kwargs,