diff --git a/doc/src/amalgamator.rst b/doc/src/amalgamator.rst index 92e944877..7adacc052 100644 --- a/doc/src/amalgamator.rst +++ b/doc/src/amalgamator.rst @@ -3,7 +3,11 @@ Amalgamator =========== -For simple problems that do not need extra specification, ``amalgamator.py`` +.. Note:: + This is an advanced topic, probably of interest only to maintainers of + confidence interfal and related code. + +The file ``amalgamator.py`` provides several drivers to solve a problem without writing a program that creates the cylinders one by one. ``amalgamator.from_module`` enables a high-level user to create a hub-and-spoke architecture using the command @@ -49,7 +53,9 @@ It must contains the following attributes for use with cylinders: Create Amalgamator from a module and command line -------------------------------------------------- Given an options +------------------------------------------------- + +Given an options ``Config`` object (typically called `cfg`) as above, ``amalgamator.Amalgamator_parser`` creates calls the appropriate functions to add the necessary information for different modules. @@ -88,7 +94,7 @@ Examples -------- As intended, the examples of use of Amalgamator are quite short. The first -example, ``farmer_ama.py``, solves directly the EF. The model can be specified, +example, ``examples.farmer.archive.farmer_ama.py``, solves directly the EF. The model can be specified, e.g. by taking an integer version of it. This specification can be made via the command line, thanks to the ``inparser_adder`` method of ``farmer.py``. diff --git a/doc/src/confidence_intervals.rst b/doc/src/confidence_intervals.rst index 37b7b5161..a2705e6c6 100644 --- a/doc/src/confidence_intervals.rst +++ b/doc/src/confidence_intervals.rst @@ -57,7 +57,7 @@ This object has a ``run`` method that returns a gap estimator and a confidence i Examples -------- -There are example scripts for sequential sampling in both ``farmer`` and ``aircond``. +There are example scripts for sequential sampling in both ``farmer/CI`` and ``aircond``. Using stand-alone ``mmw_conf.py`` --------------------------------- @@ -79,7 +79,7 @@ to be able to pass problem-specific args down without knowing what they are. Once a model satisfies the requirement for amalgamator, next a ``.npy`` file should be constructed from the given model. This can be accomplished, for example, by adding the line ``sputils.ef_ROOT_nonants_npy_serializer(instance, 'xhat.npy')`` after solving the ef ``instance``. When using ``Amalgamator`` to solve the program, this can be done by adding the line -``sputils.ef_ROOT_nonants_npy_serializer(ama_object.ef, "xhat.npy")`` to your existing program (see the example in ``farmer.py`` for an example of this). +``sputils.ef_ROOT_nonants_npy_serializer(ama_object.ef, "xhat.npy")`` to your existing program (see the example in ``examples/farmer/archive/farmer.py`` for an example of this). Once this is accomplished, on the command line, run ``python -m mpisppy.confidence_intervals.mmw_conf my_model.py xhat.npy gurobi --num-scens n --alpha 0.95``. Note that ``xhat.npy`` is assumed to be in the same directory as ``my_model.py`` in this case. If the file is saved elsewhere then the corresponding path should be called on the command line. diff --git a/doc/src/examples.rst b/doc/src/examples.rst index bdaa76199..2bd7cd3c2 100644 --- a/doc/src/examples.rst +++ b/doc/src/examples.rst @@ -47,7 +47,7 @@ constraints are the number of tons per acre that each crop will yield (2.5 for wheat, 3 for corn, and 20 for sugar beets). -The following code creates an instance of the farmer's model: +The following code in ``examples/farmer/archive/farmer.py`` (with similar code in ``examples/farmer/farmer.py``) creates an instance of the farmer's model: .. testcode:: @@ -109,7 +109,8 @@ yields for each crop. We can solve the model: The optimal objective value is: .. testoutput:: - + :options: +SKIP + -118600.0 In practice, the farmer does not know the number of tons that each crop will @@ -206,136 +207,39 @@ farmer's stochastic program. Solving the Extensive Form ^^^^^^^^^^^^^^^^^^^^^^^^^^ -The simplest approach is to solve the extensive form of the model directly. -MPI-SPPy makes this quite simple: - -.. testcode:: +The simplest approach is to solve the extensive form of the model directly. Assuming you are in the directory +``examples/farmer`` the following unix command will work. - from mpisppy.opt.ef import ExtensiveForm +.. code-block:: bash - options = {"solver": "cplex_direct"} - all_scenario_names = ["good", "average", "bad"] - ef = ExtensiveForm(options, all_scenario_names, scenario_creator) - results = ef.solve_extensive_form() + python ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --EF --EF-solver-name gurobi - objval = ef.get_objective_value() - print(f"{objval:.1f}") +We can extract the optimal solution itself using the ``--solution-base-name`` option: +.. code-block:: bash -.. testoutput:: - - ... - -108390.0 + python ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --EF --EF-solver-name gurobi --solution-base-name farmersol -We can extract the optimal solution itself using the ``get_root_solution`` -method of the ``ExtensiveForm`` object: - -.. testcode:: - - soln = ef.get_root_solution() - for (var_name, var_val) in soln.items(): - print(var_name, var_val) - -.. testoutput:: - - X[BEETS] 250.0 - X[CORN] 80.0 - X[WHEAT] 170.0 +This command writes solution data for nonanticipative variables to two files with the base name farmersol and full scenario solutions to a directory named farmersol_soldir. +.. note:: + Most command line options relevant to the EF start with --EF. Most other command line options will be silently ignored + if ``--EF`` is specified (one exception is ``--solution-base-name``). + Solving Using Progressive Hedging (PH) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -We can also solve the model using the progressive hedging (PH) algorithm. -First, we must construct a PH object: - -.. testcode:: - - from mpisppy.opt.ph import PH - - options = { - "solver_name": "cplex_persistent", - "PHIterLimit": 5, - "defaultPHrho": 10, - "convthresh": 1e-7, - "verbose": False, - "display_progress": False, - "display_timing": False, - "iter0_solver_options": dict(), - "iterk_solver_options": dict(), - } - all_scenario_names = ["good", "average", "bad"] - ph = PH( - options, - all_scenario_names, - scenario_creator, - ) - - -.. testoutput:: - :hide: - - ... - -Note that all of the options in the ``options`` dict must be specified in order -to construct the PH object. Once the PH object is constructed, we can execute -the algorithm with a call to the ``ph_main`` method: - -.. testcode:: - - ph.ph_main() - -.. testoutput:: - :hide: - - ... - - -.. testoutput:: - :options: +SKIP - - - [ 0.00] Start SPBase.__init__ - [ 0.01] Start PHBase.__init__ - [ 0.01] Creating solvers - [ 0.01] Entering solve loop in PHBase.Iter0 - [ 2.80] Reached user-specified limit=5 on number of PH iterations +Here is a simple command that uses PH as the hub algorithm and +computes lower bounds using a Lagrangian spoke (``--lagrangian``) with +upper bounds computed by randomly trying scenario solutions to fix the nonanticipative variables (``--xhatshuffle``). -Note that precise timing results may differ. In this toy example, we only -execute 5 iterations of the algorithm. Although the algorithm does not converge -completely, we can see that the first-stage variables already exhibit -relatively good agreement: -.. testcode:: +.. code-block:: bash + + mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --solver-name gurobi_persistent --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01 - variables = ph.gather_var_values_to_rank0() - for (scenario_name, variable_name) in variables: - variable_value = variables[scenario_name, variable_name] - print(scenario_name, variable_name, variable_value) - -.. testoutput:: - :hide: - - ... - average X[BEETS] - ... - -.. testoutput:: - :options: +SKIP - good X[BEETS] 280.6489711937925 - good X[CORN] 85.26131687116064 - good X[WHEAT] 134.0897119350402 - average X[BEETS] 283.2796296293019 - average X[CORN] 80.00000000014425 - average X[WHEAT] 136.72037037055298 - bad X[BEETS] 280.64897119379475 - bad X[CORN] 85.26131687116226 - bad X[WHEAT] 134.08971193504266 - -The function ``gather_var_values_to_rank0`` can be used in parallel to collect -the values of all non-anticipative variables at the root. In this (serial) -example, it simply returns the values of the first-stage variables. Solving Using Benders' Decomposition ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/src/quick_start.rst b/doc/src/quick_start.rst index 7799a23c7..07e6ef9c8 100644 --- a/doc/src/quick_start.rst +++ b/doc/src/quick_start.rst @@ -7,9 +7,9 @@ If you installed from github, install from source using pip in editable mode fro pip install -e . -This step is not needed if you installed using pip. -You can also include the extras flags ``mpi`` to install a compliant version -of mpi4py or ``docs`` to install documentation dependencies from pip. +This step is not needed if you installed using pip, but we recommend against installing using pip because +the software is under active development and the pip version is almost always out of date. +You can also include the extras flag ``docs`` to install documentation dependencies from pip. Verify installation @@ -25,13 +25,8 @@ terminal commands: cd examples cd farmer - python farmer_ef 1 3 solver_name + python ../../mpisppy/generic_cylinders.py --module-name farmer --help -but replace `solver_name` with the name of the solver you have installed, e.g., if you have installed glpk, use - -:: - - python farmer_ef 1 3 glpk If you intend to use any parallel features, you should verify that you have a *proper* installation of MPI and ``mpi4py``; see the section @@ -50,15 +45,10 @@ first thing is to code a scenario creation function. See If you create a few more helper functions (see :ref:`helper_functions`), you can make use of the ``generic_cylinders`` program (see :ref:`generic_cylinders`) to use the hub and spoke system or to solve the the EF directly. - -Alternatively, once you have the scenario creation function, -you can mimic the code in ``examples.farmer.farmer_ef`` to -solve the extensive form directly. If you want to use the hub -and spoke system to solve your problem via decomposition, you -should proceed to the section on writing :ref:`Drivers`, or to -the :ref:`Examples` section, or to the :ref:`generic_cylinders` section. - +You might want to start with the farmer example. See the `farmer` directory in the `examples` directory for the +files `farmer.py` and `farmer_generic.bash`. For more information see :ref:`Examples`. + PySP Users ---------- @@ -74,12 +64,6 @@ Here are the general steps: # Construct a ``PySPModel`` object giving its constructor information about your PySP model. -# Create an options dictionary. - -# Create a PH or EF ``mpi-sppy`` object. - -# Call its main function. - These steps alone will not result in use of the hub-spoke features of `mpi-sppy`, but they will get your PySP model running in ``mpi-sppy``. See ``examples/farmer/from_pysp`` for some @@ -95,4 +79,4 @@ The quickest thing to do is to run one of the canned examples that comes with ``mpi-sppy``. They are in subdirectories of ``examples`` and sample commands can be obtained by looking at the code in ``examples.runall.py``. There is a table in the -mpi-sppy paper that gives references for some of the examples. +mpi-sppy paper in MPC that gives references for some of the examples. diff --git a/doc/src/seqsamp.rst b/doc/src/seqsamp.rst index ac92726e7..3fdb0acd5 100644 --- a/doc/src/seqsamp.rst +++ b/doc/src/seqsamp.rst @@ -27,8 +27,8 @@ For multi-stage, use `multi_seqsampling.py`. Examples -------- -There is sample code for two-stage, sequential sampling in ``examples.farmer.farmer_seqsampling.py`` and -a bash scrip to test drive it is ``examples.farmer.farmer_sequential.bash``. +There is sample code for two-stage, sequential sampling in ``examples.farmer.CI.farmer_seqsampling.py`` and +a bash scrip to test drive it is ``examples.farmer.CI.farmer_sequential.bash``. There is sample code for multi-stage, sequential sampling in ``examples.aircond.aircond_seqsampling.py`` and a bash scrip to test drive it is ``examples.aircond.aircond_sequential.bash``. diff --git a/doc/src/spokes.rst b/doc/src/spokes.rst index bc82b66a2..6c3f47c69 100644 --- a/doc/src/spokes.rst +++ b/doc/src/spokes.rst @@ -68,9 +68,9 @@ so as to obtain better outer (lower when minimizing) bounds. It can also provide a Lagrangian bound even if the hub does not provide lagrangian multipliers. The easiest way to use ``ph_ob`` is via the vanilla ``ph_ob_spoke`` method -as illustrated in ``examples.farmer_cylinders.py``. This method takes values +as illustrated in ``examples.farmer.archive.farmer_cylinders.py``. This method takes values from the config object (assuming the config object's ``ph_ob`` method -was called as shown in the function ``examples.farmer_cylinders._parse_args``) +was called as shown in the function ``examples.farmer.archive.farmer_cylinders._parse_args``) and sets up the options for the spoke. The option ``ph-ob-initial-rho-rescale-factor`` defaults to 0.1, so if nothing diff --git a/doc/src/zhat.rst b/doc/src/zhat.rst index 72f7a30c1..cdd077f1f 100644 --- a/doc/src/zhat.rst +++ b/doc/src/zhat.rst @@ -13,14 +13,14 @@ Unless you are directly using mid-level functionality, your code may need to write the root node nonanticipative variable values (called `xhat` or `xhat_one`) to a file for later processing. This is typically done using ``sputils.ef_ROOT_nonants_npy_serializer``, which -is shown in various examples, e.g., ``examples.farmer.farmer.py`` +is shown in various examples, e.g., ``examples.farmer.archive.farmer.py`` zhat4xhat --------- The program ``zhat4xhat`` estimates approximate confidence intervals for the objective function value, zhat, given an xhat. See -``examples.farmer.farmer_zhat.bash`` for a bash script that first +``examples.farmer.CI.farmer_zhat.bash`` for a bash script that first creates the xhat file, then computes an out-of-sample confidence interval for it. Note: this program does not compute a confidence interval for zstar, which is done using software documented in diff --git a/examples/afew.py b/examples/afew.py index 04aaf341e..d5d72cecd 100644 --- a/examples/afew.py +++ b/examples/afew.py @@ -41,11 +41,14 @@ def do_one(dirname, progname, np, argstring): badguys[dirname] = [runstring] else: badguys[dirname].append(runstring) - os.chdir("..") + if '/' not in dirname: + os.chdir("..") + else: + os.chdir("../..") # hack for one level of subdirectories # for farmer, the first arg is num_scens and is required -do_one("farmer", "farmer_cylinders.py", 3, +do_one("farmer/archive", "farmer_cylinders.py", 3, "--num-scens=3 --bundles-per-rank=0 --max-iterations=50 " "--default-rho=1 --sep-rho --display-convergence-detail " "--solver-name={} --xhatshuffle --lagrangian --use-norm-rho-updater".format(solver_name)) diff --git a/examples/farmer/CI/Readme.rst b/examples/farmer/CI/Readme.rst new file mode 100644 index 000000000..315c5dba3 --- /dev/null +++ b/examples/farmer/CI/Readme.rst @@ -0,0 +1,4 @@ +CI +-- + +This directory contains code related to getting confidence intervals for the farmer example. diff --git a/examples/farmer/CI/farmer.py b/examples/farmer/CI/farmer.py new file mode 100644 index 000000000..6fe141a55 --- /dev/null +++ b/examples/farmer/CI/farmer.py @@ -0,0 +1,281 @@ +############################################################################### +# 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. +############################################################################### +# special for ph debugging DLW Dec 2018 +# unlimited crops +# ALL INDEXES ARE ZERO-BASED +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright 2018 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ +# +# special scalable farmer for stress-testing + +import pyomo.environ as pyo +import numpy as np +import mpisppy.utils.sputils as sputils + +# Use this random stream: +farmerstream = np.random.RandomState() + +def scenario_creator( + scenario_name, use_integer=False, sense=pyo.minimize, crops_multiplier=1, + num_scens=None, seedoffset=0 +): + """ Create a scenario for the (scalable) farmer example. + + Args: + scenario_name (str): + Name of the scenario to construct. + use_integer (bool, optional): + If True, restricts variables to be integer. Default is False. + sense (int, optional): + Model sense (minimization or maximization). Must be either + pyo.minimize or pyo.maximize. Default is pyo.minimize. + crops_multiplier (int, optional): + Factor to control scaling. There will be three times this many + crops. Default is 1. + num_scens (int, optional): + Number of scenarios. We use it to compute _mpisppy_probability. + Default is None. + seedoffset (int): used by confidence interval code to create replicates + """ + # scenario_name has the form e.g. scen12, foobar7 + # The digits are scraped off the right of scenario_name using regex then + # converted mod 3 into one of the below avg./avg./above avg. scenarios + scennum = sputils.extract_num(scenario_name) + basenames = ['BelowAverageScenario', 'AverageScenario', 'AboveAverageScenario'] + basenum = scennum % 3 + groupnum = scennum // 3 + scenario_name = basenames[basenum]+str(groupnum) + + # The RNG is seeded with the scenario number so that it is + # reproducible when used with multiple threads. + farmerstream.seed(scennum+seedoffset) + + # Check for minimization vs. maximization + if sense not in [pyo.minimize, pyo.maximize]: + raise ValueError("Model sense Not recognized") + + # Create the concrete model object + scengroupnum = sputils.extract_num(scenario_name) + scenario_base_name = scenario_name.rstrip("0123456789") + + model = pyo.ConcreteModel(scenario_name) + + def crops_init(m): + retval = [] + for i in range(crops_multiplier): + retval.append("WHEAT"+str(i)) + retval.append("CORN"+str(i)) + retval.append("SUGAR_BEETS"+str(i)) + return retval + + model.CROPS = pyo.Set(initialize=crops_init) + + # + # Parameters + # + + model.TOTAL_ACREAGE = 500.0 * crops_multiplier + + def _scale_up_data(indict): + outdict = {} + for i in range(crops_multiplier): + for crop in ['WHEAT', 'CORN', 'SUGAR_BEETS']: + outdict[crop+str(i)] = indict[crop] + return outdict + + model.PriceQuota = _scale_up_data( + {'WHEAT':100000.0,'CORN':100000.0,'SUGAR_BEETS':6000.0}) + + model.SubQuotaSellingPrice = _scale_up_data( + {'WHEAT':170.0,'CORN':150.0,'SUGAR_BEETS':36.0}) + + model.SuperQuotaSellingPrice = _scale_up_data( + {'WHEAT':0.0,'CORN':0.0,'SUGAR_BEETS':10.0}) + + model.CattleFeedRequirement = _scale_up_data( + {'WHEAT':200.0,'CORN':240.0,'SUGAR_BEETS':0.0}) + + model.PurchasePrice = _scale_up_data( + {'WHEAT':238.0,'CORN':210.0,'SUGAR_BEETS':100000.0}) + + model.PlantingCostPerAcre = _scale_up_data( + {'WHEAT':150.0,'CORN':230.0,'SUGAR_BEETS':260.0}) + + # + # Stochastic Data + # + Yield = {} + Yield['BelowAverageScenario'] = \ + {'WHEAT':2.0,'CORN':2.4,'SUGAR_BEETS':16.0} + Yield['AverageScenario'] = \ + {'WHEAT':2.5,'CORN':3.0,'SUGAR_BEETS':20.0} + Yield['AboveAverageScenario'] = \ + {'WHEAT':3.0,'CORN':3.6,'SUGAR_BEETS':24.0} + + def Yield_init(m, cropname): + # yield as in "crop yield" + crop_base_name = cropname.rstrip("0123456789") + if scengroupnum != 0: + return Yield[scenario_base_name][crop_base_name]+farmerstream.rand() + else: + return Yield[scenario_base_name][crop_base_name] + + model.Yield = pyo.Param(model.CROPS, + within=pyo.NonNegativeReals, + initialize=Yield_init, + mutable=True) + + # + # Variables + # + + if (use_integer): + model.DevotedAcreage = pyo.Var(model.CROPS, + within=pyo.NonNegativeIntegers, + bounds=(0.0, model.TOTAL_ACREAGE)) + else: + model.DevotedAcreage = pyo.Var(model.CROPS, + bounds=(0.0, model.TOTAL_ACREAGE)) + + model.QuantitySubQuotaSold = pyo.Var(model.CROPS, bounds=(0.0, None)) + model.QuantitySuperQuotaSold = pyo.Var(model.CROPS, bounds=(0.0, None)) + model.QuantityPurchased = pyo.Var(model.CROPS, bounds=(0.0, None)) + + # + # Constraints + # + + def ConstrainTotalAcreage_rule(model): + return pyo.sum_product(model.DevotedAcreage) <= model.TOTAL_ACREAGE + + model.ConstrainTotalAcreage = pyo.Constraint(rule=ConstrainTotalAcreage_rule) + + def EnforceCattleFeedRequirement_rule(model, i): + return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i] + + model.EnforceCattleFeedRequirement = pyo.Constraint(model.CROPS, rule=EnforceCattleFeedRequirement_rule) + + def LimitAmountSold_rule(model, i): + return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0 + + model.LimitAmountSold = pyo.Constraint(model.CROPS, rule=LimitAmountSold_rule) + + def EnforceQuotas_rule(model, i): + return (0.0, model.QuantitySubQuotaSold[i], model.PriceQuota[i]) + + model.EnforceQuotas = pyo.Constraint(model.CROPS, rule=EnforceQuotas_rule) + + # Stage-specific cost computations; + + def ComputeFirstStageCost_rule(model): + return pyo.sum_product(model.PlantingCostPerAcre, model.DevotedAcreage) + model.FirstStageCost = pyo.Expression(rule=ComputeFirstStageCost_rule) + + def ComputeSecondStageCost_rule(model): + expr = pyo.sum_product(model.PurchasePrice, model.QuantityPurchased) + expr -= pyo.sum_product(model.SubQuotaSellingPrice, model.QuantitySubQuotaSold) + expr -= pyo.sum_product(model.SuperQuotaSellingPrice, model.QuantitySuperQuotaSold) + return expr + model.SecondStageCost = pyo.Expression(rule=ComputeSecondStageCost_rule) + + def total_cost_rule(model): + if (sense == pyo.minimize): + return model.FirstStageCost + model.SecondStageCost + return -model.FirstStageCost - model.SecondStageCost + model.Total_Cost_Objective = pyo.Objective(rule=total_cost_rule, + sense=sense) + + + # Create the list of nodes associated with the scenario. + varlist = [model.DevotedAcreage] + sputils.attach_root_node(model, model.FirstStageCost, varlist) + + #Add the probability of the scenario + if num_scens is not None : + model._mpisppy_probability = 1/num_scens + else: + model._mpisppy_probability = "uniform" + return model + + +# begin helper functions +#========= +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 + if (start is None) : + start=0 + return [f"scen{i}" for i in range(start,start+num_scens)] + + +#========= +def inparser_adder(cfg): + # add options unique to farmer + cfg.num_scens_required() + cfg.add_to_config("crops_multiplier", + description="number of crops will be three times this (default 1)", + domain=int, + default=1) + + cfg.add_to_config("farmer_with_integers", + description="make the version that has integers (default False)", + domain=bool, + default=False) + + +#========= +def kw_creator(cfg): + # (for Amalgamator): linked to the scenario_creator and inparser_adder + kwargs = {"use_integer": cfg.get('farmer_with_integers', False), + "crops_multiplier": cfg.get('crops_multiplier', 1), + "num_scens" : cfg.get('num_scens', None), + } + return kwargs + +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + """ Create a scenario within a sample tree. Mainly for multi-stage and simple for two-stage. + (this function supports zhat and confidence interval code) + Args: + sname (string): scenario name to be created + stage (int >=1 ): for stages > 1, fix data based on sname in earlier stages + sample_branching_factors (list of ints): branching factors for the sample tree + seed (int): To allow random sampling (for some problems, it might be scenario offset) + given_scenario (Pyomo concrete model): if not None, use this to get data for ealier stages + scenario_creator_kwargs (dict): keyword args for the standard scenario creator funcion + Returns: + scenario (Pyomo concrete model): A scenario for sname with data in stages < stage determined + by the arguments + """ + # Since this is a two-stage problem, we don't have to do much. + sca = scenario_creator_kwargs.copy() + sca["seedoffset"] = seed + sca["num_scens"] = sample_branching_factors[0] # two-stage problem + return scenario_creator(sname, **sca) + + +# end functions not needed by farmer_cylinders + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + sname = scenario_name + s = scenario + if sname == 'scen0': + print("Arbitrary sanity checks:") + print ("SUGAR_BEETS0 for scenario",sname,"is", + pyo.value(s.DevotedAcreage["SUGAR_BEETS0"])) + print ("FirstStageCost for scenario",sname,"is", pyo.value(s.FirstStageCost)) diff --git a/examples/farmer/farmer_ef.py b/examples/farmer/CI/farmer_ef.py similarity index 93% rename from examples/farmer/farmer_ef.py rename to examples/farmer/CI/farmer_ef.py index df4edb58b..9ba6974df 100644 --- a/examples/farmer/farmer_ef.py +++ b/examples/farmer/CI/farmer_ef.py @@ -6,6 +6,9 @@ # All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for # full copyright and license information. ############################################################################### +# NOTE: You probably don't want to look at this file. You probably should be +# using generic_cylinders.py. See the farmer_generic.bash for examples. +# (This file is mainly of interest for confidence intervals) # This shows two ways to get parameters to the EF for solution; both are fairly short. import sys diff --git a/examples/farmer/farmer_mmw.bash b/examples/farmer/CI/farmer_mmw.bash similarity index 100% rename from examples/farmer/farmer_mmw.bash rename to examples/farmer/CI/farmer_mmw.bash diff --git a/examples/farmer/farmer_rho_demo.bash b/examples/farmer/CI/farmer_rho_demo.bash similarity index 100% rename from examples/farmer/farmer_rho_demo.bash rename to examples/farmer/CI/farmer_rho_demo.bash diff --git a/examples/farmer/farmer_rho_demo.py b/examples/farmer/CI/farmer_rho_demo.py similarity index 100% rename from examples/farmer/farmer_rho_demo.py rename to examples/farmer/CI/farmer_rho_demo.py diff --git a/examples/farmer/farmer_seqsampling.py b/examples/farmer/CI/farmer_seqsampling.py similarity index 100% rename from examples/farmer/farmer_seqsampling.py rename to examples/farmer/CI/farmer_seqsampling.py diff --git a/examples/farmer/farmer_sequential.bash b/examples/farmer/CI/farmer_sequential.bash similarity index 100% rename from examples/farmer/farmer_sequential.bash rename to examples/farmer/CI/farmer_sequential.bash diff --git a/examples/farmer/farmer_zhat.bash b/examples/farmer/CI/farmer_zhat.bash similarity index 100% rename from examples/farmer/farmer_zhat.bash rename to examples/farmer/CI/farmer_zhat.bash diff --git a/examples/farmer/mmw_slurm.bash b/examples/farmer/CI/mmw_slurm.bash similarity index 100% rename from examples/farmer/mmw_slurm.bash rename to examples/farmer/CI/mmw_slurm.bash diff --git a/examples/farmer/cs_farmer.py b/examples/farmer/archive/cs_farmer.py similarity index 100% rename from examples/farmer/cs_farmer.py rename to examples/farmer/archive/cs_farmer.py diff --git a/examples/farmer/archive/farmer.py b/examples/farmer/archive/farmer.py new file mode 100644 index 000000000..8f40726d8 --- /dev/null +++ b/examples/farmer/archive/farmer.py @@ -0,0 +1,303 @@ +############################################################################### +# 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. +############################################################################### +# special for ph debugging DLW Dec 2018 +# unlimited crops +# ALL INDEXES ARE ZERO-BASED +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright 2018 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ +# +# special scalable farmer for stress-testing + +import pyomo.environ as pyo +import numpy as np +import mpisppy.utils.sputils as sputils + +# Use this random stream: +farmerstream = np.random.RandomState() + +def scenario_creator( + scenario_name, use_integer=False, sense=pyo.minimize, crops_multiplier=1, + num_scens=None, seedoffset=0 +): + """ Create a scenario for the (scalable) farmer example. + + Args: + scenario_name (str): + Name of the scenario to construct. + use_integer (bool, optional): + If True, restricts variables to be integer. Default is False. + sense (int, optional): + Model sense (minimization or maximization). Must be either + pyo.minimize or pyo.maximize. Default is pyo.minimize. + crops_multiplier (int, optional): + Factor to control scaling. There will be three times this many + crops. Default is 1. + num_scens (int, optional): + Number of scenarios. We use it to compute _mpisppy_probability. + Default is None. + seedoffset (int): used by confidence interval code + """ + # scenario_name has the form e.g. scen12, foobar7 + # The digits are scraped off the right of scenario_name using regex then + # converted mod 3 into one of the below avg./avg./above avg. scenarios + scennum = sputils.extract_num(scenario_name) + basenames = ['BelowAverageScenario', 'AverageScenario', 'AboveAverageScenario'] + basenum = scennum % 3 + groupnum = scennum // 3 + scenname = basenames[basenum]+str(groupnum) + + # The RNG is seeded with the scenario number so that it is + # reproducible when used with multiple threads. + # NOTE: if you want to do replicates, you will need to pass a seed + # as a kwarg to scenario_creator then use seed+scennum as the seed argument. + farmerstream.seed(scennum+seedoffset) + + # Check for minimization vs. maximization + if sense not in [pyo.minimize, pyo.maximize]: + raise ValueError("Model sense Not recognized") + + # Create the concrete model object + model = pysp_instance_creation_callback( + scenname, + use_integer=use_integer, + sense=sense, + crops_multiplier=crops_multiplier, + ) + + # Create the list of nodes associated with the scenario (for two stage, + # there is only one node associated with the scenario--leaf nodes are + # ignored). + varlist = [model.DevotedAcreage] + sputils.attach_root_node(model, model.FirstStageCost, varlist) + + #Add the probability of the scenario + if num_scens is not None : + model._mpisppy_probability = 1/num_scens + else: + model._mpisppy_probability = "uniform" + return model + +def pysp_instance_creation_callback( + scenario_name, use_integer=False, sense=pyo.minimize, crops_multiplier=1 +): + # long function to create the entire model + # scenario_name is a string (e.g. AboveAverageScenario0) + # + # Returns a concrete model for the specified scenario + + # scenarios come in groups of three + scengroupnum = sputils.extract_num(scenario_name) + scenario_base_name = scenario_name.rstrip("0123456789") + + model = pyo.ConcreteModel(scenario_name) + + def crops_init(m): + retval = [] + for i in range(crops_multiplier): + retval.append("WHEAT"+str(i)) + retval.append("CORN"+str(i)) + retval.append("SUGAR_BEETS"+str(i)) + return retval + + model.CROPS = pyo.Set(initialize=crops_init) + + # + # Parameters + # + + model.TOTAL_ACREAGE = 500.0 * crops_multiplier + + def _scale_up_data(indict): + outdict = {} + for i in range(crops_multiplier): + for crop in ['WHEAT', 'CORN', 'SUGAR_BEETS']: + outdict[crop+str(i)] = indict[crop] + return outdict + + model.PriceQuota = _scale_up_data( + {'WHEAT':100000.0,'CORN':100000.0,'SUGAR_BEETS':6000.0}) + + model.SubQuotaSellingPrice = _scale_up_data( + {'WHEAT':170.0,'CORN':150.0,'SUGAR_BEETS':36.0}) + + model.SuperQuotaSellingPrice = _scale_up_data( + {'WHEAT':0.0,'CORN':0.0,'SUGAR_BEETS':10.0}) + + model.CattleFeedRequirement = _scale_up_data( + {'WHEAT':200.0,'CORN':240.0,'SUGAR_BEETS':0.0}) + + model.PurchasePrice = _scale_up_data( + {'WHEAT':238.0,'CORN':210.0,'SUGAR_BEETS':100000.0}) + + model.PlantingCostPerAcre = _scale_up_data( + {'WHEAT':150.0,'CORN':230.0,'SUGAR_BEETS':260.0}) + + # + # Stochastic Data + # + Yield = {} + Yield['BelowAverageScenario'] = \ + {'WHEAT':2.0,'CORN':2.4,'SUGAR_BEETS':16.0} + Yield['AverageScenario'] = \ + {'WHEAT':2.5,'CORN':3.0,'SUGAR_BEETS':20.0} + Yield['AboveAverageScenario'] = \ + {'WHEAT':3.0,'CORN':3.6,'SUGAR_BEETS':24.0} + + def Yield_init(m, cropname): + # yield as in "crop yield" + crop_base_name = cropname.rstrip("0123456789") + if scengroupnum != 0: + return Yield[scenario_base_name][crop_base_name]+farmerstream.rand() + else: + return Yield[scenario_base_name][crop_base_name] + + model.Yield = pyo.Param(model.CROPS, + within=pyo.NonNegativeReals, + initialize=Yield_init, + mutable=True) + + # + # Variables + # + + if (use_integer): + model.DevotedAcreage = pyo.Var(model.CROPS, + within=pyo.NonNegativeIntegers, + bounds=(0.0, model.TOTAL_ACREAGE)) + else: + model.DevotedAcreage = pyo.Var(model.CROPS, + bounds=(0.0, model.TOTAL_ACREAGE)) + + model.QuantitySubQuotaSold = pyo.Var(model.CROPS, bounds=(0.0, None)) + model.QuantitySuperQuotaSold = pyo.Var(model.CROPS, bounds=(0.0, None)) + model.QuantityPurchased = pyo.Var(model.CROPS, bounds=(0.0, None)) + + # + # Constraints + # + + def ConstrainTotalAcreage_rule(model): + return pyo.sum_product(model.DevotedAcreage) <= model.TOTAL_ACREAGE + + model.ConstrainTotalAcreage = pyo.Constraint(rule=ConstrainTotalAcreage_rule) + + def EnforceCattleFeedRequirement_rule(model, i): + return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i] + + model.EnforceCattleFeedRequirement = pyo.Constraint(model.CROPS, rule=EnforceCattleFeedRequirement_rule) + + def LimitAmountSold_rule(model, i): + return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0 + + model.LimitAmountSold = pyo.Constraint(model.CROPS, rule=LimitAmountSold_rule) + + def EnforceQuotas_rule(model, i): + return (0.0, model.QuantitySubQuotaSold[i], model.PriceQuota[i]) + + model.EnforceQuotas = pyo.Constraint(model.CROPS, rule=EnforceQuotas_rule) + + # Stage-specific cost computations; + + def ComputeFirstStageCost_rule(model): + return pyo.sum_product(model.PlantingCostPerAcre, model.DevotedAcreage) + model.FirstStageCost = pyo.Expression(rule=ComputeFirstStageCost_rule) + + def ComputeSecondStageCost_rule(model): + expr = pyo.sum_product(model.PurchasePrice, model.QuantityPurchased) + expr -= pyo.sum_product(model.SubQuotaSellingPrice, model.QuantitySubQuotaSold) + expr -= pyo.sum_product(model.SuperQuotaSellingPrice, model.QuantitySuperQuotaSold) + return expr + model.SecondStageCost = pyo.Expression(rule=ComputeSecondStageCost_rule) + + def total_cost_rule(model): + if (sense == pyo.minimize): + return model.FirstStageCost + model.SecondStageCost + return -model.FirstStageCost - model.SecondStageCost + model.Total_Cost_Objective = pyo.Objective(rule=total_cost_rule, + sense=sense) + + return model + +# begin functions not needed by farmer_cylinders +# (but needed by special codes such as confidence intervals) +#========= +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 + if (start is None) : + start=0 + return [f"scen{i}" for i in range(start,start+num_scens)] + + + +#========= +def inparser_adder(cfg): + # add options unique to farmer + cfg.num_scens_required() + cfg.add_to_config("crops_multiplier", + description="number of crops will be three times this (default 1)", + domain=int, + default=1) + + cfg.add_to_config("farmer_with_integers", + description="make the version that has integers (default False)", + domain=bool, + default=False) + + +#========= +def kw_creator(cfg): + # (for Amalgamator): linked to the scenario_creator and inparser_adder + kwargs = {"use_integer": cfg.get('farmer_with_integers', False), + "crops_multiplier": cfg.get('crops_multiplier', 1), + "num_scens" : cfg.get('num_scens', None), + } + return kwargs + +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + """ Create a scenario within a sample tree. Mainly for multi-stage and simple for two-stage. + (this function supports zhat and confidence interval code) + Args: + sname (string): scenario name to be created + stage (int >=1 ): for stages > 1, fix data based on sname in earlier stages + sample_branching_factors (list of ints): branching factors for the sample tree + seed (int): To allow random sampling (for some problems, it might be scenario offset) + given_scenario (Pyomo concrete model): if not None, use this to get data for ealier stages + scenario_creator_kwargs (dict): keyword args for the standard scenario creator funcion + Returns: + scenario (Pyomo concrete model): A scenario for sname with data in stages < stage determined + by the arguments + """ + # Since this is a two-stage problem, we don't have to do much. + sca = scenario_creator_kwargs.copy() + sca["seedoffset"] = seed + sca["num_scens"] = sample_branching_factors[0] # two-stage problem + return scenario_creator(sname, **sca) + + +# end functions not needed by farmer_cylinders + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + sname = scenario_name + s = scenario + if sname == 'scen0': + print("Arbitrary sanity checks:") + print ("SUGAR_BEETS0 for scenario",sname,"is", + pyo.value(s.DevotedAcreage["SUGAR_BEETS0"])) + print ("FirstStageCost for scenario",sname,"is", pyo.value(s.FirstStageCost)) diff --git a/examples/farmer/farmer_ama.py b/examples/farmer/archive/farmer_ama.py similarity index 100% rename from examples/farmer/farmer_ama.py rename to examples/farmer/archive/farmer_ama.py diff --git a/examples/farmer/farmer_cylinders.py b/examples/farmer/archive/farmer_cylinders.py similarity index 100% rename from examples/farmer/farmer_cylinders.py rename to examples/farmer/archive/farmer_cylinders.py diff --git a/examples/farmer/farmer_ph.py b/examples/farmer/archive/farmer_ph.py similarity index 100% rename from examples/farmer/farmer_ph.py rename to examples/farmer/archive/farmer_ph.py diff --git a/examples/farmer/farmer.py b/examples/farmer/farmer.py index 8f40726d8..6fe141a55 100644 --- a/examples/farmer/farmer.py +++ b/examples/farmer/farmer.py @@ -48,7 +48,7 @@ def scenario_creator( num_scens (int, optional): Number of scenarios. We use it to compute _mpisppy_probability. Default is None. - seedoffset (int): used by confidence interval code + seedoffset (int): used by confidence interval code to create replicates """ # scenario_name has the form e.g. scen12, foobar7 # The digits are scraped off the right of scenario_name using regex then @@ -57,12 +57,10 @@ def scenario_creator( basenames = ['BelowAverageScenario', 'AverageScenario', 'AboveAverageScenario'] basenum = scennum % 3 groupnum = scennum // 3 - scenname = basenames[basenum]+str(groupnum) + scenario_name = basenames[basenum]+str(groupnum) # The RNG is seeded with the scenario number so that it is # reproducible when used with multiple threads. - # NOTE: if you want to do replicates, you will need to pass a seed - # as a kwarg to scenario_creator then use seed+scennum as the seed argument. farmerstream.seed(scennum+seedoffset) # Check for minimization vs. maximization @@ -70,35 +68,6 @@ def scenario_creator( raise ValueError("Model sense Not recognized") # Create the concrete model object - model = pysp_instance_creation_callback( - scenname, - use_integer=use_integer, - sense=sense, - crops_multiplier=crops_multiplier, - ) - - # Create the list of nodes associated with the scenario (for two stage, - # there is only one node associated with the scenario--leaf nodes are - # ignored). - varlist = [model.DevotedAcreage] - sputils.attach_root_node(model, model.FirstStageCost, varlist) - - #Add the probability of the scenario - if num_scens is not None : - model._mpisppy_probability = 1/num_scens - else: - model._mpisppy_probability = "uniform" - return model - -def pysp_instance_creation_callback( - scenario_name, use_integer=False, sense=pyo.minimize, crops_multiplier=1 -): - # long function to create the entire model - # scenario_name is a string (e.g. AboveAverageScenario0) - # - # Returns a concrete model for the specified scenario - - # scenarios come in groups of three scengroupnum = sputils.extract_num(scenario_name) scenario_base_name = scenario_name.rstrip("0123456789") @@ -229,10 +198,20 @@ def total_cost_rule(model): model.Total_Cost_Objective = pyo.Objective(rule=total_cost_rule, sense=sense) + + # Create the list of nodes associated with the scenario. + varlist = [model.DevotedAcreage] + sputils.attach_root_node(model, model.FirstStageCost, varlist) + + #Add the probability of the scenario + if num_scens is not None : + model._mpisppy_probability = 1/num_scens + else: + model._mpisppy_probability = "uniform" return model -# begin functions not needed by farmer_cylinders -# (but needed by special codes such as confidence intervals) + +# begin helper functions #========= def scenario_names_creator(num_scens,start=None): # (only for Amalgamator): return the full list of num_scens scenario names @@ -242,7 +221,6 @@ def scenario_names_creator(num_scens,start=None): return [f"scen{i}" for i in range(start,start+num_scens)] - #========= def inparser_adder(cfg): # add options unique to farmer diff --git a/examples/farmer/farmer_generic.bash b/examples/farmer/farmer_generic.bash new file mode 100644 index 000000000..9decf4894 --- /dev/null +++ b/examples/farmer/farmer_generic.bash @@ -0,0 +1,18 @@ +#!/bin/bash +# Illustrate the use of generic cylinders for the farmer example. +# We assumng the current directory is examples.farmer + +SOLVER=gurobi + +# Note that to get help, you still need to specify a module-name +python ../../mpisppy/generic_cylinders.py --module-name farmer --help + +# Here is a simple ef command, that writes solution data to two files with the base name farmersol +# and to a directory named farmersol_soldir +echo "Starting EF" +python ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --EF --EF-solver-name gurobi --solution-base-name farmersol + +# Here is a very simple PH command that also computes bounds +echo "Starting PH" +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 + diff --git a/examples/farmer/schur_complement.py b/examples/farmer/schur/schur_complement.py similarity index 100% rename from examples/farmer/schur_complement.py rename to examples/farmer/schur/schur_complement.py diff --git a/examples/generic_cylinders.bash b/examples/generic_cylinders.bash index 92bcf1b4a..6c0262405 100644 --- a/examples/generic_cylinders.bash +++ b/examples/generic_cylinders.bash @@ -1,14 +1,15 @@ #!/bin/bash # This runs a few command lines to illustrate the use of generic_cylinders.py +##set -e -SOLVER="cplex" -SPB=1 +SOLVER="gurobi" +SPB=2 echo "^^^ hub only with w-writer (smoke) ^^^" -python -m mpi4py ../mpisppy/generic_cylinders.py --module-name farmer/farmer --num-scens 3 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --W-writer --W-fname w_values.csv +python -m mpi4py ../mpisppy/generic_cylinders.py --module-name farmer/farmer --num-scens 3 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --W-and-xbar-writer --W-fname w_values.csv echo "^^^ hub only with w-reader (smoke) ^^^" -python -m mpi4py ../mpisppy/generic_cylinders.py --module-name farmer/farmer --num-scens 3 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --W-reader --init-W-fname w_values.csv +python -m mpi4py ../mpisppy/generic_cylinders.py --module-name farmer/farmer --num-scens 3 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --W-and-xbar-reader --init-W-fname w_values.csv echo "^^^ Multi-stage AirCond ^^^" mpiexec -np 3 python -m mpi4py ../mpisppy/generic_cylinders.py --module-name mpisppy.tests.examples.aircond --branching-factors "3 3 3" --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatxbar --rel-gap 0.01 --solution-base-name aircond_nonants @@ -27,7 +28,7 @@ mpiexec -np 3 python -m mpi4py ../mpisppy/generic_cylinders.py --module-name mpi 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 +python ../../mpisppy/generic_cylinders.py --module-name sizes --num-scens 3 --default-rho 1 --solver-name ${SOLVER} --max-iterations 0 --write-scenario-lp-mps-files cd .. echo "^^^ pickle sizes bundles ^^^" @@ -39,7 +40,7 @@ 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 +python ../../mpisppy/generic_cylinders.py --module-name sizes --num-scens 10 --default-rho 1 --solver-name ${SOLVER} --max-iterations 0 --write-scenario-lp-mps-files --unpickle-bundles-dir sizes_pickles --scenarios-per-bundle 5 cd .. echo "^^^ pickle the scenarios ^^^" @@ -59,11 +60,6 @@ cd farmer mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 10 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01 --unpickle-scenarios-dir farmer_pickles cd .. -echo "^^^ use proper bundles without writing ^^^" -cd sslp -mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name sslp --sslp-data-path ./data --instance-name sslp_15_45_10 --scenarios-per-bundle $SPB --default-rho 1 --solver-name ${SOLVER} --max-iterations 5 --lagrangian --xhatshuffle --rel-gap 0.001 -cd .. - echo "^^^ write pickle bundles ^^^" cd sslp python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name sslp --sslp-data-path ./data --instance-name sslp_15_45_10 --pickle-bundles-dir sslp_pickles --scenarios-per-bundle $SPB --default-rho 1 @@ -91,11 +87,11 @@ 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 -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 +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 --rel-gap 0.001 # now do it again, but this time using dynamic rho echo "^^^ grad 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 --grad-rho-setter --rel-gap 0.001 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1 +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 --rel-gap 0.001 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1 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 @@ -173,6 +169,10 @@ mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name cd .. +echo "^^^ use proper bundles without writing (this sslp run will take some time)^^^" +cd sslp +mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name sslp --sslp-data-path ./data --instance-name sslp_15_45_10 --scenarios-per-bundle $SPB --default-rho 1 --solver-name ${SOLVER} --max-iterations 5 --lagrangian --xhatshuffle --rel-gap 0.001 +cd .. exit @@ -188,11 +188,11 @@ exit 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 -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 +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 --rel-gap 0.001 # now do it again, but this time using dynamic rho -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 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1 +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 --rel-gap 0.001 --dynamic-rho-dual-crit --dynamic-rho-dual-thresh 0.1 cd .. # end gradient based rho demo diff --git a/examples/run_all.py b/examples/run_all.py index 8e8d1544c..88cf8da83 100644 --- a/examples/run_all.py +++ b/examples/run_all.py @@ -14,7 +14,7 @@ # python run_all.py cplex # python run_all.py gurobi_persistent --oversubscribe # python run_all.py gurobi_persistent -envall nouc -# (envall does nothing; it is just a place-holder) +# (envall does nothing; it is just a place-holder; might not work with your mpiexec) import os import sys @@ -133,7 +133,7 @@ def time_one(ID, dirname, progname, np, argstring): timelistdf.to_csv(listfname, index=False) time_one.ID_check = list() -def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): +def do_one_mmw(dirname, modname, runefstring, npyfile, mmwargstring): # assume that the dirname matches the module name os.chdir(dirname) @@ -148,7 +148,7 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): # run mmw, remove .npy file else: runstring = "python -m mpisppy.confidence_intervals.mmw_conf {} --xhatpath {} {}".\ - format(dirname, npyfile, mmwargstring) + format(modname, npyfile, mmwargstring) code = os.system("echo {} && {}".format(runstring, runstring)) if code != 0: if dirname not in badguys: @@ -159,20 +159,20 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): os.remove(npyfile) os.chdir("..") -do_one("farmer", "farmer_ef.py", 1, +do_one("farmer/CI", "farmer_ef.py", 1, "1 3 {}".format(solver_name)) # for farmer_cylinders, the first arg is num_scens and is required -do_one("farmer", "farmer_cylinders.py", 3, +do_one("farmer/archive", "farmer_cylinders.py", 3, "--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", 5, +do_one("farmer/archive", "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 --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, +do_one("farmer/archive", "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)) @@ -180,15 +180,15 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): "--num-scens 3 --bundles-per-rank=0 --max-iterations=50 " "--solver-name={} --rel-gap=0.0 " "--xhatlshaped --max-solver-threads=1".format(solver_name)) -do_one("farmer", "farmer_cylinders.py", 3, +do_one("farmer/archive", "farmer_cylinders.py", 3, "--num-scens 3 --bundles-per-rank=0 --max-iterations=50 " "--default-rho=1 " "--solver-name={} --lagranger --xhatlooper".format(solver_name)) -do_one("farmer", "farmer_cylinders.py", 3, +do_one("farmer/archive", "farmer_cylinders.py", 3, "--num-scens 6 --bundles-per-rank=2 --max-iterations=50 " "--default-rho=1 --lagrangian --xhatshuffle " "--solver-name={}".format(solver_name)) -do_one("farmer", "farmer_cylinders.py", 4, +do_one("farmer/archive", "farmer_cylinders.py", 4, "--num-scens 6 --bundles-per-rank=2 --max-iterations=50 " "--fwph-stop-check-tol 0.1 " "--default-rho=1 --solver-name={} --lagrangian --xhatshuffle --fwph".format(solver_name)) @@ -197,15 +197,15 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): "--num-scens 6 --bundles-per-rank=2 --max-iterations=50 " "--ph-primal-hub --ph-dual --ph-dual-rescale-rho-factor=0.1 " "--default-rho=1 --solver-name={} --lagrangian --xhatshuffle".format(solver_name)) -do_one("farmer", "farmer_cylinders.py", 2, +do_one("farmer/archive", "farmer_cylinders.py", 2, "--num-scens 6 --bundles-per-rank=2 --max-iterations=50 " "--default-rho=1 " "--solver-name={} --xhatshuffle".format(solver_name)) -do_one("farmer", "farmer_cylinders.py", 3, +do_one("farmer/archive", "farmer_cylinders.py", 3, "--num-scens 3 --bundles-per-rank=0 --max-iterations=1 " "--default-rho=1 --tee-rank0-solves " "--solver-name={} --lagrangian --xhatshuffle".format(solver_name)) -time_one("FarmerLinProx", "farmer", "farmer_cylinders.py", 3, +time_one("FarmerLinProx", "farmer/archive", "farmer_cylinders.py", 3, "--num-scens 3 --default-rho=1.0 --max-iterations=50 " "--display-progress --rel-gap=0.0 --abs-gap=0.0 " "--linearize-proximal-terms --proximal-linearization-tolerance=1.e-6 " @@ -214,29 +214,29 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): do_one("farmer/from_pysp", "concrete_ampl.py", 1, solver_name) do_one("farmer/from_pysp", "abstract.py", 1, solver_name) -do_one("farmer", +do_one("farmer/archive", "farmer_cylinders.py", 2, f"--num-scens 3 --max-iterations=10 --default-rho=1.0 --display-progress --bundles-per-rank=0 --xhatshuffle --aph-gamma=1.0 --aph-nu=1.0 --aph-frac-needed=1.0 --aph-dispatch-frac=1.0 --abs-gap=1 --aph-sleep-seconds=0.01 --run-async --solver-name={solver_name}") -do_one("farmer", +do_one("farmer/archive", "farmer_cylinders.py", 2, f"--num-scens 3 --max-iterations=10 --default-rho=1.0 --display-progress --bundles-per-rank=0 --xhatlooper --aph-gamma=1.0 --aph-nu=1.0 --aph-frac-needed=1.0 --aph-dispatch-frac=0.25 --abs-gap=1 --display-convergence-detail --aph-sleep-seconds=0.01 --run-async --solver-name={solver_name}") -do_one("farmer", +do_one("farmer/archive", "farmer_cylinders.py", 2, f"--num-scens 30 --max-iterations=10 --default-rho=1.0 --display-progress --bundles-per-rank=0 --xhatlooper --aph-gamma=1.0 --aph-nu=1.0 --aph-frac-needed=1.0 --aph-dispatch-frac=1 --abs-gap=1 --aph-sleep-seconds=0.01 --run-async --bundles-per-rank=5 --solver-name={solver_name}") -do_one("farmer", +do_one("farmer/archive", "farmer_cylinders.py", 4, f"--num-scens 3 --bundles-per-rank=0 --max-iterations=50 --default-rho=1 --solver-name={solver_name} --lagrangian --xhatshuffle --fwph --max-stalled-iters 1") -do_one("farmer", +do_one("farmer/archive", "farmer_cylinders.py", 2, f"--num-scens 30 --max-iterations=10 --default-rho=1.0 --display-progress --bundles-per-rank=0 --xhatshuffle --aph-gamma=1.0 --aph-nu=1.0 --aph-frac-needed=1.0 --aph-dispatch-frac=0.5 --abs-gap=1 --aph-sleep-seconds=0.01 --run-async --bundles-per-rank=5 --solver-name={solver_name}") -do_one("farmer", - "../../mpisppy/generic_cylinders.py", +do_one("farmer/archive", + "../../../mpisppy/generic_cylinders.py", 4, "--module-name farmer --farmer-with-integer " "--num-scens=3 " @@ -248,18 +248,13 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): "--rel-gap=0.0 " "--solver-name={}".format(solver_name)) -do_one("farmer", - "farmer_ama.py", - 3, - f"--num-scens=10 --crops-multiplier=3 --farmer-with-integer --EF-solver-name={solver_name}") - -do_one("farmer", +do_one("farmer/CI", "farmer_seqsampling.py", 1, f"--num-scens 3 --crops-multiplier=1 --EF-solver-name={solver_name} " "--BM-h 2 --BM-q 1.3 --confidence-level 0.95 --BM-vs-BPL BM") -do_one("farmer", +do_one("farmer/CI", "farmer_seqsampling.py", 1, f"--num-scens 3 --crops-multiplier=1 --EF-solver-name={solver_name} " @@ -345,7 +340,7 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): #=========MMW TESTS========== # do_one_mmw is special -do_one_mmw("farmer", f"python farmer_ef.py 3 3 {solver_name}", "farmer_cyl_nonants.npy", f"--MMW-num-batches=5 --confidence-level 0.95 --MMW-batch-size=10 --start-scen 4 --EF-solver-name={solver_name}") +do_one_mmw("farmer/CI", "farmer", f"python farmer_ef.py 1 3 0 {solver_name}", "farmer_root_nonants.npy", f"--MMW-num-batches=5 --confidence-level 0.95 --MMW-batch-size=10 --start-scen 4 --EF-solver-name={solver_name}") #============================ diff --git a/mpisppy/generic_cylinders.py b/mpisppy/generic_cylinders.py index c48b8572e..716059596 100644 --- a/mpisppy/generic_cylinders.py +++ b/mpisppy/generic_cylinders.py @@ -143,6 +143,8 @@ def _name_lists(module, cfg, bundle_wrapper=None): #========== def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement, bundle_wrapper=None): + if cfg.get("scenarios_per_bundle") is not None and cfg.scenarios_per_bundle == 1: + raise RuntimeError("To get one scenarios-per-bundle=1, you need to write then read the 'bundles'") 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: @@ -529,7 +531,7 @@ def _write_bundles(module, inum = sputils.extract_num(module.scenario_names_creator(1)[0]) local_bundle_names = [f"Bundle_{bn*bsize+inum}_{(bn+1)*bsize-1+inum}" for bn in local_slice] - + if my_rank == 0: if os.path.exists(cfg.pickle_bundles_dir): shutil.rmtree(cfg.pickle_bundles_dir) @@ -650,6 +652,7 @@ def _proper_bundles(cfg): bundle_wrapper = None # the default if _proper_bundles(cfg): bundle_wrapper = proper_bundler.ProperBundler(module) + bundle_wrapper.set_bunBFs(cfg) 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) diff --git a/mpisppy/tests/straight_tests.py b/mpisppy/tests/straight_tests.py index 5f525b6ef..72cca9dbe 100644 --- a/mpisppy/tests/straight_tests.py +++ b/mpisppy/tests/straight_tests.py @@ -23,7 +23,7 @@ def _doone(cmdstr): ##################################################### # farmer -example_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', '..', 'examples', 'farmer') +example_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', '..', 'examples', 'farmer', 'archive') fpath = os.path.join(example_dir, 'farmer_cylinders.py') cmdstr = f"python {fpath} --help" diff --git a/mpisppy/tests/test_sc.py b/mpisppy/tests/test_sc.py index a1dfea073..1b87eb68b 100644 --- a/mpisppy/tests/test_sc.py +++ b/mpisppy/tests/test_sc.py @@ -20,7 +20,7 @@ class TestSC(unittest.TestCase): def setUp(self): self.original_path = sys.path - example_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', '..', 'examples', 'farmer') + example_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', '..', 'examples', 'farmer', 'schur') sys.path.append(example_dir) def tearDown(self): diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index a2ae9df0c..02d7a1b2d 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -149,7 +149,10 @@ def _bad_rho_setters(msg): if self.get("grad_rho") and self.get("sensi_rho"): _bad_rho_setters("Only one rho setter can be active.") - if not self.get("grad_rho") or self.get("sensi_rho") or self.get("sep_rho") or self.get("reduced_costs_rho"): + if not (self.get("grad_rho") + or self.get("sensi_rho") + or self.get("sep_rho") + or self.get("reduced_costs_rho")): if self.get("dynamic_rho_primal_crit") or self.get("dynamic_rho_dual_crit"): _bad_rho_setters("dynamic rho only works with grad-, sensi-, and sep-rho") if self.get("rc_fixer") and not self.get("reduced_costs"): diff --git a/mpisppy/utils/proper_bundler.py b/mpisppy/utils/proper_bundler.py index 7ac02985f..48ee3578e 100644 --- a/mpisppy/utils/proper_bundler.py +++ b/mpisppy/utils/proper_bundler.py @@ -30,7 +30,6 @@ # the caller needs to worry about what is in what rank # (local_scenarios might have bundle names, e.g.) - class ProperBundler(): """ Wrap model file functions so as to create proper bundles it might pickle them or read them from a pickle file @@ -54,29 +53,32 @@ def scenario_names_creator(self, num_scens, start=None, cfg=None): assert cfg is not None, "ProperBundler needs cfg for scenario names" return cfg.model.scenario_names_creator(num_scens, start=start) - def bundle_names_creator(self, num_buns, start=None, cfg=None): - - def _multistage_check(bunsize): - # returns bunBFs as a side-effect + def set_bunBFs(self, cfg): + # utility for bundle objects. Might throw if it doesn't like the branching factors. + if cfg.get("branching_factors") is None: + self.bunBFs = None + else: BFs = cfg.branching_factors beyond2size = np.prod(BFs[1:]) + bunsize = cfg.scenarios_per_bundle if bunsize % beyond2size!= 0: raise RuntimeError(f"Bundles must consume the same number of entire second stage nodes: {beyond2size=} {bunsize=}") # we need bunBFs for EF formulation self.bunBFs = [bunsize // beyond2size] + BFs[1:] - + + + def bundle_names_creator(self, num_buns, start=None, cfg=None): + # Sets self.bunBFs, which is needed by scenario_creator, as a side effect. + # start refers to the bundle number; bundles are always zero-based if start is None: start = 0 assert cfg is not None, "ProperBundler needs cfg for bundle names" assert cfg.get("num_scens") is not None assert cfg.get("scenarios_per_bundle") is not None - assert cfg.num_scens % cfg.scenarios_per_bundle == 0, "Bundles must consume the same number of entire second stage nodes: {cfg.num_scens=} {bunsize=}" + assert cfg.num_scens % cfg.scenarios_per_bundle == 0, "Bundles must consume the same number of entire second stage nodes: {cfg.num_scens=} {cfg.scenarios_per_bundle=}" bsize = cfg.scenarios_per_bundle # typing aid - if cfg.get("branching_factors") is not None: - _multistage_check(bsize) - else: - self.bunBFs = None + self.set_bunBFs(cfg) # We need to know if scenarios (not bundles) are one-based. inum = sputils.extract_num(self.module.scenario_names_creator(1)[0]) names = [f"Bundle_{bn*bsize+inum}_{(bn+1)*bsize-1+inum}" for bn in range(start+num_buns)]