Dlwdoc

.github/workflows/testconfint.yml

@@ -33,8 +33,14 @@ jobs:
         run: |
           python setup.py develop
 
-      - name: run tests
+      - name: run farmer tests
         timeout-minutes: 10
         run: |
           cd mpisppy/tests
-          python test_conf_int.py 
+          python test_conf_int_farmer.py 
+
+      - name: run aircond tests
+        timeout-minutes: 10
+        run: |
+          cd mpisppy/tests
+          python test_conf_int_aircond.py	  

doc/src/confidence_intervals.rst

@@ -4,9 +4,9 @@ MMW confidence interval
 =======================
 
 If we want to assess the quality of a given candidate solution ``xhat_one`` 
-(a first stage policy), we could try and evaluate the optimality gap, i.e. 
+(a first stage solution), we could try and evaluate the optimality gap, i.e. 
 the gap between the value of the objective function
-at ``xhat_one`` and the value of the solution to our problem.
+at ``xhat_one`` and the value of the solution to the problem.
 The class ``MMWConfidenceIntervals`` compute an estimator of the optimality gap
 as described in [mmw1999]_ (Section 3.2) and an asymptotic confidence interval for
 this gap. 
@@ -16,9 +16,9 @@ We will document two steps in the process : finding a candidate solution
 
 Sequential sampling is also supported; see `Sequential Sampling Confidence Intervals`_
 
-..note ::
-At the time of this writing, the confidence interval software assumes that the scenarios
-are presented in random order (i.e., they are not ordered by scenario number).
+.. note :: At the time of this writing, the confidence interval
+   software assumes that the scenarios are presented in random order
+   (i.e., they are not ordered by scenario number).
 
 Finding a candidate solution
 ----------------------------
@@ -40,10 +40,10 @@ Evaluating a candidate solution
 To evaluate a candidate solution with some scenarios, one might
 create a ``Xhat_Eval`` object and call its ``evaluate`` method 
 (resp. ``evaluate_one`` for a single scenario). It takes as
-an argument ``xhats``, a dictionnary of noon-anticipative policies for all 
+an argument ``xhats``, a dictionary of noon-anticipative policies for all 
 non-leaf nodes of a scenario tree. While for a 2-stage problem, ``xhats`` is
 just the candidate solution ``xhat_one``, for multistage problem the 
-dictionnary can be computed using the function ``walking_tree_xhats`` 
+dictionary can be computed using the function ``walking_tree_xhats`` 
 (resp. ``feasible_solution``).
 
 
@@ -68,7 +68,7 @@ Using stand alone ``mmw_conf.py``
 
 To use the stand along program a model compatible with ``Amalgomator`` and ``.npy`` file with a candidate solution to an instance of the model are required.
 
-First, ensure that the model to be used is compatable with the
+First, ensure that the model to be used is compatible with the
 ``Amalgomator`` class. This requires the model to have each of the
 following: a ``scenario_names_creator``, a ``scenario_creator``, an
 ``inparser_adder``, and a ``kw_creator``. See ``afarmer.py`` in
@@ -82,7 +82,7 @@ Once a model satisfies the requirement for amalgomator, next a ``.npy`` file sho
 ``sputils.ef_ROOT_nonants_npy_serializer(ama_object.ef, "xhat.npy")`` to your existing program (see the example in ``afarmer.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 corresponing path should be called on the command line.
+``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.
 
 Additional solver options can be specified with the ``--solver-options`` option.
 
@@ -95,5 +95,3 @@ realizations of the objective function at the candidate solution to
 construct an additional confidence interval about the mean of the
 realizations computed.
 
-
-

doc/src/drivers.rst

@@ -7,7 +7,7 @@ To make use of the hub and spoke system, you must come up with a
 driver that instantiates the objects and calls them. Nearly the
 last step in most drivers is a call to ``mpisppy.utils.sputils.spin_the_wheel``
 that calls the hub and spokes. Many of the example drivers take
-advanatage of shared code in the ``mpisppy.utils`` directory.
+advantage of shared code in the ``mpisppy.utils`` directory.
 
 We now proceed to consider driver examples that do not use the shared
 services provided with the examples; however, many developers will prefer
@@ -25,7 +25,7 @@ are
   creation. These dictionaries are ultimately fed to
   ``sputils.spin_the_wheel``.
 
-The contructors for the vanilla spokes take arguments that vary slightly depending
+The constructors for the vanilla spokes take arguments that vary slightly depending
 on the spoke, but all want the args passed in by the args parser,
 followed by ``scenario_creator`` function, a ``scenario_denoument`` function
 (that can be ``None``), a list of scenario names as ``all_scenario_names``,
@@ -83,8 +83,8 @@ they like. The disadvantage is that developers who use ``mpi-sppy``
 cannot count on it to detect spelling errors in options names.
 
 
-Contrasting ``_mpisppy_node_list` and ``all_node_names``
---------------------------------------------------------
+Contrasting ``_mpisppy_node_list`` and ``all_node_names``
+---------------------------------------------------------
 
 Note that ``_mpisppy_node_list``that is attached to scenarios does not have the leaf nodes, but ``all_node_names``
 that is part of the hub dictionary does.

doc/src/index.rst

@@ -18,6 +18,7 @@ MPI is used.
    hubs.rst
    spokes.rst
    extensions.rst
+   zhat.rst
    confidence_intervals.rst
    seqsamp.rst
    pysp.rst

doc/src/seqsamp.rst

@@ -3,19 +3,21 @@
 Sequential sampling
 ===================
 
-Similarly, given an confidence interval, one can try to find a candidate solution
+Given a confidence interval, one can try to find a candidate solution
 ``xhat_one`` such that its optimality gap has this confidence interval.
-The class ``SeqSampling`` implements three procedures described in 
+The class ``SeqSampling`` in ``mpisppy.confiden_intervals.seqsampling.py`` implements three procedures described in 
 [bm2011]_ and [bpl2012]_. It takes as an input a method to generate
-candidate solutions and options, and returns a ``xhat_one`` and a confidence interval on its optimality gap.
+candidate solutions and options, and its ``run`` method returns a ``xhat_one`` and a confidence interval on its optimality gap.
 
 There are two stopping criterion supported with names based on the initials of
-the authors: "BM" and "BPL".
+the authors who defined them: "BM" and "BPL".
 
 Examples of use with the ``farmer`` problem and several options can be found in the main of ``seqsampling.py``. The following options dictionaries are illustrated:
 
-- Relative Width
+- relative Width;
 
-- fixed width, sequential
+- fixed width, sequential;
 
 - fixed width with stochastic samples.
+
+  The keys used in the options dictionaries are taken directly from the corresponding paper, perhaps abbreviated in an obvious way. For example, the key `eps` corresponds to epsilon in the papers. 

doc/src/zhat.rst

@@ -0,0 +1,49 @@
+.. _zhat introduction:
+
+zhat and Background for Confidence Intervals
+============================================
+
+In this section, we provide some reference material useful for
+writing programs that compute confidence intervals.
+
+xhat
+----
+
+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.afarmer.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
+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
+`MMWConfidence Intervals`_ and
+`Sequential Sampling Confidence Intervals`_. Note: at the time of this writing, `zhat4xhat` does
+not support a starting scenario other than the first scenario, so
+some care might be needed if you want to avoid including scenarios
+used to compute xhat.
+
+
+seedoffset
+----------
+
+Most of the confidence interval code assumes that is can pass in a
+seed, particularly to the ``scenario_creator`` function so that
+replicates can be obtained. See ``examples.farmer.afarmer.py`` for an
+example.
+
+When only a small number of scenarios are available, the
+``scenario_creator`` function may need to take care to avoid
+attempting to access non-existent scenarios. If the data are provided
+in files, then `seedoffset` is a bit of a misnomer, but it needs to be
+added to scenarios numbers to get the scenario number. Note: to date,
+there is no code in the confidence interval code-base that rescales
+probabilities, so unlike the rest of `mpi-sppy`, present confidence
+interval code assumes equally likely scenarios.

examples/aircond/aaircond.py

@@ -1,6 +1,8 @@
 #for use with amalgomator mmw testing, so __main__ should be 2-stage for now,
 # and kw_creator also defaults to 2-stage
 #ReferenceModel for full set of scenarios for AirCond; June 2021
+# As of Oct 2021, aircond_submodels.py in the test directory has diverged
+# (it includes startups, but that needs to be handled with more elegance.)
 
 import pyomo.environ as pyo
 import numpy as np
@@ -12,6 +14,34 @@
 # Use this random stream:
 aircondstream = np.random.RandomState()
 
+def demands_creator(sname,branching_factors,start_seed, mudev, sigmadev, 
+                    starting_d=200,root_name="ROOT"):
+
+    max_d = 400
+    min_d = 0
+
+    if branching_factors is None:
+        raise RuntimeError("scenario_creator for aircond needs branching_factors")
+    scennum   = sputils.extract_num(sname)
+    # Find the right path and the associated seeds (one for each node) using scennum
+    prod = np.prod(branching_factors)
+    s = int(scennum % prod)
+    d = starting_d
+    demands = [d]
+    nodenames = [root_name]
+    for bf in branching_factors:
+        assert prod%bf == 0
+        prod = prod//bf
+        nodenames.append(str(s//prod))
+        s = s%prod
+
+    stagelist = [int(x) for x in nodenames[1:]]
+    for t in range(1,len(nodenames)):
+        aircondstream.seed(start_seed+sputils.node_idx(stagelist[:t],branching_factors))
+        d = min(max_d,max(min_d,d+aircondstream.normal(mudev,sigmadev)))
+        demands.append(d)
+
+    return demands,nodenames
 
 
 def StageModel_creator(time, demand, last_stage=False):
@@ -96,33 +126,54 @@ def total_cost_rule(m):
 
     return model
 
-def MakeNodesforScen(model,nodenames,branching_factors):
+
+def MakeNodesforScen(model,nodenames,branching_factors,starting_stage=1):
     #Create all nonleaf nodes used by the scenario
+    #Compatible with sample scenario creation
     TreeNodes = []
     for stage in model.T:
+
+        nonant_list=[model.stage_models[stage].RegularProd,
+                     model.stage_models[stage].OvertimeProd]
+        """
+        if model.start_ups:
+            nonant_list.append(model.stage_models[stage].StartUp)
+        """
         if stage ==1:
             ndn="ROOT"
             TreeNodes.append(scenario_tree.ScenarioNode(name=ndn,
                                                        cond_prob=1.0,
                                                        stage=stage,
                                                        cost_expression=model.stage_models[stage].StageObjective,
                                                        scen_name_list=None, # Not maintained
-                                                       nonant_list=[model.stage_models[stage].RegularProd,
-                                                                    model.stage_models[stage].OvertimeProd],
+                                                       nonant_list=nonant_list,
                                                        scen_model=model,
                                                        nonant_ef_suppl_list = [model.stage_models[stage].Inventory],
                                                        )
                              )
+        elif stage <=starting_stage:
+            parent_ndn = ndn
+            ndn = parent_ndn+"_0" #Only one node per stage before starting stage
+            TreeNodes.append(scenario_tree.ScenarioNode(name=ndn,
+                                                       cond_prob=1.0,
+                                                       stage=stage,
+                                                       cost_expression=model.stage_models[stage].StageObjective,
+                                                       scen_name_list=None, # Not maintained
+                                                       nonant_list=nonant_list,
+                                                       scen_model=model,
+                                                       nonant_ef_suppl_list = [model.stage_models[stage].Inventory],
+                                                       parent_name = parent_ndn
+                                                       )
+                             )
         elif stage < max(model.T): #We don't add the leaf node
             parent_ndn = ndn
-            ndn = parent_ndn+"_"+nodenames[stage-1]
+            ndn = parent_ndn+"_"+nodenames[stage-starting_stage]
             TreeNodes.append(scenario_tree.ScenarioNode(name=ndn,
-                                                       cond_prob=1.0/branching_factors[stage-2],
+                                                       cond_prob=1.0/branching_factors[stage-starting_stage-1],
                                                        stage=stage,
                                                        cost_expression=model.stage_models[stage].StageObjective,
                                                        scen_name_list=None, # Not maintained
-                                                       nonant_list=[model.stage_models[stage].RegularProd,
-                                                                    model.stage_models[stage].OvertimeProd],
+                                                       nonant_list=nonant_list,
                                                        scen_model=model,
                                                        nonant_ef_suppl_list = [model.stage_models[stage].Inventory],
                                                        parent_name = parent_ndn
@@ -131,8 +182,10 @@ def MakeNodesforScen(model,nodenames,branching_factors):
     return(TreeNodes)
 
 
-def scenario_creator(sname, BFs, num_scens=None, mudev=0, sigmadev=40, start_seed=0):
+def scenario_creator(sname, branching_factors, num_scens=None,
+                     mudev=0, sigmadev=40, start_seed=0):
     scennum   = sputils.extract_num(sname)
+    BFs = branching_factors
     # Find the right path and the associated seeds (one for each node) using scennum
     prod = np.prod(BFs)
     s = int(scennum % prod)
@@ -163,6 +216,56 @@ def scenario_creator(sname, BFs, num_scens=None, mudev=0, sigmadev=40, start_see
     return(model)
 
 
+#=========
+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.
+    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 randome 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
+    """
+    # Finding demands from stage 1 to t
+    if given_scenario is None:
+        if stage == 1:
+            past_demands = [200]
+        else:
+            raise RuntimeError(f"sample_tree_scen_creator for aircond needs a 'given_scenario' argument if the starting stage is greater than 1")
+    else:
+        past_demands = [given_scenario.stage_models[t].Demand for t in given_scenario.T if t<=stage]
+    optional_things = ['mudev','sigmadev','start_ups']
+    default_values = [0,40,False]
+    for thing,value in zip(optional_things,default_values):
+        if thing not in scenario_creator_kwargs:
+            scenario_creator_kwargs[thing] = value
+
+    #Finding demands for stages after t
+    future_demands,nodenames = demands_creator(sname, sample_branching_factors, 
+                                               start_seed = seed, 
+                                               mudev = scenario_creator_kwargs['mudev'], 
+                                               sigmadev = scenario_creator_kwargs['sigmadev'],
+                                               starting_d=past_demands[stage-1],
+                                               root_name='ROOT'+'_0'*(stage-1))
+
+    demands = past_demands+future_demands[1:] #The demand at the starting stage is in both past and future demands
+
+    model = aircond_model_creator(demands)
+
+    model._mpisppy_probability = 1/np.prod(sample_branching_factors)
+
+    #Constructing the nodes used by the scenario
+    model._mpisppy_node_list = MakeNodesforScen(model, nodenames, sample_branching_factors,
+                                                starting_stage=stage)
+
+    return model
+
+
 #=========
 def scenario_names_creator(num_scens,start=None):
     # (only for Amalgomator): return the full list of num_scens scenario names
@@ -194,7 +297,7 @@ def inparser_adder(inparser):
 #=========
 def kw_creator(options):
 
-    def _kw_arg(option_name, default = None, arg_name=None):
+    def _kwarg(option_name, default = None, arg_name=None):
         # options trumps args
         retval = options.get(option_name)
         if retval is not None:
@@ -212,7 +315,7 @@ def _kw_arg(option_name, default = None, arg_name=None):
     sigmadev = _kwarg("sigmadev", 40.)
     start_seed = _kwarg("start_seed", 1134)
     kwargs = {"num_scens" : options['num_scens'] if 'num_scens' in options else None,
-              "BFs": BFs,
+              "branching_factors": BFs,
               "mudev": mudev,
               "sigmadev": sigmadev,
               "start_seed": start_seed,