diff --git a/examples/farmer/bad_farmer.py b/examples/farmer/bad_farmer.py new file mode 100644 index 00000000..55777487 --- /dev/null +++ b/examples/farmer/bad_farmer.py @@ -0,0 +1,297 @@ +# 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.scenario_tree as scenario_tree +import mpisppy.utils.sputils as sputils +from mpisppy.utils import config + +# 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 = [var for var in model.DevotedAcreage.values()] + if scennum == 2: + varlist.reverse() + 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 + 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/bad_farmer_cylinders.py b/examples/farmer/bad_farmer_cylinders.py new file mode 100644 index 00000000..9129362a --- /dev/null +++ b/examples/farmer/bad_farmer_cylinders.py @@ -0,0 +1,161 @@ +# Copyright 2020 by B. Knueven, D. Mildebrath, C. Muir, J-P Watson, and D.L. Woodruff +# This software is distributed under the 3-clause BSD License. +# general example driver for farmer with cylinders + +import bad_farmer as farmer +import mpisppy.cylinders + +# Make it all go +from mpisppy.spin_the_wheel import WheelSpinner +import mpisppy.utils.sputils as sputils + +from mpisppy.utils import config +import mpisppy.utils.cfg_vanilla as vanilla + +from mpisppy.extensions.norm_rho_updater import NormRhoUpdater +from mpisppy.convergers.norm_rho_converger import NormRhoConverger +from mpisppy.convergers.primal_dual_converger import PrimalDualConverger +from mpisppy.utils.cfg_vanilla import extension_adder + +write_solution = True + +def _parse_args(): + # create a config object and parse + cfg = config.Config() + + cfg.num_scens_required() + cfg.popular_args() + cfg.two_sided_args() + cfg.ph_args() + cfg.aph_args() + cfg.xhatlooper_args() + cfg.fwph_args() + cfg.lagrangian_args() + cfg.lagranger_args() + cfg.xhatshuffle_args() + cfg.converger_args() + cfg.wxbar_read_write_args() + cfg.tracking_args() + cfg.add_to_config("crops_mult", + description="There will be 3x this many crops (default 1)", + domain=int, + default=1) + cfg.add_to_config("use_norm_rho_updater", + description="Use the norm rho updater extension", + domain=bool, + default=False) + cfg.add_to_config("run_async", + description="Run with async projective hedging instead of progressive hedging", + domain=bool, + default=False) + + cfg.parse_command_line("farmer_cylinders") + return cfg + + +def main(): + + cfg = _parse_args() + + num_scen = cfg.num_scens + crops_multiplier = cfg.crops_mult + + rho_setter = farmer._rho_setter if hasattr(farmer, '_rho_setter') else None + if cfg.default_rho is None and rho_setter is None: + raise RuntimeError("No rho_setter so a default must be specified via --default-rho") + + if cfg.use_norm_rho_converger: + if not cfg.use_norm_rho_updater: + raise RuntimeError("--use-norm-rho-converger requires --use-norm-rho-updater") + else: + ph_converger = NormRhoConverger + elif cfg.primal_dual_converger: + ph_converger = PrimalDualConverger + else: + ph_converger = None + + scenario_creator = farmer.scenario_creator + scenario_denouement = farmer.scenario_denouement + all_scenario_names = ['scen{}'.format(sn) for sn in range(num_scen)] + scenario_creator_kwargs = { + 'use_integer': False, + "crops_multiplier": crops_multiplier, + } + scenario_names = [f"Scenario{i+1}" for i in range(num_scen)] + + # Things needed for vanilla cylinders + beans = (cfg, scenario_creator, scenario_denouement, all_scenario_names) + + if cfg.run_async: + # Vanilla APH hub + hub_dict = vanilla.aph_hub(*beans, + scenario_creator_kwargs=scenario_creator_kwargs, + ph_extensions=None, + rho_setter = rho_setter) + else: + # Vanilla PH hub + hub_dict = vanilla.ph_hub(*beans, + scenario_creator_kwargs=scenario_creator_kwargs, + ph_extensions=None, + ph_converger=ph_converger, + rho_setter = rho_setter) + + if cfg.primal_dual_converger: + hub_dict['opt_kwargs']['options']\ + ['primal_dual_converger_options'] = { + 'verbose': True, + 'tol': cfg.primal_dual_converger_tol, + 'tracking': True} + + ## hack in adaptive rho + if cfg.use_norm_rho_updater: + extension_adder(hub_dict, NormRhoUpdater) + hub_dict['opt_kwargs']['options']['norm_rho_options'] = {'verbose': True} + + # FWPH spoke + if cfg.fwph: + fw_spoke = vanilla.fwph_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs) + + # Standard Lagrangian bound spoke + if cfg.lagrangian: + lagrangian_spoke = vanilla.lagrangian_spoke(*beans, + scenario_creator_kwargs=scenario_creator_kwargs, + rho_setter = rho_setter) + + # Special Lagranger bound spoke + if cfg.lagranger: + lagranger_spoke = vanilla.lagranger_spoke(*beans, + scenario_creator_kwargs=scenario_creator_kwargs, + rho_setter = rho_setter) + + # xhat looper bound spoke + if cfg.xhatlooper: + xhatlooper_spoke = vanilla.xhatlooper_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs) + + # xhat shuffle bound spoke + if cfg.xhatshuffle: + xhatshuffle_spoke = vanilla.xhatshuffle_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs) + + list_of_spoke_dict = list() + if cfg.fwph: + list_of_spoke_dict.append(fw_spoke) + if cfg.lagrangian: + list_of_spoke_dict.append(lagrangian_spoke) + if cfg.lagranger: + list_of_spoke_dict.append(lagranger_spoke) + if cfg.xhatlooper: + list_of_spoke_dict.append(xhatlooper_spoke) + if cfg.xhatshuffle: + list_of_spoke_dict.append(xhatshuffle_spoke) + + wheel = WheelSpinner(hub_dict, list_of_spoke_dict) + wheel.spin() + + if write_solution: + wheel.write_first_stage_solution('farmer_plant.csv') + wheel.write_first_stage_solution('farmer_cyl_nonants.npy', + first_stage_solution_writer=sputils.first_stage_nonant_npy_serializer) + wheel.write_tree_solution('farmer_full_solution') + +if __name__ == "__main__": + main() diff --git a/examples/run_all.py b/examples/run_all.py index f47aabc1..5548c02f 100644 --- a/examples/run_all.py +++ b/examples/run_all.py @@ -228,7 +228,7 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring): do_one("farmer", "farmer_cylinders.py", - 2, + 4, 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", diff --git a/mpisppy/spbase.py b/mpisppy/spbase.py index 3516e83e..6b5daafe 100644 --- a/mpisppy/spbase.py +++ b/mpisppy/spbase.py @@ -13,6 +13,7 @@ import logging import weakref import numpy as np +import array import pyomo.environ as pyo import mpisppy.utils.sputils as sputils from mpisppy import global_toc @@ -111,7 +112,7 @@ def __init__( self._attach_nonant_indices() self._attach_varid_to_nonant_index() self._create_communicators() - self._verify_nonant_lengths() + self._verify_nonants() self._set_sense() self._use_variable_probability_setter() @@ -151,7 +152,7 @@ def _set_sense(self, comm=None): "model sense (minimize or maximize)" ) - def _verify_nonant_lengths(self): + def _verify_nonants(self): local_node_nonant_lengths = {} # keys are tree node names # we need to accumulate all local contributions before the reduce @@ -177,7 +178,37 @@ def _verify_nonant_lengths(self): if val != int(max_val[0]): raise RuntimeError(f"Tree node {ndn} has scenarios with different numbers of non-anticipative " f"variables: {val} vs. max {max_val[0]}") - + + # compare node names (reduction) + local_nonant_char_array = {} + for k,s in self.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + for node in s._mpisppy_node_list: + ndn = node.name + if ndn not in local_nonant_char_array: + local_nonant_char_array[ndn] = {} + for i, var in enumerate(node.nonant_vardata_list): + if i in local_nonant_char_array[ndn]: + if var.name != local_nonant_char_array[ndn][i]: + raise RuntimeError(f"[rank {self.global_rank}] Tree node {ndn} has different non-anticipative " + f"variables in position {i}, scenario {s} has name " + f"{var.name}, some other scenario has name " + f"{local_nonant_char_array[ndn][i]}") + else: + local_nonant_char_array[ndn][i] = var.name + + for ndn, var_names in local_nonant_char_array.items(): + a = array.array("u") + for i, vname in var_names.items(): + a.extend(vname) + localnames = np.array( array.array("i", a.tobytes()) ) + globalnames = np.zeros( len(localnames), "i" ) + self.comms[ndn].Allreduce([localnames, MPI.INT], + [globalnames, MPI.INT], + op=MPI.MAX) + if not np.all(localnames == globalnames): + raise RuntimeError(f"For node {ndn}, non-anticipative variable lists do not match! ") + def _check_nodenames(self): for ndn in self.all_nodenames: if ndn != 'ROOT' and sputils.parent_ndn(ndn) not in self.all_nodenames: