diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index 73ad89d0..3a88edda 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -245,7 +245,7 @@ jobs: pip install -e . - name: run tests - timeout-minutes: 100 + timeout-minutes: 10 run: | cd mpisppy/tests # envall does nothing @@ -274,7 +274,7 @@ jobs: pip install -e . - name: run tests - timeout-minutes: 100 + timeout-minutes: 10 run: | cd mpisppy/tests # envall does nothing @@ -290,7 +290,7 @@ jobs: - uses: conda-incubator/setup-miniconda@v2 with: activate-environment: test_env - python-version: 3.8 + python-version: 3.9 auto-activate-base: false - name: Install dependencies run: | @@ -318,7 +318,7 @@ jobs: - uses: conda-incubator/setup-miniconda@v2 with: activate-environment: test_env - python-version: 3.8 + python-version: 3.9 auto-activate-base: false - name: Install dependencies run: | @@ -416,7 +416,7 @@ jobs: - uses: conda-incubator/setup-miniconda@v2 with: activate-environment: test_env - python-version: 3.8 + python-version: 3.9 auto-activate-base: false - name: Install dependencies run: | @@ -455,14 +455,14 @@ jobs: pip install -e . - name: run pysp model tests - timeout-minutes: 100 + timeout-minutes: 10 run: | cd mpisppy/tests # envall does nothing python test_pysp_model.py - name: run pysp unit tests - timeout-minutes: 100 + timeout-minutes: 10 run: | cd mpisppy/utils/pysp_model pytest -v . @@ -477,7 +477,7 @@ jobs: - uses: conda-incubator/setup-miniconda@v2 with: activate-environment: test_env - python-version: 3.8 + python-version: 3.9 auto-activate-base: false - name: Install dependencies run: | @@ -493,3 +493,41 @@ jobs: run: | cd mpisppy/tests mpiexec -np 2 python -m mpi4py test_with_cylinders.py + + test-agnostic: + name: tests on agnostic + runs-on: ubuntu-latest + needs: [ruff] + + steps: + - uses: actions/checkout@v3 + - uses: conda-incubator/setup-miniconda@v2 + with: + activate-environment: test_env + python-version: 3.11 + auto-activate-base: false + - name: Install dependencies + run: | + conda install mpi4py pandas setuptools + pip install pyomo xpress cplex + pip install numpy + python -m pip install amplpy --upgrade + python -m amplpy.modules install highs cbc gurobi + python -m pip install gamspy + # license? + + - name: setup the program + run: | + pip install -e . + + - name: run agnostic tests + timeout-minutes: 10 + run: | + cd mpisppy/tests + python test_agnostic.py + + - name: run agnostic cylinders + timeout-minutes: 10 + run: | + cd mpisppy/agnostic/examples + python afew_agnostic.py --oversubscribe diff --git a/agnostic.txt b/agnostic.txt new file mode 100644 index 00000000..417e6061 --- /dev/null +++ b/agnostic.txt @@ -0,0 +1,138 @@ +newest notes on top +---------------------------------------- + +host means mpisppy +guest means whatever the guest language is (which can even be Pyomo, of course) + + +Jan 3 2024: We are going to require that the guest takes care of + bundling and presents the host with "proper" bundles + that look to the host like regular scenarios. + We have AMPL and Pyomo working as guests and GAMS seems to work OK, + except that there are file contention issues that can probably be solved + easily (I think they have an example with parallel execution) + +------------------------- + + +Aug 27 + + +The example is now farmer_agnostic.py and agnostic_cylinders.py (ag.bash) +in the farmer example directory. + +HEY::: who updates W? Let's do it right before the solve so we don't have to care... + this is a little dangerous at present because we are being silent when + they are not there because of xhatters + + +== nonant communication needs to be on the host side, of course, but + then callouts need to get values copied to/from the guest + +== we need to enforce the assumption that the nonants match in every way possible + between host and guest at all times (e.g. _restore_nonants needs to callout, + but _save_nonants does not. + +== bundling is left hanging in dangerous ways: e.g. Eobjective in two places + +Aug 26 + +Work proceeds on three fronts: + +- Addiing callouts to mpisppy +- creating the example guest language functions (in Pyomo) +- creating examples of the guest language funtions in other languages + += note: the callouts are almost always inside the local scenario loop in mpisspy + += note: after the solve we are going update the host model x values so +it can update the xbars + += The host will update the w's then a callout will send the new values to + the guest + += I am not even thinking about extensions... + += For bundling, we need to be able to solve EFs (to state the obvious) + +circa Aug 25 + +- no changes needed in spbase so long as agnostic scenario creator is passed in + +- EF will be basically a rewrite because we use blocks + +- ? for the phbase constructor: pass cfg that contains Ag or pass Ag?; use the options dict for now + +- working on an example, which is farmer_agnostic.py run from the __main__ of agnostic.py + (farmer_agnostic is presently run from the utils directory and I have a copy of farmer.py there as well) + + +=============================================================== +Thoughts about AML-agnostic extensions to mpi-sppy +(started by DLW 18 Dec 2022) +Just thinking about support for straight PH for now. Bundles are probably the first thing to add. + +-1. BTW: both GAMS and AMPL have python APIs. + +0. On the Python side each scenario is still a Pyomo "model" (that perhaps has only the nonant Vars) +with _mpisppy_data and _mpisppy_model attached. + - it might result in a little speed-up to cut Pyomo out some day, but we should start with Pyomo, I think + +1. Some functions in spbase, most functions in phbase, and maybe some functions in spopt will call this function: + + def callout_agnostic(cfg, name, **kwargs): + """ callout for AML-agnostic support + Args: + cfg (Config): the field "AML_agnostic" might contain a module with callouts + name (str): the function name to call + kwargs (dict): the keywords args for the callout function + Calls: + a callout function that presumably has side-effects + Returns: + True if the callout was done and False if not + """ + + if cfg.get(AML_agnostic, ifmissing=None) is not None: + if not hasattr(cfg.AML_agnostic, name): + raise RuntimeError(f"AML-agnostic module is missing function {name}") + fct = getattr(cfg.AML_agnostic, name) + fct(**kwargs) + return True + else: + return False + + The function that is called by fct(**kwargs) will do the work of + interacting with the AML and updating mpi-sppy structures that are + passed in as kwargs. Note that cfg is going to need to be attached + some some objects that don't presently have it (e.g. SPOpt). Some + functions in mpi-sppy will want to return immediately if + callout_agnostic returns True (i.e., have the callout function do + all the work). + + +2 in spbase: + - don't support variable_prob + - _set_scense needs a callout + +3 The scenario_creator function will be in Python and needs to call an AML scenario creator function + +4. In solve_one it might be easiest to just do everything in the callout function. Maybe that will be the +case for many callouts. But from a maintenance perspective, it would best to have mpi-sppy code +do as much as possible and the callout Python do as little as possible. + +5. Think about Compute_Xbar and Update_W. They probably need to do all their processing then do +a callout at the end so the AML model can be updated. + +======================================= +Notes about callouts + +- There are going to be a lot of them! + +- Maybe it would be better to drop the name argument and use +inspect.stack()[1][3] in callout_agnostic so the name of the cfg field +is the name of the calling function in mpi-sppy. + This would + o standardize names, + o eliminate some cut-and-paste errors, + o make it a little wierd if there were ever two callouts from one function, + but not that wierd (just handled with kwarg flag) and two callouts should be rare. diff --git a/doc/src/agnostic.rst b/doc/src/agnostic.rst new file mode 100644 index 00000000..ea397a7d --- /dev/null +++ b/doc/src/agnostic.rst @@ -0,0 +1,171 @@ +AML Agnosticism +=============== + +The mpi-sppy package provides callouts so that algebraic modeling languages +(AMLs) other than Pyomo can be used. A growing number of AMLs are supported +as `guest` languages (we refer to mpi-sppy as the `host`). This code is +in an alpha-release state; use with extreme caution. + +From the end-user's perspective +------------------------------- + +When mpi-sppy is used for a model developed in an AML for which support +has been added, the end-user runs the ``mpisppy.agnostic.agnostic_cylinders.py`` +program which serves as a driver that takes command line arguments and +launches the requested cylinders. The file +``mpisppy.agnostic.go.bash`` provides examples of a few command lines. + + +From the modeler's perspective +------------------------------ + +Assuming support has been added for the desired AML, the modeler supplies +two files: + +- a model file with the model written in the guest AML (AMPL example: ``mpisppy.agnostic.examples.farmer.mod``) +- a thin model wrapper for the model file written in Python (AMPL example: ``mpisppy.agnostic.examples.farmer_ampl_model.py``). This thin python wrapper is model specific. + +There can be a little confusion if there are error messages because +both files are sometimes refered to as the `model file.` + +Most modelers will probably want to import the deterministic guest model into their +python wrapper for the model and the scenario_creator function in the wrapper +modifies the stochastic paramaters to have values that depend on the scenario +name argument to the scenario_creator function. + +(An exception is when the guest is in Pyomo, then the wrapper +file might as well contain the model specification as well so +there typically is only one file. However, there is not particularly +good reason to use the agnostic machinery for a Pyomo model.) + + +From the developers perspective +------------------------------- + +If support has not yet been added for an AML, it is almost easier to +add support than to write a guest interface for a particular model. To +add support for a language, you need to write a general guest +interface in Python for it (see, e.g., ampl_guest.py or +pyomo_guest.py) and you need to add/edit a few lines in +``mpisppy.agnostic.agnostic_cylinders.py`` to allow end-users to +access it. + + +Special Note for developers +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The general-purpose guest interfaces might not be the fastest possible +for many guest languages because they don't use indexes from the +original model when updating the objective function. If this is an issue, +you might want to write a problem-specific module to replace the guest +interface and the model wrapper with a single module. For an example, see +``examples.farmer.agnostic.farmer_xxxx_agnostic``, where xxxx is replaced, +e.g., by ampl. + +Architecture +^^^^^^^^^^^^ +The following picture presents the architecture of the files. + +.. image:: images/agnostic_architecture.png + :alt: Architecture of the agnostic files + :width: 700px + :align: center + +We note "xxxx" the specific problem, for instance farmer. We note "yyyy" the guest language, for instance "ampl". +Two methods are presented. Either a method specific to the problem, or a generic method. +Regardless of the method, the file ``agnostic.py`` and ``xxxx.yyyy`` need to be used. +``agnostic.py`` is already implemented and must not be modified as all the files presented above the line "developer". +``xxxx.yyyy`` is the model in the guest language and must be given by the modeler such as all the files under the line "modeler". + +The files ``agnostic_yyyy_cylinders.py`` and ``agnostic_cylinders.py`` are equivalent. +The file ``xxxx_yyyy_agnostic.py`` for the specific case is split into ``yyyy_guest.py`` and ``xxxx_yyyy_model.py`` for the generic case. + + +It is worth noting that the scenario creator is defined in 3 files. +It is first defined in the file specific to the problem and the guest language ``xxxx_yyyy_model.py``. At this point it may not return a scenario. +It is then wrapped in a file only specific to the language ``yyyy_guest.py``. At chich point it returns the dictionary ``gd`` which indludes the scenario. +Finally the tree structure is attached in ``agnostic.py``. + + +Bundles +------- + +The use of scenario bundles can dramatically improve the performance +of scenario decomposition algorithms such as PH and APH. Although mpi-sppy +has facitilites for forming bundles, the mpi-sppy +``agnostic`` package assumes that bundles will be completely handled +by the guest. Bundles will be returned by the scenario creator function +as if they are a scenario. Although it seems sort of like a trick, it is +really the way bundles are intended to operate so we sometimes refer to +`true` bundles, which are used in non-agnostic way as briefly +described in section :ref:`Pickled-Bundles`. + +Overview of Recommended Bundle Practices +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Modify the scenario creator function so that if the scenario name +starts with the string "scen" it returns a single scenario, but if the +name starts with "bundle" it returns the full extensive formulation for +a group of scenarios (i.e. a bundle). We typically number scenarios +and the scenario or bundle number is at the end of the first +positional argument for the scenario creator function (i.e. at +the end of the scenario name). + +If the name starts with bundle, the scenario creator function can call +itself with the proper list of scenario names to get the scenarios +to form the EF that will be returned. We recommend names for +bundles such as "bundle_xxx_yyy" where xxx and yyy give the +first and last scenario number in the bundle. +You could also pass in a dictionary that maps bundle numbers to lists of +scenario numbers as a keyword argument to the scenario_creator function +and then append the bundle number to "bundle" and pass it as the positional +scenario name argument to the scenario creator function. + +Some notes +^^^^^^^^^^ + +- The helper function called ``scenario_names_creator`` needs to be co-opted +to instead create bundle names and the code in the scenario_creator function +then needs to create its own scenario names for bundles. At the time +of this writing this results in a major hack being needed in order to +get bundle information to the names creator in the Pyomo example described +below. You need to supply a function called ``bundle_hack`` in your python model file that +does whatever needs to be done to alert the names creator that there +bundles. The function takes the config object as an argument. +See ``mpisppy.agnostic.farmer4agnostic.py`` +- There is a heavy bias toward uniform probabilities in the examples and in + the mpi-sppy utilities. Scenario probabilities are attached to the scenario + as ``_mpisppy_probability`` so if your probabilities are not uniform, you will + need to calculate them for each bundle (your EF maker code can do that for you). Note that even if probabilities are uniform for the scenarios, they won't + be uniform for the bundles unless you require that the bundle size divides + the number of scenarios. +- There is a similar bias toward two stage problems, which is + extreme for the agnostic package. If you have a multi-stage + problem, you can make things a lot easier for yourself if you require + that the bundles contain all scenarios emanating from each second stage node + (e.g., on bundle per some integer number of second stage nodes). This + is what is done in (non-agnostic) :ref:`Pickled-Bundles`. The result of this + is that your multi-stage problem will look like a two-stage problem to + mpi-sppy. + +Example +^^^^^^^ + +The example ``mpisppy.agnostic.farmer4agnostic.py`` contains example code. + +.. Note:: + In order to get information from the command line about bundles into the + ``scenario_names_creator`` the ``bundle_hack`` function is called + called by the cylinders driver program very early. For this example, + function sets global variables called ``bunsize`` and ``numbuns``. + +The script ``mpisppy.agnostic.examples.go.bash`` runs the example (and maybe some +other examples). + + +Notes about Gurobipy +-------------------- + +The current implementation of gurobipy assumes that nonants that are in +the objective function appear direclty there (not via some other +variable constrained in some way to represent them). diff --git a/doc/src/images/agnostic_architecture.png b/doc/src/images/agnostic_architecture.png new file mode 100644 index 00000000..055ccaff Binary files /dev/null and b/doc/src/images/agnostic_architecture.png differ diff --git a/doc/src/images/agnostic_architecture.pptx b/doc/src/images/agnostic_architecture.pptx new file mode 100644 index 00000000..faeeb6b2 Binary files /dev/null and b/doc/src/images/agnostic_architecture.pptx differ diff --git a/doc/src/index.rst b/doc/src/index.rst index 4b9ce74e..abfcdb5f 100644 --- a/doc/src/index.rst +++ b/doc/src/index.rst @@ -30,6 +30,7 @@ MPI is used. properbundles.rst grad_rho.rst w_rho.rst + agnostic.rst admmWrapper.rst stoch_admmWrapper.rst helper_functions.rst diff --git a/examples/farmer/AMPL/README.txt b/examples/farmer/AMPL/README.txt new file mode 100644 index 00000000..9a611df9 --- /dev/null +++ b/examples/farmer/AMPL/README.txt @@ -0,0 +1 @@ +The directory has a few miscellaneous files related to getting started with AMPL. diff --git a/examples/farmer/AMPL/farmer.run b/examples/farmer/AMPL/farmer.run new file mode 100644 index 00000000..d2d95aa4 --- /dev/null +++ b/examples/farmer/AMPL/farmer.run @@ -0,0 +1,5 @@ +# include farmer.run +model farmer_test.ampl; +option solver xpress; +solve; + diff --git a/examples/farmer/AMPL/farmer_stochastic.ampl b/examples/farmer/AMPL/farmer_stochastic.ampl new file mode 100644 index 00000000..174660b5 --- /dev/null +++ b/examples/farmer/AMPL/farmer_stochastic.ampl @@ -0,0 +1,98 @@ +# The farmer's problem in AMPL +# +# Reference: +# John R. Birge and Francois Louveaux. Introduction to Stochastic Programming. +# +# AMPL coding by Victor Zverovich. + +function expectation; +function random; + +suffix stage IN; + +set Crops; + +set Scen; +param P{Scen}; # probabilities + +param TotalArea; # acre +param PlantingCost{Crops}; # $/acre +param SellingPrice{Crops}; # $/T +param ExcessSellingPrice; # $/T +param PurchasePrice{Crops}; # $/T +param MinRequirement{Crops}; # T +param BeetsQuota; # T + +# Area in acres devoted to crop c. +var area{c in Crops} >= 0; + +# Tons of crop c sold (at favourable price) under scenario s. +var sell{c in Crops} >= 0, suffix stage 2; + +# Tons of sugar beets sold in excess of the quota under scenario s. +var sell_excess >= 0, suffix stage 2; + +# Tons of crop c bought under scenario s +var buy{c in Crops} >= 0, suffix stage 2; + +# The random variable (parameter) representing the yield of crop c. +var RandomYield{c in Crops}; + +# Realizations of the yield of crop c. +param Yield{c in Crops, s in Scen}; # T/acre + +maximize profit: + expectation( + ExcessSellingPrice * sell_excess + + sum{c in Crops} (SellingPrice[c] * sell[c] - + PurchasePrice[c] * buy[c])) - + sum{c in Crops} PlantingCost[c] * area[c]; + +s.t. totalArea: sum {c in Crops} area[c] <= TotalArea; + +s.t. requirement{c in Crops}: + RandomYield[c] * area[c] - sell[c] + buy[c] >= MinRequirement[c]; + +s.t. quota: sell['beets'] <= BeetsQuota; + +s.t. sellBeets: + sell['beets'] + sell_excess <= RandomYield['beets'] * area['beets']; + +yield: random({c in Crops} (RandomYield[c], {s in Scen} Yield[c, s])); + +data; + +set Crops := wheat corn beets; +set Scen := below average above; + +param TotalArea := 500; + +param Yield: + below average above := + wheat 2.0 2.5 3.0 + corn 2.4 3.0 3.6 + beets 16.0 20.0 24.0; + +param PlantingCost := + wheat 150 + corn 230 + beets 260; + +param SellingPrice := + wheat 170 + corn 150 + beets 36; + +param ExcessSellingPrice := 10; + +param PurchasePrice := + wheat 238 + corn 210 + beets 100; + +param MinRequirement := + wheat 200 + corn 240 + beets 0; + +param BeetsQuota := 6000; \ No newline at end of file diff --git a/examples/farmer/AMPL/install.txt b/examples/farmer/AMPL/install.txt new file mode 100644 index 00000000..5b27ec70 --- /dev/null +++ b/examples/farmer/AMPL/install.txt @@ -0,0 +1,7 @@ +0. follow instructions at https://pypi.org/project/amplpy +1. install asl: + get asl from coin-or tools third party ASL + $ ./get.asl + $ ./configure + $ make + $ make install diff --git a/examples/farmer/GAMS/farmer_augmented.gms b/examples/farmer/GAMS/farmer_augmented.gms new file mode 100644 index 00000000..9553e056 --- /dev/null +++ b/examples/farmer/GAMS/farmer_augmented.gms @@ -0,0 +1,101 @@ +$title The Farmer's Problem formulated for GAMS/DECIS (FARM,SEQ=199) + +$onText +This model helps a farmer to decide how to allocate +his or her land. The yields are uncertain. + + +Birge, R, and Louveaux, F V, Introduction to Stochastic Programming. +Springer, 1997. + +Keywords: linear programming, stochastic programming, agricultural cultivation, + farming, cropping +$offText + +*$if not set decisalg $set decisalg decism + +Set + crop / wheat, corn, sugarbeets / + cropr(crop) 'crops required for feeding cattle' / wheat, corn / + cropx / wheat + corn + beets1 'up to 6000 ton' + beets2 'in excess of 6000 ton' /; + +Parameter + yield(crop) 'tons per acre' / wheat 2.5 + corn 3 + sugarbeets 20 / + plantcost(crop) 'dollars per acre' / wheat 150 + corn 230 + sugarbeets 260 / + sellprice(cropx) 'dollars per ton' / wheat 170 + corn 150 + beets1 36 + beets2 10 / + purchprice(cropr) 'dollars per ton' / wheat 238 + corn 210 / + minreq(cropr) 'minimum requirements in ton' / wheat 200 + corn 240 / + ph_W(crop) 'ph weight' / wheat 0 + corn 0 + sugarbeets 0 / + xbar(crop) 'ph average' / wheat 0 + corn 0 + sugarbeets 0 / + rho(crop) 'ph rho' / wheat 0 + corn 0 + sugarbeets 0 /; + +Scalar + land 'available land' / 500 / + maxbeets1 'max allowed' / 6000 / + W_on 'activate w term' / 0 / + prox_on 'activate prox term' / 0 /; + +*-------------------------------------------------------------------------- +* First a non-stochastic version +*-------------------------------------------------------------------------- +Variable + x(crop) 'acres of land' + w(cropx) 'crops sold' + y(cropr) 'crops purchased' + yld(crop) 'yield' + negprofit 'objective variable'; + +Positive Variable x, w, y; + +Equation + profitdef 'objective function' + landuse 'capacity' + req(cropr) 'crop requirements for cattle feed' + ylddef 'calc yields' + beets 'total beet production'; + +$onText +The YLD variable and YLDDEF equation isolate the stochastic +YIELD parameter into one equation, making the DECIS setup +somewhat easier than if we would substitute YLD out of +the model. +$offText + +profitdef.. negprofit =e= + sum(crop, plantcost(crop)*x(crop)) + + sum(cropr, purchprice(cropr)*y(cropr)) + - sum(cropx, sellprice(cropx)*w(cropx)) + + W_on * sum(crop, ph_W(crop)*x(crop)) + + prox_on * sum(crop, rho(crop)*(x(crop) - xbar(crop))*(x(crop) - xbar(crop))); + +landuse.. sum(crop, x(crop)) =l= land; + +ylddef(crop).. yld(crop) =e= yield(crop)*x(crop); + +req(cropr).. yld(cropr) + y(cropr) - sum(sameas(cropx,cropr),w(cropx)) =g= minreq(cropr); + +beets.. w('beets1') + w('beets2') =l= yld('sugarbeets'); + +w.up('beets1') = maxbeets1; + +Model simple / profitdef, landuse, req, beets, ylddef /; + +Option QCP = Cplex; +solve simple using qcp minimizing negprofit; diff --git a/examples/farmer/GAMS/farmer_average.gms b/examples/farmer/GAMS/farmer_average.gms new file mode 100644 index 00000000..10af14a2 --- /dev/null +++ b/examples/farmer/GAMS/farmer_average.gms @@ -0,0 +1,91 @@ +$title The Farmer s Problem formulated for GAMS/DECIS (FARM,SEQ=199) + +$onText +This model helps a farmer to decide how to allocate +his or her land. The yields are uncertain. + + +Birge, R, and Louveaux, F V, Introduction to Stochastic Programming. +Springer, 1997. + +Keywords: linear programming, stochastic programming, agricultural cultivation, + farming, cropping +$offText + +*$if not set decisalg $set decisalg decism + +Set + crop / wheat, corn, sugarbeets / + cropr(crop) 'crops required for feeding cattle' / wheat, corn / + cropx / wheat + corn + beets1 'up to 6000 ton' + beets2 'in excess of 6000 ton' /; + +Parameter + yield(crop) 'tons per acre' / wheat 2.5 + corn 3 + sugarbeets 20 / + plantcost(crop) 'dollars per acre' / wheat 150 + corn 230 + sugarbeets 260 / + sellprice(cropx) 'dollars per ton' / wheat 170 + corn 150 + beets1 36 + beets2 10 / + purchprice(cropr) 'dollars per ton' / wheat 238 + corn 210 / + minreq(cropr) 'minimum requirements in ton' / wheat 200 + corn 240 /; + +Scalar + land 'available land' / 500 / + maxbeets1 'max allowed' / 6000 /; + +*-------------------------------------------------------------------------- +* First a non-stochastic version +*-------------------------------------------------------------------------- +Variable + x(crop) 'acres of land' + w(cropx) 'crops sold' + y(cropr) 'crops purchased' + yld(crop) 'yield' + profit 'objective variable'; + +Positive Variable x, w, y; + +Equation + profitdef 'objective function' + landuse 'capacity' + req(cropr) 'crop requirements for cattle feed' + ylddef 'calc yields' + beets 'total beet production'; + +$onText +The YLD variable and YLDDEF equation isolate the stochastic +YIELD parameter into one equation, making the DECIS setup +somewhat easier than if we would substitute YLD out of +the model. +$offText + +profitdef.. profit =e= - sum(crop, plantcost(crop)*x(crop)) + - sum(cropr, purchprice(cropr)*y(cropr)) + + sum(cropx, sellprice(cropx)*w(cropx)); + +landuse.. sum(crop, x(crop)) =l= land; + +ylddef(crop).. yld(crop) =e= yield(crop)*x(crop); + +req(cropr).. yld(cropr) + y(cropr) - sum(sameas(cropx,cropr),w(cropx)) =g= minreq(cropr); + +beets.. w('beets1') + w('beets2') =l= yld('sugarbeets'); + +x.up(crop) = land; +w.up('beets1') = maxbeets1; +$onText +__InsertPH__here_Model_defined_three_lines_later +$offText + +Model simple / profitdef, landuse, req, beets, ylddef /; + +solve simple using lp maximizing profit; \ No newline at end of file diff --git a/examples/farmer/GAMS/farmer_average.py b/examples/farmer/GAMS/farmer_average.py new file mode 100644 index 00000000..b53479e0 --- /dev/null +++ b/examples/farmer/GAMS/farmer_average.py @@ -0,0 +1,28 @@ +############################################################################### +# 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. +############################################################################### +# This is a program for testing and development +import os +import sys +import gams +import gamspy_base + +this_dir = os.path.dirname(os.path.abspath(__file__)) + +gamspy_base_dir = gamspy_base.__path__[0] + +w = gams.GamsWorkspace(working_directory=this_dir, system_directory=gamspy_base_dir) + +#model = w.add_job_from_file("farmer_average_ph.gms") +#model = w.add_job_from_file("farmer_augmented.gms") +#model = w.add_job_from_file("farmer_linear_augmented.gms") +#model = w.add_job_from_file("farmer_average_ph_quadratic") +model = w.add_job_from_file("farmer_average_ph_linearized") + + +model.run(output=sys.stdout) diff --git a/examples/farmer/GAMS/farmer_linear_augmented.gms b/examples/farmer/GAMS/farmer_linear_augmented.gms new file mode 100644 index 00000000..d085152e --- /dev/null +++ b/examples/farmer/GAMS/farmer_linear_augmented.gms @@ -0,0 +1,120 @@ +$title The Farmer's Problem formulated for GAMS/DECIS (FARM,SEQ=199) + +$onText +Linearize the prox term: + B is xbar + L is lower bound on x (which is zero for farmer) + U is upper bound on x (which is land=500 for farmer) + penalty >= B^2 - BL - Bx + Lx + penalty >= B^2 - BU - Bx + Ux + +This model helps a farmer to decide how to allocate +his or her land. The yields are uncertain. + + +Birge, R, and Louveaux, F V, Introduction to Stochastic Programming. +Springer, 1997. + +Keywords: linear programming, stochastic programming, agricultural cultivation, + farming, cropping +$offText + +*$if not set decisalg $set decisalg decism + +Set + crop / wheat, corn, sugarbeets / + cropr(crop) 'crops required for feeding cattle' / wheat, corn / + cropx / wheat + corn + beets1 'up to 6000 ton' + beets2 'in excess of 6000 ton' /; + +Parameter + yield(crop) 'tons per acre' / wheat 2.5 + corn 3 + sugarbeets 20 / + plantcost(crop) 'dollars per acre' / wheat 150 + corn 230 + sugarbeets 260 / + sellprice(cropx) 'dollars per ton' / wheat 170 + corn 150 + beets1 36 + beets2 10 / + purchprice(cropr) 'dollars per ton' / wheat 238 + corn 210 / + minreq(cropr) 'minimum requirements in ton' / wheat 200 + corn 240 / + ph_W(crop) 'ph weight' / wheat 0 + corn 0 + sugarbeets 0 / + xbar(crop) 'ph average' / wheat 0 + corn 0 + sugarbeets 0 / + rho(crop) 'ph rho' / wheat 0 + corn 0 + sugarbeets 0 /; + +Scalar + land 'available land' / 500 / + maxbeets1 'max allowed' / 6000 / + W_on 'activate w term' / 0 / + prox_on 'activate prox term' / 0 /; + +*-------------------------------------------------------------------------- +* First a non-stochastic version +*-------------------------------------------------------------------------- +Variable + x(crop) 'acres of land' + w(cropx) 'crops sold' + y(cropr) 'crops purchased' + yld(crop) 'yield' + PHpenalty(crop) 'linearized prox penalty' + negprofit 'objective variable'; + +Positive Variable x, w, y; + +Equation + profitdef 'objective function' + landuse 'capacity' + req(cropr) 'crop requirements for cattle feed' + ylddef 'calc yields' + PenLeft(crop) 'left side of linearized PH penalty' + PenRight(crop) 'right side of linearized PH penalty' + beets 'total beet production'; + +$onText +The YLD variable and YLDDEF equation isolate the stochastic +YIELD parameter into one equation, making the DECIS setup +somewhat easier than if we would substitute YLD out of +the model. +$offText + +profitdef.. negprofit =e= + sum(crop, plantcost(crop)*x(crop)) + + sum(cropr, purchprice(cropr)*y(cropr)) + - sum(cropx, sellprice(cropx)*w(cropx)) + + W_on * sum(crop, ph_W(crop)*x(crop)) + + prox_on * sum(crop, 0.5 * rho(crop) * PHpenalty(crop)); + +landuse.. sum(crop, x(crop)) =l= land; + +ylddef(crop).. yld(crop) =e= yield(crop)*x(crop); + +req(cropr).. yld(cropr) + y(cropr) - sum(sameas(cropx,cropr),w(cropx)) =g= minreq(cropr); + +beets.. w('beets1') + w('beets2') =l= yld('sugarbeets'); + +PenLeft(crop).. sqr(xbar(crop)) + xbar(crop)*0 + xbar(crop) * x(crop) + land * x(crop) =g= PHpenalty(crop); +PenRight(crop).. sqr(xbar(crop)) - xbar(crop)*land - xbar(crop)*x(crop) + land * x(crop) =g= PHpenalty(crop); + +PHpenalty.lo(crop) = 0; +PHpenalty.up(crop) = max(sqr(xbar(crop) - 0), sqr(land - xbar(crop))); + +w.up('beets1') = maxbeets1; + +x.lo(crop) = 0; +x.up(crop) = land; + +Model simple / profitdef, landuse, req, beets, ylddef, PenLeft, PenRight /; + +Option LP = Cplex; +solve simple using lp minimizing negprofit; diff --git a/examples/farmer/GAMS/farmer_linear_augmented.py b/examples/farmer/GAMS/farmer_linear_augmented.py new file mode 100644 index 00000000..f4d36647 --- /dev/null +++ b/examples/farmer/GAMS/farmer_linear_augmented.py @@ -0,0 +1,164 @@ +############################################################################### +# 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. +############################################################################### +# to assist with debugging (a lot of stuff here is really needed) +import os +import sys +import gams +import gamspy_base + +def create_model(scennum = 2): + # create a model and return a dictionary describing it + + this_dir = os.path.dirname(os.path.abspath(__file__)) + + gamspy_base_dir = gamspy_base.__path__[0] + + ws = gams.GamsWorkspace(working_directory=this_dir, system_directory=gamspy_base_dir) + + job = ws.add_job_from_file("farmer_linear_augmented.gms") + + job.run(output=sys.stdout) + + cp = ws.add_checkpoint() + mi = cp.add_modelinstance() + + job.run(checkpoint=cp) + + crop = mi.sync_db.add_set("crop", 1, "crop type") + + y = mi.sync_db.add_parameter_dc("yield", [crop,], "tons per acre") + + ph_W = mi.sync_db.add_parameter_dc("ph_W", [crop,], "ph weight") + xbar = mi.sync_db.add_parameter_dc("xbar", [crop,], "ph average") + rho = mi.sync_db.add_parameter_dc("rho", [crop,], "ph rho") + + W_on = mi.sync_db.add_parameter("W_on", 0, "activate w term") + prox_on = mi.sync_db.add_parameter("prox_on", 0, "activate prox term") + + mi.instantiate("simple min negprofit using lp", + [ + gams.GamsModifier(y), + gams.GamsModifier(ph_W), + gams.GamsModifier(xbar), + gams.GamsModifier(rho), + gams.GamsModifier(W_on), + gams.GamsModifier(prox_on), + ], + ) + + # initialize W, rho, xbar, W_on, prox_on + crops = [ "wheat", "corn", "sugarbeets" ] + for c in crops: + ph_W.add_record(c).value = 0 + xbar.add_record(c).value = 0 + rho.add_record(c).value = 0 + W_on.add_record().value = 0 + prox_on.add_record().value = 0 + + # scenario specific data applied + assert scennum < 3, "three scenarios hardwired for now" + if scennum == 0: # below + y.add_record("wheat").value = 2.0 + y.add_record("corn").value = 2.4 + y.add_record("sugarbeets").value = 16.0 + elif scennum == 1: # average + y.add_record("wheat").value = 2.5 + y.add_record("corn").value = 3.0 + y.add_record("sugarbeets").value = 20.0 + elif scennum == 2: # above + y.add_record("wheat").value = 3.0 + y.add_record("corn").value = 3.6 + y.add_record("sugarbeets").value = 24.0 + + mi.solve() + areaVarDatas = list(mi.sync_db["x"]) + + print(f"after simple solve for scenario {scennum}") + for i,v in enumerate(areaVarDatas): + print(f"{i =}, {v.get_level() =}") + + # In general, be sure to process variables in the same order has the guest does (so indexes match) + nonant_names_dict = {("ROOT",i): ("x", v.key(0)) for i, v in enumerate(areaVarDatas)} + gd = { + "scenario": mi, + "nonants": {("ROOT",i): v for i,v in enumerate(areaVarDatas)}, + "nonant_fixedness": {("ROOT",i): v.get_lower() == v.get_upper() for i,v in enumerate(areaVarDatas)}, + "nonant_start": {("ROOT",i): v.get_level() for i,v in enumerate(areaVarDatas)}, + "nonant_names": nonant_names_dict, + "nameset": {nt[0] for nt in nonant_names_dict.values()}, + "probability": "uniform", + "sense": None, + "BFs": None, + "ph" : { + "ph_W" : {("ROOT",i): p for i,p in enumerate(ph_W)}, + "xbar" : {("ROOT",i): p for i,p in enumerate(xbar)}, + "rho" : {("ROOT",i): p for i,p in enumerate(rho)}, + "W_on" : W_on.first_record(), + "prox_on" : prox_on.first_record(), + "obj" : mi.sync_db["negprofit"].find_record(), + "nonant_lbs" : {("ROOT",i): v.get_lower() for i,v in enumerate(areaVarDatas)}, + "nonant_ubs" : {("ROOT",i): v.get_upper() for i,v in enumerate(areaVarDatas)}, + }, + } + + return gd + + +if __name__ == "__main__": + + scennum = 2 + gd = create_model(scennum=scennum) + mi = gd["scenario"] + mi.solve() + print(f"iter 0 {mi.model_status =}") + print(f"after solve in main for scenario {scennum}") + for ndn_i,gxvar in gd["nonants"].items(): + print(f"{ndn_i =}, {gxvar.get_level() =}") + print("That was bad, but if we do sync_db in the right way, it will be cool") + ###mi.sync_db["x"] # not good enough, we need to make the list for some reason + # the set of names will be just x for farmer + print(f"{gd['nameset'] =}") + for n in gd["nameset"]: + list(mi.sync_db[n]) + # I don't really understand this sync_db thing; the iterator seems to have a side-effect + # But maybe the objects are the objects, so the sync has been done... yes! + for ndn_i,gxvar in gd["nonants"].items(): + print(f"{ndn_i =}, {gxvar.get_level() =}") + print("was that good?") + print("now solve again, get and display levels") + mi.solve() + print(f" iter 0 repeat solve {mi.model_status =}") + for n in gd["nameset"]: + list(mi.sync_db[n]) + for ndn_i,gxvar in gd["nonants"].items(): + print(f"{ndn_i =}, {gxvar.get_level() =}") + print(f" after repeat iter 0 solve {mi.model_status =}") + + print("Here is where the trouble starts") + print("\n Now let's try to simulate an iter 1 solve") + print(f'{gd["ph"]["prox_on"] =}') + print(f'{mi.sync_db["prox_on"].find_record() =}') + #gd["ph"]["prox_on"].set_value(1) + gd["ph"]["prox_on"].value = 1 + for ndn_i in gd["nonants"]: + print(f'{gd["ph"]["rho"][ndn_i] =}') + #gd["ph"]["rho"][ndn_i].set_value(1) + #gd["ph"]["rho"][ndn_i].value = 1 + ###gd["ph"]["xbar"][ndn_i].set_value(100) + mi.solve() + print(f" regular iter {mi.model_status =}") + print("Note that the levels do not update with status of 19") + """ + for n in gd["nameset"]: + list(mi.sync_db[n]) + for ndn_i,gxvar in gd["nonants"].items(): + print(f"{ndn_i =}, {gxvar.get_level() =}") + """ + + diff --git a/examples/farmer/GAMS/installation.txt b/examples/farmer/GAMS/installation.txt new file mode 100644 index 00000000..7bbbb350 --- /dev/null +++ b/examples/farmer/GAMS/installation.txt @@ -0,0 +1 @@ +pip install gamspy_base gamsapi diff --git a/examples/farmer/agnostic/README.txt b/examples/farmer/agnostic/README.txt new file mode 100644 index 00000000..fce34a1e --- /dev/null +++ b/examples/farmer/agnostic/README.txt @@ -0,0 +1,2 @@ +Everything in the directory is unique to farmer. See mpisppy.agnostic +and it examples directory for more general code. diff --git a/examples/farmer/agnostic/ag_gurobipy.bash b/examples/farmer/agnostic/ag_gurobipy.bash new file mode 100755 index 00000000..5f5354fd --- /dev/null +++ b/examples/farmer/agnostic/ag_gurobipy.bash @@ -0,0 +1,9 @@ +#!/bin/bash + +SOLVERNAME=gurobi + +# mpiexec -np 1 python -m mpi4py agnostic_gurobipy_cylinders.py --num-scens 3 --default-rho 1 --solver-name $SOLVERNAME --max-iterations=10 --rel-gap 0.01 --display-progress + +# mpiexec -np 2 python -m mpi4py agnostic_gurobipy_cylinders.py --num-scens 3 --default-rho 1 --solver-name $SOLVERNAME --max-iterations=10 --xhatshuffle --rel-gap 0.01 --display-progress + +mpiexec -np 3 python -m mpi4py agnostic_gurobipy_cylinders.py --num-scens 3 --default-rho 1 --solver-name $SOLVERNAME --max-iterations=10 --xhatshuffle --lagrangian --rel-gap 0.01 --display-progress diff --git a/examples/farmer/agnostic/ag_pyomo.bash b/examples/farmer/agnostic/ag_pyomo.bash new file mode 100644 index 00000000..8e30597c --- /dev/null +++ b/examples/farmer/agnostic/ag_pyomo.bash @@ -0,0 +1,15 @@ +#!/bin/bash + + + +#python agnostic_cylinders.py --help + +#mpiexec -np 3 python -m mpi4py farmer_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=5 --xhatshuffle --lagrangian --rel-gap 0.01 + +#mpiexec -np 3 python -m mpi4py agnostic_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=5 --xhatshuffle --lagrangian --rel-gap 0.01 + +#python agnostic_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=2 + +#mpiexec -np 2 python -m mpi4py agnostic_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=10 --lagrangian --rel-gap 0.01 + +mpiexec -np 2 python -m mpi4py agnostic_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=10 --xhatshuffle --rel-gap 0.01 diff --git a/examples/farmer/agnostic/ag_pyomo_ph.bash b/examples/farmer/agnostic/ag_pyomo_ph.bash new file mode 100644 index 00000000..7b4e0332 --- /dev/null +++ b/examples/farmer/agnostic/ag_pyomo_ph.bash @@ -0,0 +1,15 @@ +#!/bin/bash + + + +#python agnostic_cylinders.py --help + +#mpiexec -np 3 python -m mpi4py farmer_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=5 --xhatshuffle --lagrangian --rel-gap 0.01 + +#mpiexec -np 3 python -m mpi4py agnostic_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=5 --xhatshuffle --lagrangian --rel-gap 0.01 + +#python agnostic_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=2 + +#mpiexec -np 2 python -m mpi4py agnostic_pyomo_cylinders.py --num-scens 3 --default-rho 1 --solver-name cplex --max-iterations=10 --lagrangian --rel-gap 0.01 + +python agnostic_pyomo_ph.py --num-scens 3 --default-rho 1 --solver-name gurobi --max-iterations=10 diff --git a/examples/farmer/agnostic/agnostic_ampl_cylinders.py b/examples/farmer/agnostic/agnostic_ampl_cylinders.py new file mode 100644 index 00000000..573d86b4 --- /dev/null +++ b/examples/farmer/agnostic/agnostic_ampl_cylinders.py @@ -0,0 +1,83 @@ +############################################################################### +# 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. +############################################################################### +# Started by dlw Aug 2023 + +import farmer_ampl_agnostic +from mpisppy.spin_the_wheel import WheelSpinner +import mpisppy.utils.cfg_vanilla as vanilla +import mpisppy.utils.config as config +import mpisppy.agnostic.agnostic as agnostic +import mpisppy.utils.sputils as sputils + +def _farmer_parse_args(): + # create a config object and parse JUST FOR TESTING + cfg = config.Config() + + farmer_ampl_agnostic.inparser_adder(cfg) + + 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.parse_command_line("farmer_ampl_agnostic_cylinders") + return cfg + + +if __name__ == "__main__": + print("begin ad hoc main for agnostic.py") + + cfg = _farmer_parse_args() + Ag = agnostic.Agnostic(farmer_ampl_agnostic, cfg) + + scenario_creator = Ag.scenario_creator + scenario_denouement = farmer_ampl_agnostic.scenario_denouement # should we go though Ag? + all_scenario_names = ['scen{}'.format(sn) for sn in range(cfg.num_scens)] + + # Things needed for vanilla cylinders + beans = (cfg, scenario_creator, scenario_denouement, all_scenario_names) + + # Vanilla PH hub + hub_dict = vanilla.ph_hub(*beans, + scenario_creator_kwargs=None, # kwargs in Ag not here + ph_extensions=None, + ph_converger=None, + rho_setter = None) + # pass the Ag object via options... + hub_dict["opt_kwargs"]["options"]["Ag"] = Ag + + # xhat shuffle bound spoke + if cfg.xhatshuffle: + xhatshuffle_spoke = vanilla.xhatshuffle_spoke(*beans, scenario_creator_kwargs=None) + xhatshuffle_spoke["opt_kwargs"]["options"]["Ag"] = Ag + if cfg.lagrangian: + lagrangian_spoke = vanilla.lagrangian_spoke(*beans, scenario_creator_kwargs=None) + lagrangian_spoke["opt_kwargs"]["options"]["Ag"] = Ag + + list_of_spoke_dict = list() + if cfg.xhatshuffle: + list_of_spoke_dict.append(xhatshuffle_spoke) + if cfg.lagrangian: + list_of_spoke_dict.append(lagrangian_spoke) + + wheel = WheelSpinner(hub_dict, list_of_spoke_dict) + wheel.spin() + + write_solution = False + 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') + diff --git a/examples/farmer/agnostic/agnostic_gurobipy_cylinders.py b/examples/farmer/agnostic/agnostic_gurobipy_cylinders.py new file mode 100644 index 00000000..fb4cb74f --- /dev/null +++ b/examples/farmer/agnostic/agnostic_gurobipy_cylinders.py @@ -0,0 +1,83 @@ +############################################################################### +# 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. +############################################################################### +# This software is distributed under the 3-clause BSD License. + +import farmer_gurobipy_agnostic +from mpisppy.spin_the_wheel import WheelSpinner +import mpisppy.utils.cfg_vanilla as vanilla +import mpisppy.utils.config as config +import mpisppy.agnostic.agnostic as agnostic +import mpisppy.utils.sputils as sputils + +def _farmer_parse_args(): + # create a config object and parse JUST FOR TESTING + cfg = config.Config() + + farmer_gurobipy_agnostic.inparser_adder(cfg) + + 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.parse_command_line("farmer_gurobipy_agnostic_cylinders") + return cfg + + +if __name__ == "__main__": + print("begin ad hoc main for agnostic.py") + + cfg = _farmer_parse_args() + Ag = agnostic.Agnostic(farmer_gurobipy_agnostic, cfg) + + scenario_creator = Ag.scenario_creator + scenario_denouement = farmer_gurobipy_agnostic.scenario_denouement # should we go though Ag? + all_scenario_names = ['scen{}'.format(sn) for sn in range(cfg.num_scens)] + + # Things needed for vanilla cylinders + beans = (cfg, scenario_creator, scenario_denouement, all_scenario_names) + + # Vanilla PH hub + hub_dict = vanilla.ph_hub(*beans, + scenario_creator_kwargs=None, # kwargs in Ag not here + ph_extensions=None, + ph_converger=None, + rho_setter = None) + # pass the Ag object via options... + hub_dict["opt_kwargs"]["options"]["Ag"] = Ag + + # xhat shuffle bound spoke + if cfg.xhatshuffle: + xhatshuffle_spoke = vanilla.xhatshuffle_spoke(*beans, scenario_creator_kwargs=None) + xhatshuffle_spoke["opt_kwargs"]["options"]["Ag"] = Ag + if cfg.lagrangian: + lagrangian_spoke = vanilla.lagrangian_spoke(*beans, scenario_creator_kwargs=None) + lagrangian_spoke["opt_kwargs"]["options"]["Ag"] = Ag + + list_of_spoke_dict = list() + if cfg.xhatshuffle: + list_of_spoke_dict.append(xhatshuffle_spoke) + if cfg.lagrangian: + list_of_spoke_dict.append(lagrangian_spoke) + + wheel = WheelSpinner(hub_dict, list_of_spoke_dict) + wheel.spin() + + write_solution = False + 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') + diff --git a/examples/farmer/agnostic/agnostic_pyomo_cylinders.py b/examples/farmer/agnostic/agnostic_pyomo_cylinders.py new file mode 100644 index 00000000..a3614a04 --- /dev/null +++ b/examples/farmer/agnostic/agnostic_pyomo_cylinders.py @@ -0,0 +1,84 @@ +############################################################################### +# 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. +############################################################################### +# This software is distributed under the 3-clause BSD License. +# Started by dlw Aug 2023 + +import farmer_pyomo_agnostic +from mpisppy.spin_the_wheel import WheelSpinner +import mpisppy.utils.cfg_vanilla as vanilla +import mpisppy.utils.config as config +import mpisppy.agnostic.agnostic as agnostic +import mpisppy.utils.sputils as sputils + +def _farmer_parse_args(): + # create a config object and parse JUST FOR TESTING + cfg = config.Config() + + farmer_pyomo_agnostic.inparser_adder(cfg) + + 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.parse_command_line("farmer_pyomo_agnostic_cylinders") + return cfg + + +if __name__ == "__main__": + print("begin ad hoc main for agnostic.py") + + cfg = _farmer_parse_args() + Ag = agnostic.Agnostic(farmer_pyomo_agnostic, cfg) + + scenario_creator = Ag.scenario_creator + scenario_denouement = farmer_pyomo_agnostic.scenario_denouement # should we go though Ag? + all_scenario_names = ['scen{}'.format(sn) for sn in range(cfg.num_scens)] + + # Things needed for vanilla cylinders + beans = (cfg, scenario_creator, scenario_denouement, all_scenario_names) + + # Vanilla PH hub + hub_dict = vanilla.ph_hub(*beans, + scenario_creator_kwargs=None, # kwargs in Ag not here + ph_extensions=None, + ph_converger=None, + rho_setter = None) + # pass the Ag object via options... + hub_dict["opt_kwargs"]["options"]["Ag"] = Ag + + # xhat shuffle bound spoke + if cfg.xhatshuffle: + xhatshuffle_spoke = vanilla.xhatshuffle_spoke(*beans, scenario_creator_kwargs=None) + xhatshuffle_spoke["opt_kwargs"]["options"]["Ag"] = Ag + if cfg.lagrangian: + lagrangian_spoke = vanilla.lagrangian_spoke(*beans, scenario_creator_kwargs=None) + lagrangian_spoke["opt_kwargs"]["options"]["Ag"] = Ag + + list_of_spoke_dict = list() + if cfg.xhatshuffle: + list_of_spoke_dict.append(xhatshuffle_spoke) + if cfg.lagrangian: + list_of_spoke_dict.append(lagrangian_spoke) + + wheel = WheelSpinner(hub_dict, list_of_spoke_dict) + wheel.spin() + + write_solution = False + 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') + diff --git a/examples/farmer/agnostic/agnostic_pyomo_ph.py b/examples/farmer/agnostic/agnostic_pyomo_ph.py new file mode 100644 index 00000000..34cc1ba8 --- /dev/null +++ b/examples/farmer/agnostic/agnostic_pyomo_ph.py @@ -0,0 +1,69 @@ +############################################################################### +# 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. +############################################################################### +# This software is distributed under the 3-clause BSD License. +# Started by dlw Aug 2023 + +import farmer_pyomo_agnostic +import mpisppy.utils.cfg_vanilla as vanilla +import mpisppy.utils.config as config +import mpisppy.agnostic.agnostic as agnostic + +def _farmer_parse_args(): + # create a config object and parse JUST FOR TESTING + cfg = config.Config() + + farmer_pyomo_agnostic.inparser_adder(cfg) + + 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.parse_command_line("farmer_pyomo_agnostic_cylinders") + return cfg + + +if __name__ == "__main__": + print("begin ad hoc main for agnostic.py") + + cfg = _farmer_parse_args() + Ag = agnostic.Agnostic(farmer_pyomo_agnostic, cfg) + + scenario_creator = Ag.scenario_creator + scenario_denouement = farmer_pyomo_agnostic.scenario_denouement # should we go though Ag? + all_scenario_names = ['scen{}'.format(sn) for sn in range(cfg.num_scens)] + + # Things needed for vanilla cylinders + beans = (cfg, scenario_creator, scenario_denouement, all_scenario_names) + + # Vanilla PH hub + hub_dict = vanilla.ph_hub(*beans, + scenario_creator_kwargs=None, # kwargs in Ag not here + ph_extensions=None, + ph_converger=None, + rho_setter = None) + # pass the Ag object via options... + hub_dict["opt_kwargs"]["options"]["Ag"] = Ag + + ph = hub_dict["opt_class"](**hub_dict["opt_kwargs"]) + + conv, obj, triv = ph.ph_main() + # Obj includes prox (only ok if we find a non-ant soln) + if (conv < 1e-8): + print(f'Objective value: {obj:.2f}') + else: + print('Did not find a non-anticipative solution ' + f'(conv = {conv:.1e})') + + ph.post_solve_bound(verbose=False) diff --git a/examples/farmer/agnostic/farmer.mod b/examples/farmer/agnostic/farmer.mod new file mode 100644 index 00000000..5a5f147e --- /dev/null +++ b/examples/farmer/agnostic/farmer.mod @@ -0,0 +1,111 @@ +# The farmer's problem in AMPL +# +# Reference: +# John R. Birge and Francois Louveaux. Introduction to Stochastic Programming. +# +# AMPL coding by Victor Zverovich; ## modifed by dlw; now *minization* + +##function expectation; +##function random; + +##suffix stage IN; + +set Crops; + +##set Scen; +##param P{Scen}; # probabilities + +param TotalArea; # acre +param PlantingCost{Crops}; # $/acre +param SellingPrice{Crops}; # $/T +param ExcessSellingPrice; # $/T +param PurchasePrice{Crops}; # $/T +param MinRequirement{Crops}; # T +param BeetsQuota; # T + +# Area in acres devoted to crop c. +var area{c in Crops} >= 0; + +# Tons of crop c sold (at favourable price) under scenario s. +var sell{c in Crops} >= 0, suffix stage 2; + +# Tons of sugar beets sold in excess of the quota under scenario s. +var sell_excess >= 0, suffix stage 2; + +# Tons of crop c bought under scenario s +var buy{c in Crops} >= 0, suffix stage 2; + +# The random variable (parameter) representing the yield of crop c. +##var RandomYield{c in Crops}; +param RandomYield{c in Crops}; + +# Realizations of the yield of crop c. +##param Yield{c in Crops, s in Scen}; # T/acre + +##maximize profit: +## expectation( +## ExcessSellingPrice * sell_excess + +## sum{c in Crops} (SellingPrice[c] * sell[c] - +## PurchasePrice[c] * buy[c])) - +## sum{c in Crops} PlantingCost[c] * area[c]; + +minimize minus_profit: + - ExcessSellingPrice * sell_excess - + sum{c in Crops} (SellingPrice[c] * sell[c] - + PurchasePrice[c] * buy[c]) + + sum{c in Crops} (PlantingCost[c] * area[c]); + +s.t. totalArea: sum {c in Crops} area[c] <= TotalArea; + +s.t. requirement{c in Crops}: + RandomYield[c] * area[c] - sell[c] + buy[c] >= MinRequirement[c]; + +s.t. quota: sell['beets'] <= BeetsQuota; + +s.t. sellBeets: + sell['beets'] + sell_excess <= RandomYield['beets'] * area['beets']; + +##yield: random({c in Crops} (RandomYield[c], {s in Scen} Yield[c, s])); + +data; + +set Crops := wheat corn beets; +#set Scen := below average above; + +param TotalArea := 500; + +##param Yield: +## below average above := +## wheat 2.0 2.5 3.0 +## corn 2.4 3.0 3.6 +## beets 16.0 20.0 24.0; + +# Average Scenario +param RandomYield := + wheat 2.5 + corn 3.0 + beets 20.0; + +param PlantingCost := + wheat 150 + corn 230 + beets 260; + +param SellingPrice := + wheat 170 + corn 150 + beets 36; + +param ExcessSellingPrice := 10; + +param PurchasePrice := + wheat 238 + corn 210 + beets 100; + +param MinRequirement := + wheat 200 + corn 240 + beets 0; + +param BeetsQuota := 6000; diff --git a/examples/farmer/agnostic/farmer.py b/examples/farmer/agnostic/farmer.py new file mode 100644 index 00000000..8f40726d --- /dev/null +++ b/examples/farmer/agnostic/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/agnostic/farmer_ampl_agnostic.py b/examples/farmer/agnostic/farmer_ampl_agnostic.py new file mode 100644 index 00000000..a4022b73 --- /dev/null +++ b/examples/farmer/agnostic/farmer_ampl_agnostic.py @@ -0,0 +1,400 @@ +############################################################################### +# 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. +############################################################################### +# +# In this example, AMPL is the guest language. +# ***This is a special example where this file serves +# the rolls of ampl_guest.py and the python model file.*** + +""" +This file tries to show many ways to do things in AMPLpy, +but not necessarily the best ways in all cases. +""" + +from amplpy import AMPL +import pyomo.environ as pyo +import mpisppy.utils.sputils as sputils +import numpy as np + +from mpisppy import MPI # for debugging +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +# If you need random numbers, 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 + + NOTE: for ampl, the names will be tuples name, index + """ + + assert crops_multiplier == 1, "for AMPL, just getting started with 3 crops" + + ampl = AMPL() + + ampl.read("farmer.mod") + + # scenario specific data applied + scennum = sputils.extract_num(scenario_name) + assert scennum < 3, "three scenarios hardwired for now" + y = ampl.get_parameter("RandomYield") + if scennum == 0: # below + y.set_values({"wheat": 2.0, "corn": 2.4, "beets": 16.0}) + elif scennum == 2: # above + y.set_values({"wheat": 3.0, "corn": 3.6, "beets": 24.0}) + + areaVarDatas = list(ampl.get_variable("area").instances()) + + # In general, be sure to process variables in the same order has the guest does (so indexes match) + gd = { + "scenario": ampl, + "nonants": {("ROOT",i): v[1] for i,v in enumerate(areaVarDatas)}, + "nonant_fixedness": {("ROOT",i): v[1].astatus()=="fixed" for i,v in enumerate(areaVarDatas)}, + "nonant_start": {("ROOT",i): v[1].value() for i,v in enumerate(areaVarDatas)}, + "nonant_names": {("ROOT",i): ("area", v[0]) for i, v in enumerate(areaVarDatas)}, + "probability": "uniform", + "sense": pyo.minimize, + "BFs": None + } + + return gd + +#========= +def scenario_names_creator(num_scens,start=None): + if (start is None) : + start=0 + return [f"scen{i}" for i in range(start,start+num_scens)] + + +#========= +def inparser_adder(cfg): + 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): + # creates keywords for scenario creator + 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 + + +# This is not needed for PH +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + # 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) + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + pass + # (the fct in farmer won't work because the Var names don't match) + #farmer.scenario_denouement(rank, scenario_name, scenario) + + +################################################################################################## +# begin callouts +# NOTE: the callouts all take the Ag object as their first argument, mainly to see cfg if needed +# the function names correspond to function names in mpisppy + +def attach_Ws_and_prox(Ag, sname, scenario): + # this is AMPL farmer specific, so we know there is not a W already, e.g. + # Attach W's and rho to the guest scenario (mutable params). + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + # (there must be some way to create and assign *mutable* params in on call to AMPL) + gs.eval("param W_on;") + gs.eval("let W_on := 0;") + gs.eval("param prox_on;") + gs.eval("let prox_on := 0;") + # we are trusting the order to match the nonant indexes + gs.eval("param W{Crops};") + # Note: we should probably use set_values instead of let + gs.eval("let {c in Crops} W[c] := 0;") + # start with rho at zero, but update before solve + gs.eval("param rho{Crops};") + gs.eval("let {c in Crops} rho[c] := 0;") + + +def _disable_prox(Ag, scenario): + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gs.get_parameter("prox_on").set(0) + + +def _disable_W(Ag, scenario): + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gs.get_parameter("W_on").set(0) + + +def _reenable_prox(Ag, scenario): + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gs.get_parameter("prox_on").set(1) + + +def _reenable_W(Ag, scenario): + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gs.get_parameter("W_on").set(1) + + +def attach_PH_to_objective(Ag, sname, scenario, add_duals, add_prox): + # Deal with prox linearization and approximation later, + # i.e., just do the quadratic version + + # The host has xbars and computes without involving the guest language + gd = scenario._agnostic_dict + gs = gd["scenario"] # guest scenario handle + gs.eval("param xbars{Crops};") + + # Dual term (weights W) + try: + profitobj = gs.get_objective("minus_profit") + except: + print("big troubles!!; we can't find the objective function") + print("doing export to export.mod") + gs.export_model("export.mod") + raise + + objstr = str(profitobj) + phobjstr = "" + if add_duals: + phobjstr += " + W_on * sum{c in Crops} (W[c] * area[c])" + + # Prox term (quadratic) + if add_prox: + """ + prox_expr = 0. + for ndn_i, xvar in gd["nonants"].items(): + # expand (x - xbar)**2 to (x**2 - 2*xbar*x + xbar**2) + # x**2 is the only qradratic term, which might be + # dealt with differently depending on user-set options + if xvar.is_binary(): + xvarsqrd = xvar + else: + xvarsqrd = xvar**2 + prox_expr += (gs.rho[ndn_i] / 2.0) * \ + (xvarsqrd - 2.0 * xbars[ndn_i] * xvar + xbars[ndn_i]**2) + """ + phobjstr += " + prox_on * sum{c in Crops} ((rho[c]/2.0) * (area[c] * area[c] "+\ + " - 2.0 * xbars[c] * area[c] + xbars[c]^2))" + objstr = objstr[:-1] + "+ (" + phobjstr + ");" + objstr = objstr.replace("minimize minus_profit", "minimize phobj") + profitobj.drop() + gs.eval(objstr) + gs.eval("delete minus_profit;") + currentobj = gs.get_current_objective() + # see _copy_Ws_... see also the gams version + WParamDatas = list(gs.get_parameter("W").instances()) + xbarsParamDatas = list(gs.get_parameter("xbars").instances()) + rhoParamDatas = list(gs.get_parameter("rho").instances()) + gd["PH"] = { + "W": {("ROOT",i): v for i,v in enumerate(WParamDatas)}, + "xbars": {("ROOT",i): v for i,v in enumerate(xbarsParamDatas)}, + "rho": {("ROOT",i): v for i,v in enumerate(rhoParamDatas)}, + "obj": currentobj, + } + + +def solve_one(Ag, s, solve_keyword_args, gripe, tee, need_solution=True): + # This needs to attach stuff to s (see solve_one in spopt.py) + # Solve the guest language version, then copy values to the host scenario + + # This function needs to W on the guest right before the solve + + # We need to operate on the guest scenario, not s; however, attach things to s (the host scenario) + # and copy to s. If you are working on a new guest, you should not have to edit the s side of things + + # To acommdate the solve_one call from xhat_eval.py, we need to attach the obj fct value to s + + # time.sleep(np.random.uniform()/10) + + _copy_Ws_xbars_rho_from_host(s) + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + + #### start debugging + if False: # True: # global_rank == 0: + try: + WParamDatas = list(gs.get_parameter("W").instances()) + print(f" ^^^ in _solve_one {WParamDatas =} {global_rank =}") + except: # noqa + print(f" ^^^^ no W for xhat {global_rank=}") + #prox_on = gs.get_parameter("prox_on").value() + #print(f" ^^^ in _solve_one {prox_on =} {global_rank =}") + #W_on = gs.get_parameter("W_on").value() + #print(f" ^^^ in _solve_one {W_on =} {global_rank =}") + #xbarsParamDatas = list(gs.get_parameter("xbars").instances()) + #print(f" in _solve_one {xbarsParamDatas =} {global_rank =}") + #rhoParamDatas = list(gs.get_parameter("rho").instances()) + #print(f" in _solve_one {rhoParamDatas =} {global_rank =}") + #### stop debugging + + solver_name = s._solver_plugin.name + gs.set_option("solver", solver_name) + if 'persistent' in solver_name: + raise RuntimeError("Persistent solvers are not currently supported in the farmer agnostic example.") + gs.set_option("presolve", 0) + + solver_exception = None + try: + gs.solve() + except Exception as e: + solver_exception = e + + # debug + #fname = f"{s.name}_{global_rank}" + #print(f"debug export to {fname}") + #gs.export_model(f"{fname}.mod") + #gs.export_data(f"{fname}.dat") + + if gs.solve_result != "solved": + s._mpisppy_data.scenario_feasible = False + if gripe: + print (f"Solve failed for scenario {s.name} on rank {global_rank}") + print(f"{gs.solve_result =}") + + if solver_exception is not None and need_solution: + raise solver_exception + + + s._mpisppy_data.scenario_feasible = True + # For AMPL mips, we need to use the gap option to compute bounds + # https://amplmp.readthedocs.io/rst/features-guide.html + objobj = gs.get_current_objective() # different for xhatters + objval = objobj.value() + if gd["sense"] == pyo.minimize: + s._mpisppy_data.outer_bound = objval + else: + s._mpisppy_data.inner_bound = objval + + # copy the nonant x values from gs to s so mpisppy can use them in s + # in general, we need more checks (see the pyomo agnostic guest example) + for ndn_i, gxvar in gd["nonants"].items(): + try: # not sure this is needed + float(gxvar.value()) + except: # noqa + raise RuntimeError( + f"Non-anticipative variable {gxvar.name} on scenario {s.name} " + "had no value. This usually means this variable " + "did not appear in any (active) components, and hence " + "was not communicated to the subproblem solver. ") + if gxvar.astatus() == "pre": + raise RuntimeError( + f"Non-anticipative variable {gxvar.name} on scenario {s.name} " + "was presolved out. This usually means this variable " + "did not appear in any (active) components, and hence " + "was not communicated to the subproblem solver. ") + + s._mpisppy_data.nonant_indices[ndn_i]._value = gxvar.value() + + # the next line ignores bundling + s._mpisppy_data._obj_from_agnostic = objval + + # TBD: deal with other aspects of bundling (see solve_one in spopt.py) + + +# local helper called right before the solve +def _copy_Ws_xbars_rho_from_host(s): + # print(f" debug copy_Ws {s.name =}, {global_rank =}") + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + + # We can't use a simple list because of indexes, we have to use a dict + # NOTE that we know that W is indexed by crops for this problem + # and the nonant_names are tuple with the index in the 1 slot + # AMPL params are tuples (index, value), which are immutable + if hasattr(s._mpisppy_model, "W"): + Wdict = {gd["nonant_names"][ndn_i][1]:\ + pyo.value(v) for ndn_i, v in s._mpisppy_model.W.items()} + gs.get_parameter("W").set_values(Wdict) + rhodict = {gd["nonant_names"][ndn_i][1]:\ + pyo.value(v) for ndn_i, v in s._mpisppy_model.rho.items()} + gs.get_parameter("rho").set_values(rhodict) + xbarsdict = {gd["nonant_names"][ndn_i][1]:\ + pyo.value(v) for ndn_i, v in s._mpisppy_model.xbars.items()} + gs.get_parameter("xbars").set_values(xbarsdict) + # debug + fname = f"{s.name}_{global_rank}" + print(f"debug export to {fname}") + gs.export_model(f"{fname}.mod") + gs.export_data(f"{fname}.dat") + + else: + pass # presumably an xhatter; we should check, I suppose + + +# local helper +def _copy_nonants_from_host(s): + # values and fixedness; + gd = s._agnostic_dict + for ndn_i, gxvar in gd["nonants"].items(): + hostVar = s._mpisppy_data.nonant_indices[ndn_i] + guestVar = gd["nonants"][ndn_i] + if guestVar.astatus() == "fixed": + guestVar.unfix() + if hostVar.is_fixed(): + guestVar.fix(hostVar._value) + else: + guestVar.set_value(hostVar._value) + + +def _restore_nonants(Ag, s): + # the host has already restored + _copy_nonants_from_host(s) + + +def _restore_original_fixedness(Ag, s): + # The host has restored already + # Note that this also takes values from the host, which should be OK + _copy_nonants_from_host(s) + + +def _fix_nonants(Ag, s): + # the host has already fixed + _copy_nonants_from_host(s) + + +def _fix_root_nonants(Ag, s): + # the host has already fixed + _copy_nonants_from_host(s) + + + diff --git a/examples/farmer/agnostic/farmer_gurobipy_agnostic.py b/examples/farmer/agnostic/farmer_gurobipy_agnostic.py new file mode 100644 index 00000000..c82ae013 --- /dev/null +++ b/examples/farmer/agnostic/farmer_gurobipy_agnostic.py @@ -0,0 +1,350 @@ +############################################################################### +# 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. +############################################################################### +# In this example, Gurobipy is the guest language +import gurobipy as gp +from gurobipy import GRB +import pyomo.environ as pyo +import mpisppy.utils.sputils as sputils +import numpy as np +# debugging +from mpisppy import MPI +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +farmerstream = np.random.RandomState() + +def scenario_creator(scenario_name, use_integer=False, sense=GRB.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): + gurobipy sense (minimization or maximization). Must be either + GRB.MINIMIZE or GRB.MAXIMIZE. Default is GRB.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 + + NOTE: + """ + scennum = sputils.extract_num(scenario_name) + basenames = ['BelowAverageScenario', 'AverageScenario', 'AboveAverageScenario'] + basenum = scennum % 3 + groupnum = scennum // 3 + scenname = basenames[basenum] + str(groupnum) + + farmerstream.seed(scennum + seedoffset) + + if sense not in [GRB.MINIMIZE, GRB.MAXIMIZE]: + raise ValueError("Model sense Not recognized") + + model = gp.Model(scenname) + # Silence gurobi output + model.setParam('OutputFlag', 0) + + crops = ["WHEAT", "CORN", "SUGAR_BEETS"] + CROPS = [f"{crop}{i}" for i in range(crops_multiplier) for crop in crops] + + # Data + TOTAL_ACREAGE = 500.0 * crops_multiplier + + def get_scaled_data(indict): + outdict = {} + for i in range(crops_multiplier): + for crop in crops: + outdict[f"{crop}{i}"] = indict[crop] + return outdict + + PriceQuota = get_scaled_data({'WHEAT': 100000.0, 'CORN': 100000.0, 'SUGAR_BEETS': 6000.0}) + SubQuotaSellingPrice = get_scaled_data({'WHEAT': 170.0, 'CORN': 150.0, 'SUGAR_BEETS': 36.0}) + SuperQuotaSellingPrice = get_scaled_data({'WHEAT': 0.0, 'CORN': 0.0, 'SUGAR_BEETS': 10.0}) + CattleFeedRequirement = get_scaled_data({'WHEAT': 200.0, 'CORN': 240.0, 'SUGAR_BEETS': 0.0}) + PurchasePrice = get_scaled_data({'WHEAT': 238.0, 'CORN': 210.0, 'SUGAR_BEETS': 100000.0}) + PlantingCostPerAcre = get_scaled_data({'WHEAT': 150.0, 'CORN': 230.0, 'SUGAR_BEETS': 260.0}) + + Yield = { + 'BelowAverageScenario': {'WHEAT': 2.0, 'CORN': 2.4, 'SUGAR_BEETS': 16.0}, + 'AverageScenario': {'WHEAT': 2.5, 'CORN': 3.0, 'SUGAR_BEETS': 20.0}, + 'AboveAverageScenario': {'WHEAT': 3.0, 'CORN': 3.6, 'SUGAR_BEETS': 24.0} + } + + yield_vals = {crop: Yield[basenames[basenum]][crop.rstrip("0123456789")] + (farmerstream.rand() if groupnum != 0 else 0) for crop in CROPS} + + # Variables + DevotedAcreage = model.addVars(CROPS, vtype=GRB.INTEGER if use_integer else GRB.CONTINUOUS, lb=0.0, ub=TOTAL_ACREAGE, name="DevotedAcreage") + QuantitySubQuotaSold = model.addVars(CROPS, lb=0.0, name="QuantitySubQuotaSold") + QuantitySuperQuotaSold = model.addVars(CROPS, lb=0.0, name="QuantitySuperQuotaSold") + QuantityPurchased = model.addVars(CROPS, lb=0.0, name="QuantityPurchased") + + # Constraints + model.addConstr(gp.quicksum(DevotedAcreage[crop] for crop in CROPS) <= TOTAL_ACREAGE, "TotalAcreage") + + for crop in CROPS: + model.addConstr(CattleFeedRequirement[crop] <= yield_vals[crop] * DevotedAcreage[crop] + QuantityPurchased[crop] - QuantitySubQuotaSold[crop] - QuantitySuperQuotaSold[crop], f"CattleFeedReq_{crop}") + model.addConstr(QuantitySubQuotaSold[crop] + QuantitySuperQuotaSold[crop] - (yield_vals[crop] * DevotedAcreage[crop]) <= 0.0, f"LimitAmountSold_{crop}") + model.addConstr(QuantitySubQuotaSold[crop] <= PriceQuota[crop], f"EnforceQuota_{crop}") + + # Objective + total_costs = gp.quicksum(PlantingCostPerAcre[crop] * DevotedAcreage[crop] for crop in CROPS) + purchase_costs = gp.quicksum(PurchasePrice[crop] * QuantityPurchased[crop] for crop in CROPS) + subquota_revenue = gp.quicksum(SubQuotaSellingPrice[crop] * QuantitySubQuotaSold[crop] for crop in CROPS) + superquota_revenue = gp.quicksum(SuperQuotaSellingPrice[crop] * QuantitySuperQuotaSold[crop] for crop in CROPS) + + total_cost = total_costs + purchase_costs - subquota_revenue - superquota_revenue + model.setObjective(total_cost, sense) + + model.optimize() + + gd = { + "scenario": model, + "nonants": {("ROOT", i): v for i, v in enumerate(DevotedAcreage.values())}, + "nonants_coeffs": {("ROOT", i): v.Obj for i, v in enumerate(DevotedAcreage.values())}, + "nonant_fixedness": {("ROOT", i): v.LB == v.UB for i, v in enumerate(DevotedAcreage.values())}, + "nonant_start": {("ROOT", i): v.Start for i, v in enumerate(DevotedAcreage.values())}, + "nonant_names": {("ROOT", i): v.VarName for i, v in enumerate(DevotedAcreage.values())}, + "probability": "uniform", + "sense": sense, + "BFs": None, + "nonant_bounds": {("ROOT", i): (v.LB, v.UB) for i, v in enumerate(DevotedAcreage.values())} + } + + return gd + +#========== +def scenario_names_creator(num_scens,start=None): + 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): + # creates keywords for scenario creator + 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 + +# This is not needed for PH +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + sca = scenario_creator_kwargs.copy() + sca["seedoffset"] = seed + sca["num_scens"] = sample_branching_factors[0] # two-stage problem + return scenario_creator(sname, **sca) + +def scenario_denouement(rank, scenario_name, scenario): + pass + +################################################################################################## +# begin callouts + +def attach_Ws_and_prox(Ag, sname, scenario): + """Gurobipy does not have symbolic data, so this function just collects some data""" + # Gurobipy is special so we are gonna need to maintain the coeffs ourselves + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gd = scenario._agnostic_dict + obj = gs.getObjective() + + obj_func_terms = [obj.getVar(i) for i in range(obj.size())] + nonant_coeffs = {} + nonants_not_in_obj = {} + + # Check to see if the nonants are in the objective function + for ndn_i, nonant in gd["nonants"].items(): + found_nonant_in_obj = False + for obj_func_term in obj_func_terms: + if obj_func_term.sameAs(nonant): + nonant_coeffs[nonant] = nonant.Obj + found_nonant_in_obj = True + break + if not found_nonant_in_obj: + print(f'No objective coeff for {gd["nonant_names"][ndn_i]=} which is bad so we will default to 0') + nonant_coeffs[nonant] = 0 + nonants_not_in_obj[ndn_i] = nonant + + # Update/attach nonants' coeffs to the dictionary + gd["nonant_coeffs"] = nonant_coeffs + gd["nonants_not_in_obj"] = nonants_not_in_obj + + +def _disable_prox(Ag, scenario): + pass + # raise RuntimeError("Did not expect _disable_prox") + +def _disable_W(Ag, scenario): + pass + # raise RuntimeError("Did not expect _disable_W") + +def _reenable_prox(Ag, scenario): + pass + # raise RuntimeError("Did not expect _reenable_prox") + +def _reenable_W(Ag, scenario): + pass + # raise RuntimeError("Did not expect _reenable_W") + +def attach_PH_to_objective(Ag, sname, scenario, add_duals, add_prox): + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gd = scenario._agnostic_dict + nonant_coeffs = gd["nonant_coeffs"] + + ''' + At this point we can assume that all that all the nonants are already in the objective function from attach_Ws_and_prox + but for the nonants that are not in the objective function we can just attach it to the objective funciton with thier respective coefficients + ''' + + # Adding all the nonants into the obj function + obj = gs.getObjective() + for nonant_not_in_obj in gd["nonants_not_in_obj"]: + obj += (nonant_coeffs[nonant_not_in_obj] * nonant_not_in_obj) + + # At this point all nonants are in the obj + nonant_sqs = {} + # Need to create a new var for x^2, and then add a constraint to the var with x val, so that we can set coeff value later on + for i, nonant in gd["nonants"].items(): + # Create a constaint that sets x * x = xsq + nonant_sq = gs.addVar(vtype=GRB.CONTINUOUS, obj=nonant.Obj**2, name=f"{nonant.VarName}sq") + gs.addConstr(nonant * nonant == nonant_sq, f'{nonant.VarName}sqconstr') + # Put the x^2 in the objective function + obj += nonant_sq + nonant_sqs[i] = nonant_sq + + # Update model and gd + gs.update() + gd["nonant_sqs"] = nonant_sqs + + _copy_Ws_xbars_rho_from_host(scenario) + +def solve_one(Ag, s, solve_keyword_args, gripe, tee): + _copy_Ws_xbars_rho_from_host(s) + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + + # Assuming gs is a Gurobi model, we can start solving + try: + gs.optimize() + except gp.GurobiError as e: + print(f"Error occurred: {str(e)}") + s._mpisppy_data.scenario_feasible = False + if gripe: + print(f"Solve failed for scenario {s.name}") + return + + if gs.status != gp.GRB.Status.OPTIMAL: + s._mpisppy_data.scenario_feasible = False + if gripe: + print(f"Solve failed for scenario {s.name}") + return + + s._mpisppy_data.scenario_feasible = True + + # Objective value extraction + objval = gs.getObjective().getValue() + + if gd["sense"] == gp.GRB.MINIMIZE: + s._mpisppy_data.outer_bound = objval + else: + s._mpisppy_data.outer_bound = objval + + # Copy the non-anticipative variable values from guest to host scenario + for ndn_i, gxvar in gd["nonants"].items(): + grb_var = gs.getVarByName(gxvar.VarName) + if grb_var is None: + raise RuntimeError( + f"Non-anticipative variable {gxvar.varname} on scenario {s.name} " + "was not found in the Gurobi model." + ) + s._mpisppy_data.nonant_indices[ndn_i]._value = grb_var.X + + # Store the objective function value in the host scenario + s._mpisppy_data._obj_from_agnostic = objval + + # Additional checks and operations for bundling if needed (depending on the problem) + # ... + +def _copy_Ws_xbars_rho_from_host(scenario): + # Calculates the coefficients of the new expanded objective function + # Regardless need to calculate coefficients for x^2 and x + gd = scenario._agnostic_dict + gs = scenario._agnostic_dict["scenario"] # guest handle + + # Decide if we are using PH or xhatter + if hasattr(scenario._mpisppy_model, "W"): + # Get our Ws, rhos, and xbars + Wdict = {gd["nonant_names"][ndn_i]:\ + pyo.value(v) for ndn_i, v in scenario._mpisppy_model.W.items()} + rhodict = {gd["nonant_names"][ndn_i]:\ + pyo.value(v) for ndn_i, v in scenario._mpisppy_model.rho.items()} + xbarsdict = {gd["nonant_names"][ndn_i]:\ + pyo.value(v) for ndn_i, v in scenario._mpisppy_model.xbars.items()} + + # Get data from host model + nonants_coeffs = gd["nonants_coeffs"] + host_model = scenario._mpisppy_model + W_on = host_model.W_on.value + prox_on = host_model.prox_on.value + # Update x coeff and x^2 coeff + for i, nonant in gd["nonants"].items(): # (Root, 1) : Nonant + new_coeff_val_xvar = nonants_coeffs[i] + W_on * (Wdict[nonant.VarName]) - prox_on * (rhodict[nonant.VarName] * xbarsdict[nonant.VarName]) + new_coeff_val_xsq = prox_on * rhodict[nonant.VarName]/2.0 + # Gurobipy does not seem to have setters/getters, instead we use attributes + nonant.Obj = new_coeff_val_xvar + gd["nonant_sqs"][i].Obj = new_coeff_val_xsq + + gs.update() + else: + pass # presumably an xhatter; we should check, I suppose + +def _copy_nonants_from_host(s): + gd = s._agnostic_dict + for ndn_i, gxvar in gd["nonants"].items(): + hostVar = s._mpisppy_data.nonant_indices[ndn_i] + guestVar = gd["nonants"][ndn_i] + if guestVar.LB == guestVar.UB: + guestVar.LB = gd["nonant_bounds"][ndn_i][0] + guestVar.UB = gd["nonant_bounds"][ndn_i][1] + if hostVar.is_fixed(): + guestVar.LB = hostVar._value + guestVar.UB = hostVar._value + else: + guestVar.Start = hostVar._value + +def _restore_nonants(Ag, s=None): + _copy_nonants_from_host(s) + +def _restore_original_fixedness(Ag, scenario): + _copy_nonants_from_host(scenario) + +def _fix_nonants(Ag, s=None): + _copy_nonants_from_host(s) + +def _fix_root_nonants(Ag, scenario): + _copy_nonants_from_host(scenario) diff --git a/examples/farmer/agnostic/farmer_pyomo_agnostic.py b/examples/farmer/agnostic/farmer_pyomo_agnostic.py new file mode 100644 index 00000000..1a423285 --- /dev/null +++ b/examples/farmer/agnostic/farmer_pyomo_agnostic.py @@ -0,0 +1,299 @@ +############################################################################### +# 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. +############################################################################### +# +# In this example, Pyomo is the guest language just for +# testing and documentation purposed. +""" +For other guest languages, the corresponding module is +still written in Python, it just needs to interact +with the guest language +""" + +import pyomo.environ as pyo +from pyomo.opt import SolutionStatus, TerminationCondition +import farmer # the native farmer (makes a few things easy) + +# for debuggig +from mpisppy import MPI +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +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, but + but treat Pyomo as a guest language. + + 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 + """ + s = farmer.scenario_creator(scenario_name, use_integer, sense, crops_multiplier, + num_scens, seedoffset) + # In general, be sure to process variables in the same order has the guest does (so indexes match) + gd = { + "scenario": s, + "nonants": {("ROOT",i): v for i,v in enumerate(s.DevotedAcreage.values())}, + "nonant_fixedness": {("ROOT",i): v.is_fixed() for i,v in enumerate(s.DevotedAcreage.values())}, + "nonant_start": {("ROOT",i): v._value for i,v in enumerate(s.DevotedAcreage.values())}, + "nonant_names": {("ROOT",i): v.name for i, v in enumerate(s.DevotedAcreage.values())}, + "probability": "uniform", + "sense": pyo.minimize, + "BFs": None + } + return gd + +#========= +def scenario_names_creator(num_scens,start=None): + return farmer.scenario_names_creator(num_scens,start) + + +#========= +def inparser_adder(cfg): + farmer.inparser_adder(cfg) + + +#========= +def kw_creator(cfg): + # creates keywords for scenario creator + return farmer.kw_creator(cfg) + +# This is not needed for PH +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + return farmer.sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario, **scenario_creator_kwargs) + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + pass + # (the fct in farmer won't work because the Var names don't match) + #farmer.scenario_denouement(rank, scenario_name, scenario) + + + +################################################################################################## +# begin callouts +# NOTE: the callouts all take the Ag object as their first argument, mainly to see cfg if needed +# the function names correspond to function names in mpisppy + +def attach_Ws_and_prox(Ag, sname, scenario): + # this is farmer specific, so we know there is not a W already, e.g. + # Attach W's and prox to the guest scenario. + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + nonant_idx = list(scenario._agnostic_dict["nonants"].keys()) + gs.W = pyo.Param(nonant_idx, initialize=0.0, mutable=True) + gs.W_on = pyo.Param(initialize=0, mutable=True, within=pyo.Binary) + gs.prox_on = pyo.Param(initialize=0, mutable=True, within=pyo.Binary) + gs.rho = pyo.Param(nonant_idx, mutable=True, default=Ag.cfg.default_rho) + + +def _disable_prox(Ag, scenario): + scenario._agnostic_dict["scenario"].prox_on._value = 0 + + +def _disable_W(Ag, scenario): + scenario._agnostic_dict["scenario"].W_on._value = 0 + + +def _reenable_prox(Ag, scenario): + scenario._agnostic_dict["scenario"].prox_on._value = 1 + + +def _reenable_W(Ag, scenario): + scenario._agnostic_dict["scenario"].W_on._value = 1 + + +def attach_PH_to_objective(Ag, sname, scenario, add_duals, add_prox): + # Deal with prox linearization and approximation later, + # i.e., just do the quadratic version + + ### The host has xbars and computes without involving the guest language + ### xbars = scenario._mpisppy_model.xbars + ### but instead, we are going to make guest xbars like other guests + + + gd = scenario._agnostic_dict + gs = gd["scenario"] # guest scenario handle + nonant_idx = list(gd["nonants"].keys()) + objfct = gs.Total_Cost_Objective # we know this is farmer... + ph_term = 0 + gs.xbars = pyo.Param(nonant_idx, mutable=True) + # Dual term (weights W) + if add_duals: + gs.WExpr = pyo.Expression(expr= sum(gs.W[ndn_i] * xvar for ndn_i,xvar in gd["nonants"].items())) + ph_term += gs.W_on * gs.WExpr + + # Prox term (quadratic) + if (add_prox): + prox_expr = 0. + for ndn_i, xvar in gd["nonants"].items(): + # expand (x - xbar)**2 to (x**2 - 2*xbar*x + xbar**2) + # x**2 is the only qradratic term, which might be + # dealt with differently depending on user-set options + if xvar.is_binary(): + xvarsqrd = xvar + else: + xvarsqrd = xvar**2 + prox_expr += (gs.rho[ndn_i] / 2.0) * \ + (xvarsqrd - 2.0 * gs.xbars[ndn_i] * xvar + gs.xbars[ndn_i]**2) + gs.ProxExpr = pyo.Expression(expr=prox_expr) + ph_term += gs.prox_on * gs.ProxExpr + + if gd["sense"] == pyo.minimize: + objfct.expr += ph_term + elif gd["sense"] == pyo.maximize: + objfct.expr -= ph_term + else: + raise RuntimeError(f"Unknown sense {gd['sense'] =}") + + +def solve_one(Ag, s, solve_keyword_args, gripe, tee=False, need_solution=True): + # This needs to attach stuff to s (see solve_one in spopt.py) + # Solve the guest language version, then copy values to the host scenario + + # This function needs to W on the guest right before the solve + + # We need to operate on the guest scenario, not s; however, attach things to s (the host scenario) + # and copy to s. If you are working on a new guest, you should not have to edit the s side of things + + # To acommdate the solve_one call from xhat_eval.py, we need to attach the obj fct value to s + + _copy_Ws_xbars_rho_from_host(s) + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + + print(f" in _solve_one {global_rank =}") + if global_rank == 0: + print(f"{gs.W.pprint() =}") + print(f"{gs.xbars.pprint() =}") + solver_name = s._solver_plugin.name + solver = pyo.SolverFactory(solver_name) + if 'persistent' in solver_name: + raise RuntimeError("Persistent solvers are not currently supported in the farmer agnostic example.") + ###solver.set_instance(ef, symbolic_solver_labels=True) + ###solver.solve(tee=True) + else: + solver_exception = None + try: + results = solver.solve(gs, tee=tee, symbolic_solver_labels=True,load_solutions=False) + except Exception as e: + results = None + solver_exception = e + + if (results is None) or (len(results.solution) == 0) or \ + (results.solution(0).status == SolutionStatus.infeasible) or \ + (results.solver.termination_condition == TerminationCondition.infeasible) or \ + (results.solver.termination_condition == TerminationCondition.infeasibleOrUnbounded) or \ + (results.solver.termination_condition == TerminationCondition.unbounded): + + s._mpisppy_data.scenario_feasible = False + + if gripe: + print (f"Solve failed for scenario {s.name} on rank {global_rank}") + if results is not None: + print ("status=", results.solver.status) + print ("TerminationCondition=", + results.solver.termination_condition) + + if solver_exception is not None and need_solution: + raise solver_exception + + else: + s._mpisppy_data.scenario_feasible = True + if gd["sense"] == pyo.minimize: + s._mpisppy_data.outer_bound = results.Problem[0].Lower_bound + else: + s._mpisppy_data.outer_bound = results.Problem[0].Upper_bound + gs.solutions.load_from(results) + # copy the nonant x values from gs to s so mpisppy can use them in s + for ndn_i, gxvar in gd["nonants"].items(): + # courtesy check for staleness on the guest side before the copy + if not gxvar.fixed and gxvar.stale: + try: + float(pyo.value(gxvar)) + except: # noqa + raise RuntimeError( + f"Non-anticipative variable {gxvar.name} on scenario {s.name} " + "reported as stale. This usually means this variable " + "did not appear in any (active) components, and hence " + "was not communicated to the subproblem solver. ") + + s._mpisppy_data.nonant_indices[ndn_i]._value = gxvar._value + + # the next line ignore bundling + s._mpisppy_data._obj_from_agnostic = pyo.value(gs.Total_Cost_Objective) + + # TBD: deal with other aspects of bundling (see solve_one in spopt.py) + + +# local helper +def _copy_Ws_xbars_rho_from_host(s): + # This is an important function because it allows us to capture whatever the host did + # print(f" {s.name =}, {global_rank =}") + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + for ndn_i, gxvar in gd["nonants"].items(): + assert hasattr(s, "_mpisppy_model"),\ + f"what the heck!! no _mpisppy_model {s.name =} {global_rank =}" + if hasattr(s._mpisppy_model, "W"): + gs.W[ndn_i] = pyo.value(s._mpisppy_model.W[ndn_i]) + gs.rho[ndn_i] = pyo.value(s._mpisppy_model.rho[ndn_i]) + gs.xbars[ndn_i] = pyo.value(s._mpisppy_model.xbars[ndn_i]) + else: + # presumably an xhatter + pass + + +# local helper +def _copy_nonants_from_host(s): + # values and fixedness; + gd = s._agnostic_dict + for ndn_i, gxvar in gd["nonants"].items(): + hostVar = s._mpisppy_data.nonant_indices[ndn_i] + guestVar = gd["nonants"][ndn_i] + if guestVar.is_fixed(): + guestVar.fixed = False + if hostVar.is_fixed(): + guestVar.fix(hostVar._value) + else: + guestVar._value = hostVar._value + + +def _restore_nonants(Ag, s): + # the host has already restored + _copy_nonants_from_host(s) + + +def _restore_original_fixedness(Ag, s): + _copy_nonants_from_host(s) + + +def _fix_nonants(Ag, s): + _copy_nonants_from_host(s) + + +def _fix_root_nonants(Ag, s): + _copy_nonants_from_host(s) + + + diff --git a/examples/farmer/farmer.py b/examples/farmer/farmer.py index 2d620b68..8f40726d 100644 --- a/examples/farmer/farmer.py +++ b/examples/farmer/farmer.py @@ -86,6 +86,8 @@ def scenario_creator( #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( diff --git a/mpisppy/agnostic/__init__.py b/mpisppy/agnostic/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/mpisppy/agnostic/agnostic.py b/mpisppy/agnostic/agnostic.py new file mode 100644 index 00000000..f1b42556 --- /dev/null +++ b/mpisppy/agnostic/agnostic.py @@ -0,0 +1,128 @@ +############################################################################### +# 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. +############################################################################### +# Agnostic.py +# This software is distributed under the 3-clause BSD License. + +""" +notes by dlw: + - The cfg will happen to have an Agnostic object added to it so mpisppy code can find it; + however, this code does not care about that. + - If a function in mpisppy has two callouts (a rarity), then the kwargs need to distinguish. +""" + +import sys +import inspect +import pyomo.environ as pyo +from mpisppy.utils import sputils + + +#======================================== +class Agnostic(): + """ + Args: + module (python module): None or the module with the callout fcts, + scenario_creator and other helper fcts. + cfg (Config): controls + """ + + def __init__(self, module, cfg): + self.module = module + self.cfg = cfg + + + def callout_agnostic(self, kwargs): + """ callout from mpi-sppy for AML-agnostic support + Args: + cfg (Config): the field "AML_agnostic" might contain a module with callouts + kwargs (dict): the keyword args for the callout function (e.g. scenario) + Calls: + a callout function that presumably has side-effects + Returns: + True if the callout was done and False if not + Note: + Throws an error if the module exists, but the fct is missing + """ + + if self.module is not None: + fname = inspect.stack()[1][3] + fct = getattr(self.module, fname, None) + if fct is None: + raise RuntimeError(f"AML-agnostic module or object is missing function {fname}") + try: + fct(Ag=self, **kwargs) + except Exception as e: + print(f"ERROR: AML-agnostic module or object had an error when calling {fname}", file=sys.stderr) + raise e + return True + else: + return False + + + def scenario_creator(self, sname): + """ create scenario sname by calling guest language, then attach stuff + Args: + sname (str): the scenario name that usually ends with a number + Returns: + scenario (Pyomo concrete model): a skeletal pyomo model with + a lot attached to it. + Note: + The python function scenario_creator in the module needs to + return a dict that we will call gd. + gd["scenario"]: the guest language model handle + gd["nonants"]: dict [(ndn,i)]: guest language Var handle + gd["nonant_names"]: dict [(ndn,i)]: str with name of variable + gd["nonant_fixedness"]: dict [(ndn,i)]: indicator of fixed variable + gd["nonant_start"]: dict [(ndn,i)]: float with starting value + gd["probability"]: float prob or str "uniform" + gd["obj_fct"]: the objective function from the guest + gd["sense"]: pyo.minimize or pyo.maximize + gd["BFs"]: scenario tree branching factors list or None + (for two stage models, the only value of ndn is "ROOT"; + i in (ndn, i) is always just an index) + """ + + crfct = getattr(self.module, "scenario_creator", None) + if crfct is None: + raise RuntimeError(f"AML-agnostic module {self.module.__name__} is missing the scenario_creator function") + kwfct = getattr(self.module, "kw_creator", None) + if kwfct is not None: + kwargs = kwfct(self.cfg) + gd = crfct(sname, **kwargs) + else: + gd = crfct(sname) + s = pyo.ConcreteModel(sname) + + ndns = [ndn for (ndn,i) in gd["nonants"].keys()] + iis = [i for (ndn,i) in gd["nonants"].keys()] # is is reserved... + s.nonantVars = pyo.Var(ndns, iis) + for idx,v in s.nonantVars.items(): + v._value = gd["nonant_start"][idx] + v.fixed = gd["nonant_fixedness"][idx] + + # we don't really need an objective, but we do need a sense + # note that other code may put W's and prox's on it + s.Obj = pyo.Objective(expr=0, sense=gd["sense"]) + s._agnostic_dict = gd + + assert gd["BFs"] is None, "We are only doing two stage for now" + # (it would not be that hard to be multi-stage; see hydro.py) + + sputils.attach_root_node(s, s.Obj, [s.nonantVars]) + + s._mpisppy_probability = gd["probability"] + + return s + + +############################################################################################################ + + +if __name__ == "__main__": + # For use by developers doing ad hoc testing + print("no main") diff --git a/mpisppy/agnostic/agnostic_cylinders.py b/mpisppy/agnostic/agnostic_cylinders.py new file mode 100644 index 00000000..7345421f --- /dev/null +++ b/mpisppy/agnostic/agnostic_cylinders.py @@ -0,0 +1,169 @@ +############################################################################### +# 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. +############################################################################### +""" +We need to get the module from the command line, then construct +the X_guest (e.g. Pyomo_guest, AMPL_guest) class to wrap the module +and feed it to the Agnostic object. + +I think it can be a pretty general cylinders program? +""" +import sys +from mpisppy.spin_the_wheel import WheelSpinner +import mpisppy.utils.cfg_vanilla as vanilla +import mpisppy.utils.config as config +import mpisppy.agnostic.agnostic as agnostic +import mpisppy.utils.sputils as sputils + +def _setup_args(m): + # m is the model file module + # NOTE: try to avoid adding features here that are not supported for agnostic + cfg = config.Config() + cfg.add_to_config(name="module_name", + description="file name that has scenario creator, etc.", + domain=str, + default=None, + argparse=True) + assert hasattr(m, "inparser_adder"), "The model file must have an inparser_adder function" + m.inparser_adder(cfg) + cfg.add_to_config(name="guest_language", + description="The language in which the modle is written (e.g. Pyomo or AMPL)", + domain=str, + default="None", + argparse=True) + cfg.add_to_config(name="ampl_model_file", + description="The .m file needed if the language is AMPL", + domain=str, + default=None, + argparse=True) + cfg.add_to_config(name="gams_model_file", + description="The original .gms file needed if the language is GAMS", + domain=str, + default=None, + argparse=True) + cfg.add_to_config(name="solution_base_name", + description="The string used fo a directory of ouput along with a csv and an npv file (default None, which means no soltion output)", + domain=str, + default=None) + 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() + return cfg + +def _parse_args(m): + # m is the model file module + cfg = _setup_args(m) + cfg.parse_command_line("agnostic cylinders") + return cfg + +def main(model_fname, module, cfg): + """ main is outside __main__ for testing + + Args: + model_fname (str): the module name from the command line) + module (Python module): from model_fname (redundant) + cfg (Pyomo config object): parsed arguments + """ + + supported_guests = {"Pyomo", "AMPL", "GAMS"} + # special hack to support bundles + if hasattr(module, "bundle_hack"): + module.bundle_hack(cfg) + # num_scens is now really numbuns + if cfg.guest_language not in supported_guests: + raise ValueError(f"Not a supported guest language: {cfg.guest_language}\n" + f" supported guests: {supported_guests}") + if cfg.guest_language == "Pyomo": + # now I need the pyomo_guest wrapper, then feed that to agnostic + from mpisppy.agnostic.pyomo_guest import Pyomo_guest + pg = Pyomo_guest(model_fname) + Ag = agnostic.Agnostic(pg, cfg) + + elif cfg.guest_language == "AMPL": + assert cfg.ampl_model_file is not None, "If the guest language is AMPL, you need ampl-model-file" + from mpisppy.agnostic.ampl_guest import AMPL_guest + guest = AMPL_guest(model_fname, cfg.ampl_model_file) + Ag = agnostic.Agnostic(guest, cfg) + + elif cfg.guest_language == "GAMS": + from mpisppy import MPI + fullcomm = MPI.COMM_WORLD + global_rank = fullcomm.Get_rank() + import mpisppy.agnostic.gams_guest as gams_guest + original_file_path = cfg.gams_model_file + new_file_path = gams_guest.file_name_creator(original_file_path) + nonants_name_pairs = module.nonants_name_pairs_creator() + if global_rank == 0: + gams_guest.create_ph_model(original_file_path, new_file_path, nonants_name_pairs) + print("Global rank 0 has created the new .gms model file") + fullcomm.Barrier() + + guest = gams_guest.GAMS_guest(model_fname, new_file_path, nonants_name_pairs, cfg) + Ag = agnostic.Agnostic(guest, cfg) + + scenario_creator = Ag.scenario_creator + assert hasattr(module, "scenario_denouement"), "The model file must have a scenario_denouement function" + scenario_denouement = module.scenario_denouement # should we go though Ag? + # note that if you are bundling, cfg.num_scens will be a fib (numbuns) + all_scenario_names = module.scenario_names_creator(cfg.num_scens) + + # Things needed for vanilla cylinders + beans = (cfg, scenario_creator, scenario_denouement, all_scenario_names) + + # Vanilla PH hub + hub_dict = vanilla.ph_hub(*beans, + scenario_creator_kwargs=None, # kwargs in Ag not here + ph_extensions=None, + ph_converger=None, + rho_setter = None) + # pass the Ag object via options... + hub_dict["opt_kwargs"]["options"]["Ag"] = Ag + + # xhat shuffle bound spoke + if cfg.xhatshuffle: + xhatshuffle_spoke = vanilla.xhatshuffle_spoke(*beans, scenario_creator_kwargs=None) + xhatshuffle_spoke["opt_kwargs"]["options"]["Ag"] = Ag + if cfg.lagrangian: + lagrangian_spoke = vanilla.lagrangian_spoke(*beans, scenario_creator_kwargs=None) + lagrangian_spoke["opt_kwargs"]["options"]["Ag"] = Ag + + list_of_spoke_dict = list() + if cfg.xhatshuffle: + list_of_spoke_dict.append(xhatshuffle_spoke) + if cfg.lagrangian: + list_of_spoke_dict.append(lagrangian_spoke) + + wheel = WheelSpinner(hub_dict, list_of_spoke_dict) + wheel.spin() + + if cfg.solution_base_name is not None: + wheel.write_first_stage_solution(f'{cfg.solution_base_name}.csv') + wheel.write_first_stage_solution(f'{cfg.solution_base_name}.npy', + first_stage_solution_writer=sputils.first_stage_nonant_npy_serializer) + wheel.write_tree_solution(f'{cfg.solution_base_name}') + + +if __name__ == "__main__": + if len(sys.argv) == 1: + print("need the python model file module name (no .py)") + print("usage, e.g.: python -m mpi4py agnostic_cylinders.py --module-name farmer4agnostic" --help) + quit() + + model_fname = sys.argv[2] + + module = sputils.module_name_to_module(model_fname) + + cfg = _parse_args(module) + + main(model_fname, module, cfg) diff --git a/mpisppy/agnostic/ampl_guest.py b/mpisppy/agnostic/ampl_guest.py new file mode 100644 index 00000000..19518507 --- /dev/null +++ b/mpisppy/agnostic/ampl_guest.py @@ -0,0 +1,369 @@ +############################################################################### +# 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. +############################################################################### +# This code sits between the guest model file wrapper and mpi-sppy +# AMPL is the guest language. Started by DLW June 2024 +""" +The guest model file (not this file) provides a scenario creator in Python +that attaches to each scenario a scenario probability (or "uniform") +and the following items to populate the guest dict (aka gd): + name + conditional probability + stage number + nonant var list + +(As of June 2024, we are going to be two-stage only...) + +The guest model file (which is in Python) also needs to provide: + scenario_creator_kwargs + scenario_names_creator + scenario_denouement + +The signature for the scenario creator in the Python model file for an +AMPL guest is the scenario_name, an ampl model file name, and then +keyword args. It is up to the function to instantiate the model +in the guest language and to make sure the data is correct for the +given scenario name and kwargs. + +For AMPL, the _nonant_varadata_list should contain objects obtained +from something like the get_variable method of an AMPL model object. + +Not concerning indexes: +To keep it simple and completely generic, we are throwing away a lot +of index information. This will probably slow down instantantion of +the objective fuction (see attach_PH_to_objective) +See farmer_ample_agnostic.py for an example where indexes are retained. + +""" + +import mpisppy.utils.sputils as sputils +import pyomo.environ as pyo + +from mpisppy import MPI # for debuggig +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +class AMPL_guest(): + """ + Provide an interface to a model file for an AMPL guest. + + Args: + model_file_name (str): name of Python file that has functions like scenario_creator + ampl_file_name (str): name of AMPL file that is passed to the model file + """ + def __init__(self, model_file_name, ampl_file_name): + self.model_file_name = model_file_name + self.model_module = sputils.module_name_to_module(model_file_name) + self.ampl_file_name = ampl_file_name + + + def scenario_creator(self, scenario_name, **kwargs): + """ Wrap the guest (AMPL in this case) scenario creator + + Args: + scenario_name (str): + Name of the scenario to construct. + """ + def _has_ints(s): + for _,v in s.getVariables(): + if "binary" in str(v) or "integer" in str(v): + return True + return False + + s, prob, nonant_vardata_list, obj_fct = self.model_module.scenario_creator(scenario_name, + self.ampl_file_name, + **kwargs) + if len(nonant_vardata_list) == 0: + raise RuntimeError(f"model file {self.model_file_name} has an empty " + f" nonant_vardata_list for {scenario_name =}") + # In general, be sure to process variables in the same order has the guest does (so indexes match) + nonant_vars = nonant_vardata_list # typing aid + def _vname(v): + return v[1].name().split('[')[0] + gd = { + "scenario": s, + "nonants": {("ROOT",i): v[1] for i,v in enumerate(nonant_vars)}, + "nonant_fixedness": {("ROOT",i): v[1].astatus()=="fixed" for i,v in enumerate(nonant_vars)}, + "nonant_start": {("ROOT",i): v[1].value() for i,v in enumerate(nonant_vars)}, + "nonant_names": {("ROOT",i): (_vname(v), v[0]) for i, v in enumerate(nonant_vars)}, + "probability": prob, + "obj_fct": obj_fct, + "sense": pyo.minimize, + "BFs": None, + "has_ints": _has_ints(s), + } + ##?xxxxx ? create nonant vars and put them on the ampl model, + ##?xxxxx create constraints to make them equal to the original nonants + + return gd + + + #========= + def scenario_names_creator(self, num_scens,start=None): + return self.model_module.scenario_names_creator(num_scens,start) + + + #========= + def inparser_adder(self, cfg): + self.model_module.inparser_adder(cfg) + + + #========= + def kw_creator(self, cfg): + # creates keywords for scenario creator + return self.model_module.kw_creator(cfg) + + # This is not needed for PH + def sample_tree_scen_creator(self, sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + return self.model_module.sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario, **scenario_creator_kwargs) + + #============================ + def scenario_denouement(self, rank, scenario_name, scenario): + pass + # (the fct in farmer won't work because the Var names don't match) + #self.model_module.scenario_denouement(rank, scenario_name, scenario) + + + ############################################################################ + # begin callouts + # NOTE: the callouts all take the Ag object as their first argument, mainly to see cfg if needed + # the function names correspond to function names in mpisppy + + def attach_Ws_and_prox(self, Ag, sname, scenario): + # Attach W's and rho to the guest scenario (mutable params). + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gd = scenario._agnostic_dict + nonants = gd["nonants"] + # (there must be some way to create and assign *mutable* params in on call to AMPL) + gs.eval("param W_on;") + gs.eval("let W_on := 0;") + gs.eval("param prox_on;") + gs.eval("let prox_on := 0;") + # we are trusting the order to match the nonant indexes + #create a nonant_indices set in AMPL + # This should exactly match the nonant_vars.keys() + gs.eval(f"set nonant_indices := {{0..{len(nonants)}-1}};") + # TBD xxxxx check for W on the model already + gs.eval("param W{nonant_indices};") + # Note: we should probably use set_values instead of let + gs.eval("let {i in nonant_indices} W[i] := 0;") + # start with rho at zero, but update before solve + # TBD xxxxx check for rho on the model already + gs.eval("param rho{nonant_indices};") + gs.eval("let {i in nonant_indices} rho[i] := 0;") + + def _disable_prox(self, Ag, scenario): + scenario._agnostic_dict["scenario"].get_parameter("prox_on").set(0) + + + def _disable_W(self, Ag, scenario): + scenario._agnostic_dict["scenario"].get_parameter("W_on").set(0) + + + def _reenable_prox(self, Ag, scenario): + scenario._agnostic_dict["scenario"].get_parameter("prox_on").set(1) + + + def _reenable_W(self, Ag, scenario): + scenario._agnostic_dict["scenario"].get_parameter("W_on").set(1) + + + def attach_PH_to_objective(self, Ag, sname, scenario, add_duals, add_prox): + # TBD: Deal with prox linearization and approximation later, + # i.e., just do the quadratic version + # Assume that nonant_indices is on the AMPL model + + # The host has xbars and computes without involving the guest language + + def _vname(i): + vtuple = gd['nonant_names'][('ROOT',i)] + return f"{vtuple[0]}" if vtuple[1] == "" else f"{vtuple[0]}['{vtuple[1]}']" + + gd = scenario._agnostic_dict + gs = gd["scenario"] # guest scenario handle + gs.eval("param xbars{nonant_indices};") + obj_fct = gd["obj_fct"] + objstr = str(obj_fct) + assert objstr.split (' ')[0] == "minimize", "We currently assume minimization" + + # Dual term (weights W) (This is where indexes are an issue) + phobjstr = "" + if add_duals: + phobjstr += "W_on * (" + for i in range(len(gd["nonants"])): + vname = _vname(i) + phobjstr += f"W[{i}] * {vname} + " + phobjstr = phobjstr[:-3] + ")" + # Prox term (quadratic) + ####### _see copy_nonants_from_host + if add_prox: + phobjstr += " + prox_on * (" + for i in range(len(gd["nonants"])): + vname = _vname(i) + phobjstr += f"(rho[{i}]/2.0) * ({vname} * {vname} - 2.0 * xbars[{i}] * {vname} + xbars[{i}]^2) + " + phobjstr = phobjstr[:-3] + ")" + objstr = objstr[:-1] + "+ (" + phobjstr + ");" + objparts = objstr.split() + objname = objparts[1] # has the colon, too + objstr = objstr.replace(f"minimize {objname}", "minimize phobj:") + obj_fct.drop() + gs.eval(objstr) + gs.eval("delete minus_profit;") + currentobj = gs.get_current_objective() + # see _copy_Ws_... see also the gams version + WParamDatas = list(gs.get_parameter("W").instances()) + xbarsParamDatas = list(gs.get_parameter("xbars").instances()) + rhoParamDatas = list(gs.get_parameter("rho").instances()) + gd["PH"] = { + "W": {("ROOT",i): v for i,v in enumerate(WParamDatas)}, + "xbars": {("ROOT",i): v for i,v in enumerate(xbarsParamDatas)}, + "rho": {("ROOT",i): v for i,v in enumerate(rhoParamDatas)}, + "obj": currentobj, + } + + + def solve_one(self, Ag, s, solve_keyword_args, gripe, tee=False, need_solution=True): + # This needs to attach stuff to s (see solve_one in spopt.py) + # Solve the guest language version, then copy values to the host scenario + + # This function needs to update W on the guest right before the solve + + # We need to operate on the guest scenario, not s; however, attach things to s (the host scenario) + # and copy to s. If you are working on a new guest, you should not have to edit the s side of things + + # To acommdate the solve_one call from xhat_eval.py, we need to attach the obj fct value to s + + self._copy_Ws_xbars_rho_from_host(s) + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + #### start debugging + if global_rank == 0 and False: + WParamDatas = list(gs.get_parameter("W").instances()) + print(f" in _solve_one {WParamDatas =} {global_rank =}") + xbarsParamDatas = list(gs.get_parameter("xbars").instances()) + print(f" in _solve_one {xbarsParamDatas =} {global_rank =}") + rhoParamDatas = list(gs.get_parameter("rho").instances()) + print(f" in _solve_one {rhoParamDatas =} {global_rank =}") + #### stop debugging + + solver_name = s._solver_plugin.name + gs.set_option("solver", solver_name) + if 'persistent' in solver_name: + raise RuntimeError("Persistent solvers are not currently supported in AMPL agnostic.") + gs.set_option("presolve", 0) + + solver_exception = None + try: + gs.solve() + except Exception as e: + solver_exception = e + + if gs.solve_result != "solved": + s._mpisppy_data.scenario_feasible = False + if gripe: + print (f"Solve failed for scenario {s.name} on rank {global_rank}") + print(f"{gs.solve_result =}") + s._mpisppy_data._obj_from_agnostic = None + return + + else: + s._mpisppy_data.scenario_feasible = True + + if solver_exception is not None and need_solution: + raise solver_exception + + + # For AMPL mips, we need to use the gap option to compute bounds + # https://amplmp.readthedocs.io/rst/features-guide.html + # As of Aug 2024, this is not tested... + mipgap = gs.getValue('_mipgap') if gd["has_ints"] else 0 + objobj = gs.get_current_objective() # different for xhatters + objval = objobj.value() + if gd["sense"] == pyo.minimize: + s._mpisppy_data.outer_bound = objval - mipgap + else: + s._mpisppy_data.inner_bound = objval + mipgap + + # copy the nonant x values from gs to s so mpisppy can use them in s + # in general, we need more checks (see the pyomo agnostic guest example) + for ndn_i, gxvar in gd["nonants"].items(): + try: # not sure this is needed + float(gxvar.value()) + except: # noqa + raise RuntimeError( + f"Non-anticipative variable {gxvar.name} on scenario {s.name} " + "had no value. This usually means this variable " + "did not appear in any (active) components, and hence " + "was not communicated to the subproblem solver. ") + if gxvar.astatus() == "pre": + raise RuntimeError( + f"Non-anticipative variable {gxvar.name} on scenario {s.name} " + "was presolved out. This usually means this variable " + "did not appear in any (active) components, and hence " + "was not communicated to the subproblem solver. ") + + s._mpisppy_data.nonant_indices[ndn_i]._value = gxvar.value() + + s._mpisppy_data._obj_from_agnostic = objval + + + # local helper + def _copy_Ws_xbars_rho_from_host(self, s): + # This is an important function because it allows us to capture whatever the host did + # print(f" {s.name =}, {global_rank =}") + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + # AMPL params are tuples (index, value), which are immutable + # ndn_i is a tuple (name, index) + if hasattr(s._mpisppy_model, "W"): + Wdict = {ndn_i[1]: pyo.value(v) for ndn_i, v in s._mpisppy_model.W.items()} + gs.get_parameter("W").set_values(Wdict) + rhodict = {ndn_i[1]: pyo.value(v) for ndn_i, v in s._mpisppy_model.rho.items()} + gs.get_parameter("rho").set_values(rhodict) + xbarsdict = {ndn_i[1]: pyo.value(v) for ndn_i, v in s._mpisppy_model.xbars.items()} + gs.get_parameter("xbars").set_values(xbarsdict) + else: + pass # presumably an xhatter; we should check, I suppose + + + # local helper + def _copy_nonants_from_host(self, s): + # values and fixedness; + gd = s._agnostic_dict + for ndn_i, gxvar in gd["nonants"].items(): + hostVar = s._mpisppy_data.nonant_indices[ndn_i] + guestVar = gd["nonants"][ndn_i] + if guestVar.astatus() == "fixed": + guestVar.unfix() + if hostVar.is_fixed(): + guestVar.fix(hostVar._value) + else: + guestVar.set_value(hostVar._value) + + + def _restore_nonants(self, Ag, s): + # the host has already restored + self._copy_nonants_from_host(s) + + + def _restore_original_fixedness(self, Ag, s): + self._copy_nonants_from_host(s) + + + def _fix_nonants(self, Ag, s): + # We are assuming the host did the fixing + self._copy_nonants_from_host(s) + + + def _fix_root_nonants(self, Ag, s): + # We are assuming the host did the fixing + self._copy_nonants_from_host(s) + + diff --git a/mpisppy/agnostic/examples/__init__.py b/mpisppy/agnostic/examples/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/mpisppy/agnostic/examples/afew_agnostic.py b/mpisppy/agnostic/examples/afew_agnostic.py new file mode 100644 index 00000000..6634ba4d --- /dev/null +++ b/mpisppy/agnostic/examples/afew_agnostic.py @@ -0,0 +1,55 @@ +############################################################################### +# 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. +############################################################################### +# Run a few examples; dlw Nov 2024; user-unfriendly +# Assumes you run from the agnostic/examples directory. +# python run_agnostic.py +# python run_agnostic.py --oversubscribe + +import os +import sys + +pyomo_solver_name = "cplex_direct" +ampl_solver_name = "gurobi" +gams_solver_name = "cplex" + +# Use oversubscribe if your computer does not have enough cores. +# Don't use this unless you have to. +# (This may not be allowed on versions of mpiexec) +mpiexec_arg = "" # "--oversubscribe" +if len(sys.argv) == 2: + mpiexec_arg = sys.argv[1] + +badguys = list() + +def do_one(np, argstring): + runstring = f"mpiexec -np {np} {mpiexec_arg} python -m mpi4py ../agnostic_cylinders.py {argstring}" + print(runstring) + code = os.system(runstring) + if code != 0: + badguys.append(runstring) + + + +print("skipping Pyomo because it is not working on github due to cplex_direct versus cplexdirect") +#do_one(3, f"--module-name farmer4agnostic --default-rho 1 --num-scens 6 --solver-name {pyomo_solver_name} --guest-language Pyomo --max-iterations 5") + +do_one(3, f"--module-name mpisppy.agnostic.examples.farmer_ampl_model --default-rho 1 --num-scens 3 --solver-name {ampl_solver_name} --guest-language AMPL --ampl-model-file farmer.mod --lagrangian --xhatshuffle --max-iterations 5") + +do_one(3, f"--module-name mpisppy.agnostic.examples.steel_ampl_model --default-rho 1 --num-scens 3 --seed 17 --solver-name {ampl_solver_name} --guest-language AMPL --ampl-model-file steel.mod --ampl-data-file steel.dat --max-iterations 10 --rel-gap 0.01 --xhatshuffle --lagrangian --solution-base-name steel") + +do_one(3, f"--module-name mpisppy.agnostic.examples.farmer_gams_model --num-scens 3 --default-rho 1 --solver-name {gams_solver_name} --max-iterations=5 --rel-gap 0.01 --display-progress --guest-language GAMS --gams-model-file farmer_average.gms") + + +if len(badguys) > 0: + print("\nBad Guys:") + for i in badguys: + print(i) + sys.exit(1) +else: + print("\nAll OK.") diff --git a/mpisppy/agnostic/examples/ag_farmer_gams.bash b/mpisppy/agnostic/examples/ag_farmer_gams.bash new file mode 100644 index 00000000..bf6d41bd --- /dev/null +++ b/mpisppy/agnostic/examples/ag_farmer_gams.bash @@ -0,0 +1,14 @@ +#!/bin/bash + +SOLVERNAME=gurobi + +python -u ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.farmer_gams_model --num-scens 3 --default-rho 1 --solver-name $SOLVERNAME --max-iterations=5 --rel-gap 0.01 --display-progress --guest-language GAMS --gams-model-file farmer_average.gms + +echo "^^^^lagrangian only" +mpiexec -np 2 python -m mpi4py ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.farmer_gams_model --num-scens 3 --default-rho 0.5 --solver-name $SOLVERNAME --max-iterations=30 --lagrangian --rel-gap 0.01 --guest-language GAMS --gams-model-file farmer_average.gms + +echo "^^^^ xhat only" +mpiexec -np 2 python -m mpi4py ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.farmer_gams_model --num-scens 3 --default-rho 0.5 --solver-name $SOLVERNAME --max-iterations=30 --xhatshuffle --rel-gap 0.01 --guest-language GAMS --gams-model-file farmer_average.gms + +echo "^^^^ lagrangian and xhat" +mpiexec -np 3 python -m mpi4py ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.farmer_gams_model --num-scens 3 --default-rho 0.5 --solver-name $SOLVERNAME --max-iterations=30 --xhatshuffle --lagrangian --rel-gap 0.01 --guest-language GAMS --gams-model-file farmer_average.gms diff --git a/mpisppy/agnostic/examples/ag_transport_gams.bash b/mpisppy/agnostic/examples/ag_transport_gams.bash new file mode 100644 index 00000000..0d74cb38 --- /dev/null +++ b/mpisppy/agnostic/examples/ag_transport_gams.bash @@ -0,0 +1,10 @@ +#!/bin/bash + +SOLVERNAME=gurobi +IT=150 + +python ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.transport_gams_model --num-scens 3 --default-rho 1 --solver-name $SOLVERNAME --max-iterations=${IT} --rel-gap 0.01 --display-progress --guest-language GAMS --gams-model-file transport_average.gms + +mpiexec -np 3 python -m mpi4py ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.transport_gams_model --num-scens 3 --default-rho 0.1 --solver-name $SOLVERNAME --max-iterations=${IT} --xhatshuffle --lagrangian --rel-gap 0.01 --guest-language GAMS --gams-model-file transport_average.gms + + diff --git a/mpisppy/agnostic/examples/executing_single_gms_file.py b/mpisppy/agnostic/examples/executing_single_gms_file.py new file mode 100644 index 00000000..09551798 --- /dev/null +++ b/mpisppy/agnostic/examples/executing_single_gms_file.py @@ -0,0 +1,33 @@ +############################################################################### +# 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. +############################################################################### +import os +import gams +import gamspy_base + +# This can be useful to execute a single gams file and print its values + +this_dir = os.path.dirname(os.path.abspath(__file__)) + +gamspy_base_dir = gamspy_base.__path__[0] + +w = gams.GamsWorkspace(working_directory=this_dir, system_directory=gamspy_base_dir) + +model = w.add_job_from_file("transport_ef") + +model.run()#output=sys.stdout) + +obj_dict = {} +for record in model.out_db.get_variable('z_stoch'): + obj_dict[tuple(record.keys)] = record.get_level() + +x_dict = {} +for record in model.out_db.get_variable('x'): + x_dict[tuple(record.keys)] = record.get_level() +print(f"{x_dict=}, {obj_dict=}") +print(f"obj_val = {model.out_db.get_variable('z_average').find_record().get_level()}") diff --git a/mpisppy/agnostic/examples/farmer.mod b/mpisppy/agnostic/examples/farmer.mod new file mode 100644 index 00000000..5a5f147e --- /dev/null +++ b/mpisppy/agnostic/examples/farmer.mod @@ -0,0 +1,111 @@ +# The farmer's problem in AMPL +# +# Reference: +# John R. Birge and Francois Louveaux. Introduction to Stochastic Programming. +# +# AMPL coding by Victor Zverovich; ## modifed by dlw; now *minization* + +##function expectation; +##function random; + +##suffix stage IN; + +set Crops; + +##set Scen; +##param P{Scen}; # probabilities + +param TotalArea; # acre +param PlantingCost{Crops}; # $/acre +param SellingPrice{Crops}; # $/T +param ExcessSellingPrice; # $/T +param PurchasePrice{Crops}; # $/T +param MinRequirement{Crops}; # T +param BeetsQuota; # T + +# Area in acres devoted to crop c. +var area{c in Crops} >= 0; + +# Tons of crop c sold (at favourable price) under scenario s. +var sell{c in Crops} >= 0, suffix stage 2; + +# Tons of sugar beets sold in excess of the quota under scenario s. +var sell_excess >= 0, suffix stage 2; + +# Tons of crop c bought under scenario s +var buy{c in Crops} >= 0, suffix stage 2; + +# The random variable (parameter) representing the yield of crop c. +##var RandomYield{c in Crops}; +param RandomYield{c in Crops}; + +# Realizations of the yield of crop c. +##param Yield{c in Crops, s in Scen}; # T/acre + +##maximize profit: +## expectation( +## ExcessSellingPrice * sell_excess + +## sum{c in Crops} (SellingPrice[c] * sell[c] - +## PurchasePrice[c] * buy[c])) - +## sum{c in Crops} PlantingCost[c] * area[c]; + +minimize minus_profit: + - ExcessSellingPrice * sell_excess - + sum{c in Crops} (SellingPrice[c] * sell[c] - + PurchasePrice[c] * buy[c]) + + sum{c in Crops} (PlantingCost[c] * area[c]); + +s.t. totalArea: sum {c in Crops} area[c] <= TotalArea; + +s.t. requirement{c in Crops}: + RandomYield[c] * area[c] - sell[c] + buy[c] >= MinRequirement[c]; + +s.t. quota: sell['beets'] <= BeetsQuota; + +s.t. sellBeets: + sell['beets'] + sell_excess <= RandomYield['beets'] * area['beets']; + +##yield: random({c in Crops} (RandomYield[c], {s in Scen} Yield[c, s])); + +data; + +set Crops := wheat corn beets; +#set Scen := below average above; + +param TotalArea := 500; + +##param Yield: +## below average above := +## wheat 2.0 2.5 3.0 +## corn 2.4 3.0 3.6 +## beets 16.0 20.0 24.0; + +# Average Scenario +param RandomYield := + wheat 2.5 + corn 3.0 + beets 20.0; + +param PlantingCost := + wheat 150 + corn 230 + beets 260; + +param SellingPrice := + wheat 170 + corn 150 + beets 36; + +param ExcessSellingPrice := 10; + +param PurchasePrice := + wheat 238 + corn 210 + beets 100; + +param MinRequirement := + wheat 200 + corn 240 + beets 0; + +param BeetsQuota := 6000; diff --git a/mpisppy/agnostic/examples/farmer.py b/mpisppy/agnostic/examples/farmer.py new file mode 100644 index 00000000..8f40726d --- /dev/null +++ b/mpisppy/agnostic/examples/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/mpisppy/agnostic/examples/farmer_ampl_model.py b/mpisppy/agnostic/examples/farmer_ampl_model.py new file mode 100644 index 00000000..6a38c839 --- /dev/null +++ b/mpisppy/agnostic/examples/farmer_ampl_model.py @@ -0,0 +1,110 @@ +############################################################################### +# 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. +############################################################################### +# In this example, AMPL is the guest language. +# This is the python model file for AMPL farmer. +# It will work with farmer.mod and slight deviations. + +from amplpy import AMPL +import pyomo.environ as pyo +import mpisppy.utils.sputils as sputils +import mpisppy.agnostic.examples.farmer as farmer +import numpy as np +from mpisppy import MPI # for debugging +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +# If you need random numbers, use this random stream: +farmerstream = np.random.RandomState() + + +# the first two args are in every scenario_creator for an AMPL model +def scenario_creator(scenario_name, ampl_file_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. + ampl_file_name (str): + The name of the ampl model file (with AMPL in it) + (This adds flexibility that maybe we don't need; it could be hardwired) + 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 + + NOTE: for ampl, the names will be tuples name, index + + Returns: + ampl_model (AMPL object): the AMPL model + prob (float or "uniform"): the scenario probability + nonant_var_data_list (list of AMPL variables): the nonants + obj_fct (AMPL Objective function): the objective function + """ + + assert crops_multiplier == 1, "for AMPL, just getting started with 3 crops" + + ampl = AMPL() + + ampl.read(ampl_file_name) + + # scenario specific data applied + scennum = sputils.extract_num(scenario_name) + assert scennum < 3, "three scenarios hardwired for now" + y = ampl.get_parameter("RandomYield") + if scennum == 0: # below + y.set_values({"wheat": 2.0, "corn": 2.4, "beets": 16.0}) + elif scennum == 2: # above + y.set_values({"wheat": 3.0, "corn": 3.6, "beets": 24.0}) + + areaVarDatas = list(ampl.get_variable("area").instances()) + + try: + obj_fct = ampl.get_objective("minus_profit") + except: + print("big troubles!!; we can't find the objective function") + raise + return ampl, "uniform", areaVarDatas, obj_fct + +#========= +def scenario_names_creator(num_scens,start=None): + return farmer.scenario_names_creator(num_scens,start) + + +#========= +def inparser_adder(cfg): + farmer.inparser_adder(cfg) + + +#========= +def kw_creator(cfg): + # creates keywords for scenario creator + return farmer.kw_creator(cfg) + +# This is not needed for PH +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + return farmer.sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario, **scenario_creator_kwargs) + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + pass + # (the fct in farmer won't work because the Var names don't match) + #farmer.scenario_denouement(rank, scenario_name, scenario) diff --git a/mpisppy/agnostic/examples/farmer_average.gms b/mpisppy/agnostic/examples/farmer_average.gms new file mode 100644 index 00000000..10af14a2 --- /dev/null +++ b/mpisppy/agnostic/examples/farmer_average.gms @@ -0,0 +1,91 @@ +$title The Farmer s Problem formulated for GAMS/DECIS (FARM,SEQ=199) + +$onText +This model helps a farmer to decide how to allocate +his or her land. The yields are uncertain. + + +Birge, R, and Louveaux, F V, Introduction to Stochastic Programming. +Springer, 1997. + +Keywords: linear programming, stochastic programming, agricultural cultivation, + farming, cropping +$offText + +*$if not set decisalg $set decisalg decism + +Set + crop / wheat, corn, sugarbeets / + cropr(crop) 'crops required for feeding cattle' / wheat, corn / + cropx / wheat + corn + beets1 'up to 6000 ton' + beets2 'in excess of 6000 ton' /; + +Parameter + yield(crop) 'tons per acre' / wheat 2.5 + corn 3 + sugarbeets 20 / + plantcost(crop) 'dollars per acre' / wheat 150 + corn 230 + sugarbeets 260 / + sellprice(cropx) 'dollars per ton' / wheat 170 + corn 150 + beets1 36 + beets2 10 / + purchprice(cropr) 'dollars per ton' / wheat 238 + corn 210 / + minreq(cropr) 'minimum requirements in ton' / wheat 200 + corn 240 /; + +Scalar + land 'available land' / 500 / + maxbeets1 'max allowed' / 6000 /; + +*-------------------------------------------------------------------------- +* First a non-stochastic version +*-------------------------------------------------------------------------- +Variable + x(crop) 'acres of land' + w(cropx) 'crops sold' + y(cropr) 'crops purchased' + yld(crop) 'yield' + profit 'objective variable'; + +Positive Variable x, w, y; + +Equation + profitdef 'objective function' + landuse 'capacity' + req(cropr) 'crop requirements for cattle feed' + ylddef 'calc yields' + beets 'total beet production'; + +$onText +The YLD variable and YLDDEF equation isolate the stochastic +YIELD parameter into one equation, making the DECIS setup +somewhat easier than if we would substitute YLD out of +the model. +$offText + +profitdef.. profit =e= - sum(crop, plantcost(crop)*x(crop)) + - sum(cropr, purchprice(cropr)*y(cropr)) + + sum(cropx, sellprice(cropx)*w(cropx)); + +landuse.. sum(crop, x(crop)) =l= land; + +ylddef(crop).. yld(crop) =e= yield(crop)*x(crop); + +req(cropr).. yld(cropr) + y(cropr) - sum(sameas(cropx,cropr),w(cropx)) =g= minreq(cropr); + +beets.. w('beets1') + w('beets2') =l= yld('sugarbeets'); + +x.up(crop) = land; +w.up('beets1') = maxbeets1; +$onText +__InsertPH__here_Model_defined_three_lines_later +$offText + +Model simple / profitdef, landuse, req, beets, ylddef /; + +solve simple using lp maximizing profit; \ No newline at end of file diff --git a/mpisppy/agnostic/examples/farmer_gams_model.py b/mpisppy/agnostic/examples/farmer_gams_model.py new file mode 100644 index 00000000..e70bac14 --- /dev/null +++ b/mpisppy/agnostic/examples/farmer_gams_model.py @@ -0,0 +1,103 @@ +############################################################################### +# 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. +############################################################################### +import os +import gamspy_base +import mpisppy.utils.sputils as sputils +import mpisppy.agnostic.farmer4agnostic as farmer + +from mpisppy import MPI # for debugging +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +LINEARIZED = True +this_dir = os.path.dirname(os.path.abspath(__file__)) +gamspy_base_dir = gamspy_base.__path__[0] + + +def nonants_name_pairs_creator(): + """Mustn't take any argument. Is called in agnostic cylinders + + Returns: + list of pairs (str, str): for each non-anticipative variable, the name of the support set must be given with the name of the variable. + If the set is a cartesian set, there should be no paranthesis when given + """ + return [("crop", "x")] + + +def stoch_param_name_pairs_creator(): + """ + Returns: + list of pairs (str, str): for each stochastic parameter, the name of the support set must be given with the name of the parameter. + If the set is a cartesian set, there should be no paranthesis when given + """ + return [("crop", "yield")] + + +def scenario_creator(scenario_name, mi, job, cfg=None): + """ Create a scenario for the (scalable) farmer example. + + Args: + scenario_name (str): + Name of the scenario to construct. + mi (gams model instance): the base model + job (gams job) : not used for farmer + cfg: pyomo config + """ + ### This part is model specific, we define the values of the stochastic parameters depending on scenario_name + scennum = sputils.extract_num(scenario_name) + assert scennum < 3, "three scenarios hardwired for now" + y = mi.sync_db.get_parameter("yield") + + if scennum == 0: # below + y.add_record("wheat").value = 2.0 + y.add_record("corn").value = 2.4 + y.add_record("sugarbeets").value = 16.0 + elif scennum == 1: # average + y.add_record("wheat").value = 2.5 + y.add_record("corn").value = 3.0 + y.add_record("sugarbeets").value = 20.0 + elif scennum == 2: # above + y.add_record("wheat").value = 3.0 + y.add_record("corn").value = 3.6 + y.add_record("sugarbeets").value = 24.0 + + return mi + + +#========= +def scenario_names_creator(num_scens,start=None): + return farmer.scenario_names_creator(num_scens,start) + + +#========= +def inparser_adder(cfg): + farmer.inparser_adder(cfg) + + +#========= +def kw_creator(cfg): + # creates keywords for scenario creator + #kwargs = farmer.kw_creator(cfg) + kwargs = {} + kwargs["cfg"] = cfg + #kwargs["nonants_name_pairs"] = nonants_name_pairs_creator() + return kwargs + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + # doesn't seem to be called + if global_rank == 1: + x_dict = {} + for x_record in scenario._agnostic_dict["scenario"].sync_db.get_variable('x'): + x_dict[x_record.get_keys()[0]] = x_record.get_level() + print(f"In {scenario_name}: {x_dict}") + pass + # (the fct in farmer won't work because the Var names don't match) + #farmer.scenario_denouement(rank, scenario_name, scenario) diff --git a/mpisppy/agnostic/examples/farmer_gurobipy_model.py b/mpisppy/agnostic/examples/farmer_gurobipy_model.py new file mode 100644 index 00000000..84945f87 --- /dev/null +++ b/mpisppy/agnostic/examples/farmer_gurobipy_model.py @@ -0,0 +1,264 @@ +############################################################################### +# 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. +############################################################################### +# In this example, Gurobipy is the guest language +import gurobipy as gp +from gurobipy import GRB +import pyomo.environ as pyo +import examples.farmer.farmer as farmer +import example.farmer.farmer_gurobipy_mod as farmer_gurobipy_mod +import numpy as np + +from mpisppy import MPI # for debugging +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +farmerstream = np.random.RandomState() + + +def scenario_creator(scenario_name, use_integer=False, sense=GRB.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): + gurobipy sense (minimization or maximization). Must be either + GRB.MINIMIZE or GRB.MAXIMIZE. Default is GRB.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 + + NOTE: + Most other people wont have the model file already written + """ + gd = farmer_gurobipy_mod.scenario_creator(scenario_name, use_integer, sense, crops_multiplier) + + return gd + +#========== +def scenario_names_creator(num_scens,start=None): + return farmer.scenario_names_creator(num_scens,start) + +#========== +def inparser_adder(cfg): + return farmer.inparser_adder(cfg) + +#========== +def kw_creator(cfg): + # creates keywords for scenario creator + return farmer.kw_creator(cfg) + +# This is not needed for PH +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + return farmer.sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario, **scenario_creator_kwargs) + +def scenario_denouement(rank, scenario_name, scenario): + pass + # (the fct in farmer won't work because the Var names don't match) + #farmer.scenario_denouement(rank, scenario_name, scenario) + +################################################################################################## +# begin callouts +# NOTE: the callouts all take the Ag object as their first argument, mainly to see cfg if needed +# the function names correspond to function names in mpisppy +def attach_Ws_and_prox(Ag, sname, scenario): + + # this is gurobipy farmer specific, so we know there is not a W already + # Gurobipy is special so we are gonna need to maintain the coeffs ourselves + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gd = scenario._agnostic_dict + obj = gs.getObjective() + + obj_func_terms = [obj.getVar(i) for i in range(obj.size())] + nonant_coeffs = {} + nonants_not_in_obj = {} + + # Check to see if the nonants are in the objective function + for i, nonant in gd["nonants"].items(): + found_nonant_in_obj = False + for obj_func_term in obj_func_terms: + if obj_func_term.sameAs(nonant): + nonant_coeffs[nonant] = nonant.Obj + found_nonant_in_obj = True + break + if not found_nonant_in_obj: + print('We dont have a coeff for this nonant which is bad so we will default to 0') + nonant_coeffs[nonant] = 0 + nonants_not_in_obj[i] = nonant + + # Update/attach nonants' coeffs to the dictionary + gd["nonant_coeffs"] = nonant_coeffs + gd["nonants_not_in_obj"] = nonants_not_in_obj + + +def _disable_prox(Ag, scenario): + pass + # raise RuntimeError("Did not expect _disable_prox") + +def _disable_W(Ag, scenario): + pass + # raise RuntimeError("Did not expect _disable_W") + +def _reenable_prox(Ag, scenario): + pass + # raise RuntimeError("Did not expect _reenable_prox") + +def _reenable_W(Ag, scenario): + pass + # raise RuntimeError("Did not expect _reenable_W") + +def attach_PH_to_objective(Ag, sname, scenario, add_duals, add_prox): + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + gd = scenario._agnostic_dict + nonant_coeffs = gd["nonant_coeffs"] + + ''' + At this point we can assume that all that all the nonants are already in the objective function from attach_Ws_and_prox + but for the nonants that are not in the objective function we can just attach it to the objective funciton with thier respective coefficients + ''' + + # Adding all the nonants into the obj function + obj = gs.getObjective() + for nonant_not_in_obj in gd["nonants_not_in_obj"]: + obj += (nonant_coeffs[nonant_not_in_obj] * nonant_not_in_obj) + + # At this point all nonants are in the obj + nonant_sqs = {} + # Need to create a new var for x^2, and then add a constraint to the var with x val, so that we can set coeff value later on + for i, nonant in gd["nonants"].items(): + # Create a constaint that sets x * x = xsq + nonant_sq = gs.addVar(vtype=GRB.CONTINUOUS, obj=nonant.Obj**2, name=f"{nonant.VarName}sq") + gs.addConstr(nonant * nonant == nonant_sq, f'{nonant.VarName}sqconstr') + # Put the x^2 in the objective function + obj += nonant_sq + nonant_sqs[i] = nonant_sq + + # Update model and gd + gs.update() + gd["nonant_sqs"] = nonant_sqs + + _copy_Ws_xbars_rho_from_host(scenario) + +def solve_one(Ag, s, solve_keyword_args, gripe, tee, need_solution=True): + _copy_Ws_xbars_rho_from_host(s) + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + + solver_exception = None + # Assuming gs is a Gurobi model, we can start solving + try: + gs.optimize() + except gp.GurobiError as e: + s._mpisppy_data.scenario_feasible = False + solver_exception = e + if gripe: + print(f"Solve failed with error {str(e)} for scenario {s.name}") + return + + # TBD: what about not optimal and need_solution? + if solver_exception is not None and need_solution: + raise solver_exception + + if gs.status != gp.GRB.Status.OPTIMAL: + s._mpisppy_data.scenario_feasible = False + if gripe: + print(f"Solve failed for scenario {s.name}") + return + + s._mpisppy_data.scenario_feasible = True + + # Objective value extraction + objval = gs.getObjective().getValue() + + if gd["sense"] == gp.GRB.MINIMIZE: + s._mpisppy_data.outer_bound = objval + else: + s._mpisppy_data.outer_bound = objval + + # Copy the non-anticipative variable values from guest to host scenario + for ndn_i, gxvar in gd["nonants"].items(): + grb_var = gs.getVarByName(gxvar.VarName) + if grb_var is None: + raise RuntimeError( + f"Non-anticipative variable {gxvar.varname} on scenario {s.name} " + "was not found in the Gurobi model." + ) + s._mpisppy_data.nonant_indices[ndn_i]._value = grb_var.X + + # Store the objective function value in the host scenario + s._mpisppy_data._obj_from_agnostic = objval + + # Additional checks and operations for bundling if needed (depending on the problem) + # ... + +def _copy_Ws_xbars_rho_from_host(scenario): + # Calculates the coefficients of the new expanded objective function + # Regardless need to calculate coefficients for x^2 and x + gd = scenario._agnostic_dict + gs = scenario._agnostic_dict["scenario"] # guest handle + + # Decide if we are using PH or xhatter + if hasattr(scenario._mpisppy_model, "W"): + # Get our Ws, rhos, and xbars + Wdict = {gd["nonant_names"][ndn_i]:\ + pyo.value(v) for ndn_i, v in scenario._mpisppy_model.W.items()} + rhodict = {gd["nonant_names"][ndn_i]:\ + pyo.value(v) for ndn_i, v in scenario._mpisppy_model.rho.items()} + xbarsdict = {gd["nonant_names"][ndn_i]:\ + pyo.value(v) for ndn_i, v in scenario._mpisppy_model.xbars.items()} + + # Get data from host model + nonants_coeffs = gd["nonants_coeffs"] + host_model = scenario._mpisppy_model + W_on = host_model.W_on.value + prox_on = host_model.prox_on.value + # Update x coeff and x^2 coeff + for i, nonant in gd["nonants"].items(): # (Root, 1) : Nonant + new_coeff_val_xvar = nonants_coeffs[i] + W_on * (Wdict[nonant.VarName] * nonants_coeffs[i]) - prox_on * (rhodict[nonant.VarName] * xbarsdict[nonant.VarName]) + new_coeff_val_xsq = prox_on * rhodict[nonant.VarName]/2.0 + nonant.Obj = new_coeff_val_xvar + gd["nonant_sqs"][i].Obj = new_coeff_val_xsq + + gs.update() + else: + pass # presumably an xhatter; we should check, I suppose + +def _copy_nonants_from_host(s): + gd = s._agnostic_dict + for ndn_i, gxvar in gd["nonants"].items(): + hostVar = s._mpisppy_data.nonant_indices[ndn_i] + guestVar = gd["nonants"][ndn_i] + if guestVar.lb == guestVar.ub: + guestVar.lb = 0 + guestVar.ub = float("inf") + if hostVar.is_fixed(): + guestVar.lb = hostVar._value + guestVar.ub = hostVar._value + else: + guestVar.Start = hostVar._value + +def _restore_nonants(Ag, s=None): + _copy_nonants_from_host(s) + +def _restore_original_fixedness(Ag, scenario): + _copy_nonants_from_host(scenario) + +def _fix_nonants(Ag, s=None): + _copy_nonants_from_host(s) + +def _fix_root_nonants(Ag, scenario): + _copy_nonants_from_host(scenario) diff --git a/mpisppy/agnostic/examples/go.bash b/mpisppy/agnostic/examples/go.bash new file mode 100644 index 00000000..ca1cd428 --- /dev/null +++ b/mpisppy/agnostic/examples/go.bash @@ -0,0 +1,12 @@ +#!/bin/bash +python ../agnostic_cylinders.py --module-name farmer4agnostic --default-rho 1 --num-scens 6 --solver-name cplex --guest-language Pyomo --max-iterations 1 + +python ../agnostic_cylinders.py --module-name farmer4agnostic --default-rho 1 --num-scens 6 --solver-name cplex --guest-language Pyomo --max-iterations 1 --bundle-size 3 + +mpiexec -np 3 python ../agnostic_cylinders.py --module-name farmer4agnostic --default-rho 1 --num-scens 6 --solver-name cplex --guest-language Pyomo --xhatshuffle --lagrangian --max-iterations 10 --rel-gap .01 + +mpiexec -np 3 python ../agnostic_cylinders.py --module-name farmer4agnostic --default-rho 1 --num-scens 6 --solver-name cplex --guest-language Pyomo --bundle-size 3 --xhatshuffle --lagrangian --max-iterations 10 --rel-gap .01 + +#python ../agnostic_cylinders.py --module-name farmer4agnostic --default-rho 1 --num-scens 3 --solver-name cplex --guest-language Pyomo +# NOTE: you need the AMPL solvers!!! +#python ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.farmer_ampl_model --default-rho 1 --num-scens 3 --solver-name gurobi --guest-language AMPL --ampl-model-file farmer.mod diff --git a/mpisppy/agnostic/examples/steel.bash b/mpisppy/agnostic/examples/steel.bash new file mode 100644 index 00000000..b41a064b --- /dev/null +++ b/mpisppy/agnostic/examples/steel.bash @@ -0,0 +1,7 @@ +#!/bin/bash +#python ../agnostic_cylinders.py --module-name farmer4agnostic --default-rho 1 --num-scens 3 --solver-name cplex --guest-language Pyomo +# NOTE: you need the AMPL solvers!!! +#python -u ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.steel_ampl_model --default-rho 1 --num-scens 3 --seed 17 --solver-name gurobi --guest-language AMPL --ampl-model-file steel.mod --ampl-data-file "$ampl_data_file" --max-iterations 10 +#could i ad a -- data file name +ampl_data_file="steel.dat" +mpiexec -np 3 python -u ../agnostic_cylinders.py --module-name mpisppy.agnostic.examples.steel_ampl_model --default-rho 1 --num-scens 3 --seed 17 --solver-name gurobi --guest-language AMPL --ampl-model-file steel.mod --ampl-data-file "$ampl_data_file" --max-iterations 10 --rel-gap 0.01 --xhatshuffle --lagrangian --solution-base-name steel diff --git a/mpisppy/agnostic/examples/steel.dat b/mpisppy/agnostic/examples/steel.dat new file mode 100644 index 00000000..da0c1754 --- /dev/null +++ b/mpisppy/agnostic/examples/steel.dat @@ -0,0 +1,7 @@ +data; +param avail; +set PROD := bands coils; +param: rate profit market := + bands 200 25 6000 + coils 140 30 4000 ; +param avail := 40; \ No newline at end of file diff --git a/mpisppy/agnostic/examples/steel.mod b/mpisppy/agnostic/examples/steel.mod new file mode 100644 index 00000000..3c8212f6 --- /dev/null +++ b/mpisppy/agnostic/examples/steel.mod @@ -0,0 +1,18 @@ +set PROD; # products + +param rate {PROD} > 0; # tons produced per hour +param avail >= 0; # hours available in week + +param profit {PROD}; # profit per ton +param market {PROD} >= 0; # limit on tons sold in week + +var Make {p in PROD} >= 0, <= market[p]; # tons produced + +minimize minus_profit: 0 - sum {p in PROD} profit[p] * Make[p]; + + # Objective: total profits from all products + +subject to Time: sum {p in PROD} (1/rate[p]) * Make[p] <= avail; + + # Constraint: total of hours used by all + # products may not exceed hours available \ No newline at end of file diff --git a/mpisppy/agnostic/examples/steel_ampl_model.py b/mpisppy/agnostic/examples/steel_ampl_model.py new file mode 100644 index 00000000..e8c2e079 --- /dev/null +++ b/mpisppy/agnostic/examples/steel_ampl_model.py @@ -0,0 +1,104 @@ +############################################################################### +# 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. +############################################################################### +# In this example, AMPL is the guest language. +# This is the python model file for AMPL farmer. +# It will work with farmer.mod and slight deviations. + +from amplpy import AMPL +import mpisppy.utils.sputils as sputils +import numpy as np + +from mpisppy import MPI # for debugging +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +# If you need random numbers, use this random stream: +steelstream = np.random.RandomState() + +# the first two args are in every scenario_creator for an AMPL model +ampl = AMPL() + +def scenario_creator(scenario_name, ampl_file_name, cfg=None): + """" + NOTE: for ampl, the names will be tuples name, index + + Returns: + ampl_model (AMPL object): the AMPL model + prob (float or "uniform"): the scenario probability + nonant_var_data_list (list of AMPL variables): the nonants + obj_fct (AMPL Objective function): the objective function + tbd finish doc string + need to add the args after completed + """ + + + ampl = AMPL() + ampl.read(ampl_file_name) + ampl.read_data(cfg.ampl_data_file) + scennum = sputils.extract_num(scenario_name) + seedoffset = cfg.seed + #steelstream.seed(scennum+seedoffset) + np.random.seed(scennum + seedoffset) + + # RANDOMIZE THE DATA****** + products = ampl.get_set("PROD") + ppu = ampl.get_parameter("profit") + if np.random.uniform() 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 seedoffset + # as a kwarg to scenario_creator + 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 a varlist, which is used to create a vardata list + # (This list needs to whatever the guest needs, not what Pyomo needs) + varlist = [model.DevotedAcreage] + model._nonant_vardata_list = sputils.build_vardatalist(model, varlist) + sputils.attach_root_node(model, 0, 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 + + elif "bund" == scenario_name[:4] or "Bund" == scenario_name[:4]: + firstnum = int(scenario_name.split("_")[1]) + lastnum = int(scenario_name.split("_")[2]) + assert (lastnum-firstnum+1) == bunsize + assert num_scens % bunsize != 0, "Due to laziness, we need equal sized bundels" + snames = [f"scen{i}" for i in range(firstnum, lastnum+1)] + + bunkwargs = {"use_integer": use_integer, + "sense": sense, + "crops_multiplier": crops_multiplier, + "num_scens":None} + bunkwargs["seedoffset"] = seedoffset + firstnum + + # it is easy to make the EF in Pyomo; see create_EF + # Note that it call scenario_creator, but this time it will be + # with scenario names. + bundle = sputils.create_EF(snames, scenario_creator, + scenario_creator_kwargs=bunkwargs, + EF_name=scenario_name, + nonant_for_fixed_vars = False) + # It simplifies things if we assume that it is a 2-stage problem, + # or that the bundles consume entire second stage nodes, + # then all we need is a root node and the only nonants that need to be reported are + # at the root node (otherwise, more coding is required here to figure out which nodes and Vars + # are shared with other bundles) + # Note: farmer is 2 stage. + nonantlist = [v for idx,v in bundle.ref_vars.items() if idx[0] =="ROOT"] + sputils.attach_root_node(bundle, 0, nonantlist) + # scenarios are equally likely so bundles are too + bundle._mpisppy_probability = 1/numbuns + return bundle + else: + raise RuntimeError (f"Scenario name does not have scen or bund: {scenario_name}") + +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 helper functions +#========= +def scenario_names_creator(num_scens, start=None): + # 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 + if bunsize == 0: + return [f"scen{i}" for i in range(start,start+num_scens)] + else: + # The hack should have changed the value of num_scens to be a fib! + # We will assume that start and and num_scens refers to bundle counts. + # Bundle numbers are zero based and scenario numbers as well. + return [f"bundle_{i*bunsize}_{(i+1)*bunsize-1}" 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) + cfg.add_to_config("bundle_size", + description="number of scenarios per bundle (default 0, which means no bundles, as does 1)", + domain=int, + default=0) + + +#========= +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 helper functions + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + sname = scenario_name + #print("denouement needs work") + #scenario.pprint() + return + 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)) + + +# special helper function (hack) for bundles +def bundle_hack(cfg): + # Hack to put bundle information in global variables to be used by + # the names creator. (only relevant for bundles) + # numbuns and bunsize are globals with default value 0 + if cfg.bundle_size > 1: + assert cfg.num_scens % cfg.bundle_size == 0,\ + "Due to laziness, the bundle size must divide the number of scenarios" + global bunsize, numbuns, original_num_scens + bunsize = cfg.bundle_size + numbuns = cfg.num_scens // cfg.bundle_size + original_num_scens = cfg.num_scens + cfg.num_scens = numbuns + + diff --git a/mpisppy/agnostic/gams_guest.py b/mpisppy/agnostic/gams_guest.py new file mode 100644 index 00000000..829d6ad7 --- /dev/null +++ b/mpisppy/agnostic/gams_guest.py @@ -0,0 +1,624 @@ +############################################################################### +# 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. +############################################################################### +# In this example, GAMS is the guest language. +# NOTE: unlike everywhere else, we are using xbar instead of xbars (no biggy) + +""" +This file tries to show many ways to do things in gams, +but not necessarily the best ways in any case. +""" + +import itertools +import os +import gams +import gamspy_base +import shutil +import tempfile +import re + +import pyomo.environ as pyo +import mpisppy.utils.sputils as sputils + +from mpisppy import MPI # for debugging +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +LINEARIZED = True +GM_NAME = "PH___Model" # to put in the new GAMS files +this_dir = os.path.dirname(os.path.abspath(__file__)) +gamspy_base_dir = gamspy_base.__path__[0] + +class GAMS_guest(): + """ + Provide an interface to a model file for a GAMS guest. + + Args: + model_file_name (str): name of Python file that has functions like scenario_creator + gams_file_name (str): name of GAMS file that is passed to the model file + nonants_name_pairs (list of (str,str)): list of (non_ant_support_set_name, non_ant_variable_name) + cfg (pyomo config object) + """ + def __init__(self, model_file_name, new_file_name, nonants_name_pairs, cfg): + self.model_file_name = model_file_name + self.model_module = sputils.module_name_to_module(model_file_name) + self.new_file_name = new_file_name + self.nonants_name_pairs = nonants_name_pairs + self.cfg = cfg + self.temp_dir = tempfile.TemporaryDirectory() + + def __del__(self): + self.temp_dir.cleanup() + + def scenario_creator(self, scenario_name, **kwargs): + """ Wrap the guest (GAMS in this case) scenario creator + + Args: + scenario_name (str): + Name of the scenario to construct. + + """ + new_file_name = self.new_file_name + assert new_file_name is not None + stoch_param_name_pairs = self.model_module.stoch_param_name_pairs_creator() + + ws = gams.GamsWorkspace(working_directory=self.temp_dir.name, system_directory=gamspy_base_dir) + + ### Calling this function is required regardless of the model + # This function creates a model instance not instantiated yet, and gathers in glist all the parameters and variables that need to be modifiable + mi, job, glist, all_ph_parameters_dicts, xlo_dict, xup_dict, x_out_dict = pre_instantiation_for_PH(ws, new_file_name, self.nonants_name_pairs, stoch_param_name_pairs) + + opt = ws.add_options() + opt.all_model_types = self.cfg.solver_name + if LINEARIZED: + mi.instantiate(f"{GM_NAME} using lp minimizing objective_ph", glist, opt) + else: + mi.instantiate(f"{GM_NAME} using qcp minimizing objective_ph", glist, opt) + + ### Calling this function is required regardless of the model + # This functions initializes, by adding records (and values), all the parameters that appear due to PH + nonant_set_sync_dict = adding_record_for_PH(self.nonants_name_pairs, self.cfg, all_ph_parameters_dicts, xlo_dict, xup_dict, x_out_dict, job) + + mi = self.model_module.scenario_creator(scenario_name, mi, job, **kwargs) + mi.solve() + nonant_variable_list = [nonant_var for (_, nonant_variables_name) in self.nonants_name_pairs for nonant_var in mi.sync_db.get_variable(nonant_variables_name)] + + gd = { + "scenario": mi, + "nonants": {("ROOT",i): 0 for i,v in enumerate(nonant_variable_list)}, + "nonant_fixedness": {("ROOT",i): v.get_lower() == v.get_upper() for i,v in enumerate(nonant_variable_list)}, + "nonant_start": {("ROOT",i): v.get_level() for i,v in enumerate(nonant_variable_list)}, + "probability": "uniform", + "sense": pyo.minimize, + "BFs": None, + "nonants_name_pairs": self.nonants_name_pairs, + "nonant_set_sync_dict": nonant_set_sync_dict + } + return gd + + + #========= + def scenario_names_creator(self, num_scens,start=None): + return self.model_module.scenario_names_creator(num_scens,start) + + + #========= + def inparser_adder(self, cfg): + self.model_module.inparser_adder(cfg) + + + #========= + def kw_creator(self, cfg): + # creates keywords for scenario creator + return self.model_module.kw_creator(cfg) + + + # This is not needed for PH + def sample_tree_scen_creator(self, sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + return self.model_module.sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario, **scenario_creator_kwargs) + + #============================ + def scenario_denouement(self, rank, scenario_name, scenario): + pass + # (the fct in farmer won't work because the Var names don't match) + #self.model_module.scenario_denouement(rank, scenario_name, scenario) + + ################################################################################################## + # begin callouts + # NOTE: the callouts all take the Ag object as their first argument, mainly to see cfg if needed + # the function names correspond to function names in mpisppy + + def attach_Ws_and_prox(self, Ag, sname, scenario): + # Done in create_ph_model + pass + + + def _disable_prox(self, Ag, scenario): + #print(f"In {global_rank=} for {scenario.name}: disabling prox") + scenario._agnostic_dict["scenario"].sync_db.get_parameter("prox_on").first_record().value = 0 + + + def _disable_W(self, Ag, scenario): + #print(f"In {global_rank=} for {scenario.name}: disabling W") + scenario._agnostic_dict["scenario"].sync_db.get_parameter("W_on").first_record().value = 0 + + + def _reenable_prox(self, Ag, scenario): + #print(f"In {global_rank=} for {scenario.name}: reenabling prox") + scenario._agnostic_dict["scenario"].sync_db.get_parameter("prox_on").first_record().value = 1 + + + def _reenable_W(self, Ag, scenario): + #print(f"In {global_rank=} for {scenario.name}: reenabling W") + scenario._agnostic_dict["scenario"].sync_db.get_parameter("W_on").first_record().value = 1 + + + def attach_PH_to_objective(self, Ag, sname, scenario, add_duals, add_prox): + # Done in create_ph_model + pass + + def solve_one(self, Ag, s, solve_keyword_args, gripe, tee, need_solution=True): + # s is the host scenario + # This needs to attach stuff to s (see solve_one in spopt.py) + # Solve the guest language version, then copy values to the host scenario + + # This function needs to put W on the guest right before the solve + + # We need to operate on the guest scenario, not s; however, attach things to s (the host scenario) + # and copy to s. If you are working on a new guest, you should not have to edit the s side of things + + # To acommdate the solve_one call from xhat_eval.py, we need to attach the obj fct value to s + self._copy_Ws_xbar_rho_from_host(s) + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + + """before_solve_bounds = {} + W_dict = {} + for record in s._agnostic_dict["scenario"].sync_db.get_variable("x"): + before_solve_bounds[tuple(record.keys)] = record.get_level(), record.get_lower(), record.get_upper() + W_dict[tuple(record.keys)] = s._agnostic_dict["scenario"].sync_db.get_parameter("ph_W_x").find_record(record.keys).get_value() + print(f"For {s.name} in {global_rank=}: {before_solve_bounds=} \n and {W_dict=}")""" + + solver_exception = None + #import sys + #print(f"SOLVING FOR {s.name}") + try: + gs.solve() + #gs.solve(output=sys.stdout)#update_type=2) + except Exception as e: + solver_exception = e + print(f"{solver_exception=}") + + solve_ok = (1, 2, 7, 8, 15, 16, 17) + + #print(f"{gs.model_status=}, {gs.solver_status=}") + if gs.model_status not in solve_ok: + s._mpisppy_data.scenario_feasible = False + if gripe: + print (f"Solve failed for scenario {s.name} on rank {global_rank}") + print(f"{gs.model_status =}") + s._mpisppy_data._obj_from_agnostic = None + return + + if solver_exception is not None and need_solution: + raise solver_exception + + s._mpisppy_data.scenario_feasible = True + + # For debugging, thisprints the weights of the PH model + """W_dict = {} + for record in s._agnostic_dict["scenario"].sync_db.get_variable("x"): + #after_solve_bounds[tuple(record.keys)] = record.get_level(), record.get_lower(), record.get_upper() + W_dict[tuple(record.keys)] = s._agnostic_dict["scenario"].sync_db.get_parameter("ph_W_x").find_record(record.keys).get_value() + print(f"For {s.name} in {global_rank=}: {W_dict=}")""" + + # For debugging + """if global_rank == 0: + x_dict = {} + for record in s._agnostic_dict["scenario"].sync_db.get_variable("x"): + #xbar = s._agnostic_dict["scenario"].sync_db.get_parameter("xbar").find_record(record.keys).get_value() + x_var_rec = s._agnostic_dict["scenario"].sync_db.get_variable("x").find_record(record.keys) + #x = x_var_rec.get_level() + x = x_var_rec.get_level() #, x_var_rec.get_lower() == x_var_rec.get_level(), x_var_rec.get_upper() == x_var_rec.get_level(), x_var_rec.get_lower()==s._agnostic_dict["scenario"].sync_db.get_parameter("xlo").find_record(record.keys).get_value(), x_var_rec.get_upper()==s._agnostic_dict["scenario"].sync_db.get_parameter("xup").find_record(record.keys).get_value() + #assert xbar == x, f"for {s.name}, {tuple(record.keys)}: {xbar=} and {x=}" + x_dict[tuple(record.keys)] = x + print(f"for {s.name}: {x_dict=}")""" + + objval = gs.sync_db.get_variable('objective_ph').find_record().get_level() + + # copy the nonant x values from gs to s so mpisppy can use them in s + # in general, we need more checks (see the pyomo agnostic guest example) + i = 0 + for nonants_set, nonants_var in gd["nonants_name_pairs"]: + for record in gs.sync_db.get_variable(nonants_var): + ndn_i = ('ROOT', i) + s._mpisppy_data.nonant_indices[ndn_i]._value = record.get_level() + i += 1 + + if gd["sense"] == pyo.minimize: + s._mpisppy_data.outer_bound = objval + else: + s._mpisppy_data.outer_bound = objval + + # the next line ignores bundling + s._mpisppy_data._obj_from_agnostic = objval + + # TBD: deal with other aspects of bundling (see solve_one in spopt.py) + #print(f"For {s.name} in {global_rank=}: {objval=}") + + + # local helper + def _copy_Ws_xbar_rho_from_host(self, s): + # special for farmer + # print(f" debug copy_Ws {s.name =}, {global_rank =}") + + gd = s._agnostic_dict + # could/should use set values, but then use a dict to get the indexes right + gs = gd["scenario"] + if hasattr(s._mpisppy_model, "W"): + i=0 + for nonants_set, nonants_var in gd["nonants_name_pairs"]: + #print(f'{gd["nonant_set_sync_dict"]}') + rho_dict = {} + for element in gd["nonant_set_sync_dict"][nonants_set]: + ndn_i = ('ROOT', i) + rho_dict[element] = s._mpisppy_model.rho[ndn_i].value + gs.sync_db[f"ph_W_{nonants_var}"].find_record(element).set_value(s._mpisppy_model.W[ndn_i].value) + gs.sync_db[f"rho_{nonants_var}"].find_record(element).set_value(s._mpisppy_model.rho[ndn_i].value) + gs.sync_db[f"{nonants_var}bar"].find_record(element).set_value(s._mpisppy_model.xbars[ndn_i].value) + i += 1 + #print(f"for {s.name} in {global_rank=}, {rho_dict=}") + + + # local helper + def _copy_nonants_from_host(self, s): + # values and fixedness; + #print(f"copiing nonants from host in {s.name}") + gs = s._agnostic_dict["scenario"] + gd = s._agnostic_dict + + i = 0 + for nonants_set, nonants_var in gd["nonants_name_pairs"]: + gs.sync_db.get_parameter(f"{nonants_var}lo").clear() + gs.sync_db.get_parameter(f"{nonants_var}up").clear() + for element in gd["nonant_set_sync_dict"][nonants_set]: + ndn_i = ("ROOT", i) + hostVar = s._mpisppy_data.nonant_indices[ndn_i] + if hostVar.is_fixed(): + #print(f"FIXING for {global_rank=}, in {s.name}, for {element=}: {hostVar._value=}") + gs.sync_db.get_parameter(f"{nonants_var}lo").add_record(element).set_value(hostVar._value) + gs.sync_db.get_parameter(f"{nonants_var}up").add_record(element).set_value(hostVar._value) + else: + gs.sync_db.get_variable(nonants_var).find_record(element).set_level(hostVar._value) + i += 1 + #print("LEAVING _copy_nonants") + + + def _restore_nonants(self, Ag, s): + # the host has already restored + self._copy_nonants_from_host(s) + + + def _restore_original_fixedness(self, Ag, s): + # the host has already restored + #This doesn't seem to be used and may not work correctly + self._copy_nonants_from_host(self, s) + + + def _fix_nonants(self, Ag, s): + # the host has already fixed? + self._copy_nonants_from_host(s) + + + def _fix_root_nonants(self, Ag, s): + #This doesn't seem to be used and may not work correctly + self._copy_nonants_from_host(s) + + +### This function creates a new gams model file including PH before anything else happens + +def create_ph_model(original_file_path, new_file_path, nonants_name_pairs): + """Copies the original file in original_file_path and creates a new file in new_file_path including the PH penalty, weights... + + Args: + original_file_path (str): + The path (including the name of the gms file) of the original file + new_file_path (str): + The path (including the name of the gms file) in which is created the gams model with the ph_objective + nonants_name_pairs (list of (str,str)): + List of (non_ant_support_set_name, non_ant_variable_name) + """ + # Copy the original file + shutil.copy2(original_file_path, new_file_path) + + # Read the content of the new file + with open(new_file_path, 'r') as file: + lines = file.readlines() + + # The keyword is used to insert everything from PH. + # For now the model was always defined three lines later to make life easier + insert_keyword = "__InsertPH__here_Model_defined_three_lines_later" + line_number = None + + # Searches through the lines for the keyword (from the end because its located in the end) + # Also captures whether the problem is a minimization or maximization + # problem and modifies the solve line (although it might not be necessary) + for i in range(len(lines)): + index = i # len(lines)-1-i + line = lines[index] + if line.startswith("Model"): + words = re.findall(r'\b\w+\b', line) + gmodel_name = words[1] + line = line.replace(gmodel_name, GM_NAME) + lines[index] = line + + if line.startswith("solve"): + # Should be in the last lines. This line + words = re.findall(r'\b\w+\b', line) + # print(f"{words=}") + if "minimizing" in words: + sense = "minimizing" + sign = "+" + elif "maximizing" in words: + print("WARNING: the objective function's sign has been changed") + sense = "maximizing" + sign = "-" + else: + raise RuntimeError(f"The line: {line}, doesn't include any sense") + assert gmodel_name == words[1], f"Model line as model name {gmodel_name} but solve line has {words[1]}" + # The word following the sense is the objective value + index_word = words.index(sense) + previous_objective = words[index_word + 1] + line = line.replace(sense, "minimizing") + line = line.replace(gmodel_name, GM_NAME) + new_obj = "objective_ph" + lines[index] = new_obj.join(line.rsplit(previous_objective, 1)) + + if insert_keyword in line: + line_number = index + + assert line_number is not None, "the insert_keyword is not used" + + # Where the text will be inserted. After the comment + insert_position = line_number + 2 + + #First modify the model to include the new equations and assert that the model is defined at the good position + model_line = lines[insert_position + 1] + model_line_stripped = model_line.strip().lower() + + ### Modify the line where the model is defined. If all the lines are included, nthing changes. Otherwise the new lines are added. + model_line_text = "" + if LINEARIZED: + for nonants_name_pair in nonants_name_pairs: + nonants_support_set, nonant_variables = nonants_name_pair + model_line_text += f", PenLeft_{nonant_variables}, PenRight_{nonant_variables}" + + assert "model" in model_line_stripped and "/" in model_line_stripped and model_line_stripped.endswith("/;"), "this is not the model line" + all_words = [" all ", "/all ", " all/", "/all/"] + all_model = False + for word in all_words: + if word in model_line: + all_model = True + if all_model: # we still use all the equations + lines[insert_position + 1] = model_line + else: # we have to specify which equations we use + lines[insert_position + 1] = model_line[:-4] + model_line_text + ", objective_ph_def" + model_line[-4:] + + ### Adds all the lines in the end of the gms file (where insert here is written). + # The lines contain the parameters, variables, equations... that are necessary for PH. Their value will be instanciated + # later in adding_record_for_PH + + parameter_definition = "" + scalar_definition = """ + W_on 'activate w term' / 0 / + prox_on 'activate prox term' / 0 /""" + variable_definition = "" + linearized_inequation_definition = "" + objective_ph_excess = "" + linearized_equation_expression = "" + parameter_initialization = "" + + for nonant_name_pair in nonants_name_pairs: + nonants_support_set, nonant_variables = nonant_name_pair + if "(" in nonants_support_set: + nonants_paranthesis_support_set = nonants_support_set + else: + nonants_paranthesis_support_set = f"({nonants_support_set})" + + parameter_definition += f""" + ph_W_{nonant_variables}({nonants_support_set}) 'ph weight' + {nonant_variables}bar({nonants_support_set}) 'ph average' + rho_{nonant_variables}({nonants_support_set}) 'ph rho'""" + + parameter_definition2 = f""" + ph_W_{nonant_variables}({nonants_support_set}) 'ph weight' /set.{nonants_support_set} 0/ + {nonant_variables}bar({nonants_support_set}) 'ph average' /set.{nonants_support_set} 0/ + rho_{nonant_variables}({nonants_support_set}) 'ph rho' /set.{nonants_support_set} 0/""" + + parameter_definition += f""" + {nonant_variables}up({nonants_support_set}) 'upper bound on {nonant_variables}' + {nonant_variables}lo({nonants_support_set}) 'lower bound on {nonant_variables}'""" + + parameter_definition2 += f""" + {nonant_variables}up({nonants_support_set}) 'upper bound on {nonant_variables}' /set.{nonants_support_set} 500/ + {nonant_variables}lo({nonants_support_set}) 'lower bound on {nonant_variables}' /set.{nonants_support_set} 0/""" + + variable_definition += f""" + PHpenalty_{nonant_variables}({nonants_support_set}) 'linearized prox penalty'""" + + if LINEARIZED: + linearized_inequation_definition += f""" + PenLeft_{nonant_variables}({nonants_support_set}) 'left side of linearized PH penalty' + PenRight_{nonant_variables}({nonants_support_set}) 'right side of linearized PH penalty'""" + + if LINEARIZED: + PHpenalty = f"PHpenalty_{nonant_variables}({nonants_support_set})" + else: + PHpenalty = f"({nonant_variables}({nonants_support_set}) - {nonant_variables}bar({nonants_support_set}))*({nonant_variables}({nonants_support_set}) - {nonant_variables}bar({nonants_support_set}))" + + # This line is the only line where a cartesian set should appear with paranthesis for the sum. That is why we use nonants_paranthesis_support_set + objective_ph_excess += f""" + + W_on * sum({nonants_paranthesis_support_set}, ph_W_{nonant_variables}({nonants_support_set})*{nonant_variables}({nonants_support_set})) + + prox_on * sum({nonants_paranthesis_support_set}, 0.5 * rho_{nonant_variables}({nonants_support_set}) * {PHpenalty})""" + + if LINEARIZED: + linearized_equation_expression += f""" +PenLeft_{nonant_variables}({nonants_support_set}).. PHpenalty_{nonant_variables}({nonants_support_set}) =g= ({nonant_variables}.up({nonants_support_set}) - {nonant_variables}bar({nonants_support_set})) * ({nonant_variables}({nonants_support_set}) - {nonant_variables}bar({nonants_support_set})); +PenRight_{nonant_variables}({nonants_support_set}).. PHpenalty_{nonant_variables}({nonants_support_set}) =g= ({nonant_variables}bar({nonants_support_set}) - {nonant_variables}.lo({nonants_support_set})) * ({nonant_variables}bar({nonants_support_set}) - {nonant_variables}({nonants_support_set})); +""" + parameter_initialization += f""" +ph_W_{nonant_variables}({nonants_support_set}) = 0; +{nonant_variables}bar({nonants_support_set}) = 0; +rho_{nonant_variables}({nonants_support_set}) = 0; +{nonant_variables}up({nonants_support_set}) = 0; +{nonant_variables}lo({nonants_support_set}) = 0; +""" + + my_text = f""" + +Parameter{parameter_definition}; + +Scalar{scalar_definition}; + +{parameter_initialization} + +Variable{variable_definition} + objective_ph 'final objective augmented with ph cost'; + +Equation{linearized_inequation_definition} + objective_ph_def 'defines objective_ph'; + +objective_ph_def.. objective_ph =e= {sign} {previous_objective} {objective_ph_excess}; + +{linearized_equation_expression} +""" + + lines.insert(insert_position, my_text) + + # Write the modified content back to the new file + with open(new_file_path, 'w') as file: + file.writelines(lines) + + #print(f"Modified file saved as: {new_file_path}") + + +def file_name_creator(original_file_path): + """ Returns the path and name of the gms file where ph will be added + Args: + original_file_path (str): the path (including the name) of the original gms path + """ + # Get the directory and filename + + directory, filename = os.path.split(os.path.abspath(original_file_path)) + name, ext = os.path.splitext(filename) + + assert ext == ".gms", "the original data file should be a gms file" + + # Create the new filename + if LINEARIZED: + new_filename = f"{name}_ph_linearized{ext}" + else: + print("WARNING: the normal quadratic PH has not been tested") + new_filename = f"{name}_ph_quadratic{ext}" + new_file_path = os.path.join(directory, new_filename) + return new_file_path + + +### Generic functions called inside the specific scenario creator +def _add_or_get_set(mi, out_set): + # Captures the set using data from the out_database. If it hasn't been added yet to the model insatnce it adds it as well + try: + return mi.sync_db.add_set(out_set.name, out_set._dim, out_set.text) + except gams.GamsException: + return mi.sync_db.get_set(out_set.name) + +def pre_instantiation_for_PH(ws, new_file_name, nonants_name_pairs, stoch_param_name_pairs): + """creates dictionaries indexed over the name of the stochastic parameters or non anticipative variables + so they can be added to the GamsModifier. Variables are obtained from the initial model database. + + Args: + ws (GamsWorkspace): the workspace to create the model instance + new_file_name (str): the gms file in which is created the gams model with the ph_objective + nonants_name_pairs (list of pairs (str, str)): for each non-anticipative variable, the name of the support set must be given with the name of the paramete + stoch_param_name_pairs (_type_): for each stochastic parameter, the name of the support set must be given with the name of the variable + + Returns: + tuple: include everything needed for creating the model instance + nonant_set_sync_dict gives the name of all the elements of the sets presented as tuples. It is useful if the set is a cartesian set: i,j then the elements + will be of the shape (element_in_i,element_in_j). Some functions iterate over this set. + + """ + ### First create the model instance + job = ws.add_job_from_file(new_file_name.replace(".gms","")) + cp = ws.add_checkpoint() + mi = cp.add_modelinstance() + + job.run(checkpoint=cp) # at this point the model (with what data?) is solved, it creates the file _gams_py_gjo0.lst + + ### Add to the elements that should be modified the stochastic parameters + # The parameters don't exist yet in the model instance, so they need to be redefined using the job + stoch_sets_out_dict = {param_name: [job.out_db.get_set(elementary_set) for elementary_set in set_name.split(",")] for set_name, param_name in stoch_param_name_pairs} + stoch_sets_sync_dict = {param_name: [_add_or_get_set(mi, out_elementary_set) for out_elementary_set in out_elementary_sets] for param_name, out_elementary_sets in stoch_sets_out_dict.items()} + glist = [gams.GamsModifier(mi.sync_db.add_parameter_dc(param_name, [sync_elementary_set for sync_elementary_set in sync_elementary_sets])) for param_name, sync_elementary_sets in stoch_sets_sync_dict.items()] + + ph_W_dict = {nonant_variables_name: mi.sync_db.add_parameter_dc(f"ph_W_{nonant_variables_name}", nonants_support_set_name.split(","), "ph weight") for nonants_support_set_name, nonant_variables_name in nonants_name_pairs} + xbar_dict = {nonant_variables_name: mi.sync_db.add_parameter_dc(f"{nonant_variables_name}bar", nonants_support_set_name.split(","), "xbar") for nonants_support_set_name, nonant_variables_name in nonants_name_pairs} + rho_dict = {nonant_variables_name: mi.sync_db.add_parameter_dc(f"rho_{nonant_variables_name}", nonants_support_set_name.split(","), "rho") for nonants_support_set_name, nonant_variables_name in nonants_name_pairs} + + # x_out is necessary to add the x variables to the database as we need the type and dimension of x + x_out_dict = {nonant_variables_name: job.out_db.get_variable(f"{nonant_variables_name}") for _, nonant_variables_name in nonants_name_pairs} + x_dict = {nonant_variables_name: mi.sync_db.add_variable(f"{nonant_variables_name}", x_out_dict[nonant_variables_name]._dim, x_out_dict[nonant_variables_name].vartype) for _, nonant_variables_name in nonants_name_pairs} + xlo_dict = {nonant_variables_name: mi.sync_db.add_parameter(f"{nonant_variables_name}lo", x_out_dict[nonant_variables_name]._dim, f"lower bound on {nonant_variables_name}") for _, nonant_variables_name in nonants_name_pairs} + xup_dict = {nonant_variables_name: mi.sync_db.add_parameter(f"{nonant_variables_name}up", x_out_dict[nonant_variables_name]._dim, f"upper bound on {nonant_variables_name}") for _, nonant_variables_name in nonants_name_pairs} + + W_on = mi.sync_db.add_parameter("W_on", 0, "activate w term") + prox_on = mi.sync_db.add_parameter("prox_on", 0, "activate prox term") + + glist += [gams.GamsModifier(ph_W_dict[nonants_name_pair[1]]) for nonants_name_pair in nonants_name_pairs] \ + + [gams.GamsModifier(xbar_dict[nonants_name_pair[1]]) for nonants_name_pair in nonants_name_pairs] \ + + [gams.GamsModifier(rho_dict[nonants_name_pair[1]]) for nonants_name_pair in nonants_name_pairs] \ + + [gams.GamsModifier(W_on)] \ + + [gams.GamsModifier(prox_on)] \ + + [gams.GamsModifier(x_dict[nonants_name_pair[1]], gams.UpdateAction.Lower, xlo_dict[nonants_name_pair[1]]) for nonants_name_pair in nonants_name_pairs] \ + + [gams.GamsModifier(x_dict[nonants_name_pair[1]], gams.UpdateAction.Upper, xup_dict[nonants_name_pair[1]]) for nonants_name_pair in nonants_name_pairs] + + all_ph_parameters_dicts = {"ph_W_dict": ph_W_dict, "xbar_dict": xbar_dict, "rho_dict": rho_dict, "W_on": W_on, "prox_on": prox_on} + return mi, job, glist, all_ph_parameters_dicts, xlo_dict, xup_dict, x_out_dict + + +def adding_record_for_PH(nonants_name_pairs, cfg, all_ph_parameters_dicts, xlo_dict, xup_dict, x_out_dict, job): + """This function adds records and repopulates the parameters in the gams instance. + + Args: + most arguments returned by pre_instanciation + + Returns: + nonant_set_sync_dict (dict{str: list of str tuples}): given a nonant_support_set_name (key), the value is the list + of all the elements of the sets presented as tuples. It is useful if the set is a cartesian set: i,j then the elements + will be of the shape (element_in_i,element_in_j). Some functions iterate over this set. + """ + + ### Gather the list of non-anticipative variables and their sets from the job, to modify them and add PH related parameters + nonant_sets_out_dict = {nonant_set_name: [job.out_db.get_set(elementary_set) for elementary_set in nonant_set_name.split(",")] for nonant_set_name, param_name in nonants_name_pairs} + + nonant_set_sync_dict = {nonant_set_name: [element for element in itertools.product(*[[rec.keys[0] for rec in out_elementary_set] for out_elementary_set in out_elementary_sets])] for nonant_set_name, out_elementary_sets in nonant_sets_out_dict.items()} + + for nonants_name_pair in nonants_name_pairs: + nonants_set_name, nonant_variables_name = nonants_name_pair + + for c in nonant_set_sync_dict[nonants_set_name]: + all_ph_parameters_dicts["ph_W_dict"][nonant_variables_name].add_record(c).value = 0 + all_ph_parameters_dicts["xbar_dict"][nonant_variables_name].add_record(c).value = 0 + all_ph_parameters_dicts["rho_dict"][nonant_variables_name].add_record(c).value = cfg.default_rho + xlo_dict[nonant_variables_name].add_record(c).value = x_out_dict[nonant_variables_name].find_record(c).lower + xup_dict[nonant_variables_name].add_record(c).value = x_out_dict[nonant_variables_name].find_record(c).upper + all_ph_parameters_dicts["W_on"].add_record().value = 0 + all_ph_parameters_dicts["prox_on"].add_record().value = 0 + return nonant_set_sync_dict diff --git a/mpisppy/agnostic/gurobipy_dilemma.md b/mpisppy/agnostic/gurobipy_dilemma.md new file mode 100644 index 00000000..fd6c4efc --- /dev/null +++ b/mpisppy/agnostic/gurobipy_dilemma.md @@ -0,0 +1,18 @@ +# +The following is the gurobipy dilemma + +Unlike other mathematical modeling languages like Pyomo, AMPL, and GAMS, Gurobipy serves as an API to the Gurobi solver and does not retain symbols. Symbols are essential for updating and modifying the model when new data is added. In Gurobipy, if changed need to be made to a constraint or an objective functio, we workaround with two approaches. + +**Recompute the Model**: Recomputing the model each time a modification occurs. This method is inefficient due to the repeated recomputation of the model. + +**Maintain Symbols Manually**: Since symbols are not retained in Gurobipy, we must manually maintain our version of the symbols by expanding the objective function expression and calculating the coefficients for each term (e.g., xx and x2x2). The useful feature in Gurobipy here is the ability to modify coefficients in linear constraints using the `changeCoeff(constraint, var, newVal)` function. This allows us to update the coefficients in linear constraints accordingly. + +The callout functions in `farmer_yyyy_agnostic.py` will behave diffrently. Instead of attaching elements to the model and changing those Param values, functions must now operate on an expanded version of the objective function. It's important to note that `changeCoeff` only modifies linear expressions, which poses a challenge since quadriatic expresssions need to be related to linear expressions; so we will need a new variable xs2 which is constrained to be equal to $x^2$ + +The diffrence between the callout functions: + - `attach_Ws_and_prox(Ag, sname, scenario)` will retrieve the current nonants coefficients from the objective function, and if the nonants aren't in the obj they are noted with 0 and put in a seperate list + - `attach_PH_to_obj()` - create list of x^2 vars and attach coefs to obj function, make sure all nonants are in the objective function, call `_copy_Ws_xbars_rho_from_host(scenario)`` + - `_copy_Ws_xbars_rho_from_host(scenario)` function that will be called to calculate the necessary coefficients for the objective function + - we should probably call the functions something else now + + diff --git a/mpisppy/agnostic/pyomo_guest.py b/mpisppy/agnostic/pyomo_guest.py new file mode 100644 index 00000000..3dae7948 --- /dev/null +++ b/mpisppy/agnostic/pyomo_guest.py @@ -0,0 +1,319 @@ +############################################################################### +# 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. +############################################################################### +# This code sits between the guest model file and mpi-sppy +# Pyomo is the guest language. Started by DLW April 2024 +""" +For other guest languages, the corresponding module is +still written in Python, it just needs to interact +with the guest language + +The guest model file (not this file) provides a scenario creator in the guest +language that attaches to each scenario a scenario probability (or "uniform") +and the following items to populate the guest dict (aka gd): + name + conditional probability + stage number + nonant var list + +The guest model file also needs to somehow (it might depend on the language) +provide hooks to: + scenario_creator_kwargs + scenario_names_creator + scenario_denouement + + Note: we already have a lot of two-stage models in Pyomo that would + be handy for testing. All that needs to be done, is to attach + the nonant varlist as _nonant_vars to the scenario when it is created. +""" +import mpisppy.utils.sputils as sputils +import pyomo.environ as pyo +from pyomo.opt import SolutionStatus, TerminationCondition + +from mpisppy import MPI # for debugging +fullcomm = MPI.COMM_WORLD +global_rank = fullcomm.Get_rank() + +class Pyomo_guest(): + """ + Provide an interface to a model file for a Pyomo guest + """ + def __init__(self, model_file_name): + self.model_module = sputils.module_name_to_module(model_file_name) + + + def scenario_creator(self, scenario_name, **kwargs): + """ Wrap the guest (Pyomo in this case) scenario creator + + Args: + scenario_name (str): + Name of the scenario to construct. + """ + s = self.model_module.scenario_creator(scenario_name, **kwargs) + ### TBD: assert that this is minimization? + if hasattr(s, "_nonant_vardata_list"): + nonant_vars = s._nonant_vardata_list # a list of vars + elif hasattr(s, "_mpisppy_node_list"): + assert len(s._mpisppy_node_list) == 1, "multi-stage agnostic with Pyomo as guest not yet supported." + nonant_vars = s._mpisppy_node_list[0].nonant_vardata_list + else: + raise RuntimeError("Scenario must have either _mpisppy_node_list or _nonant_vardata_list") + # In general, be sure to process variables in the same order has the guest does (so indexes match) + gd = { + "scenario": s, + "nonants": {("ROOT",i): v for i,v in enumerate(nonant_vars)}, + "nonant_fixedness": {("ROOT",i): v.is_fixed() for i,v in enumerate(nonant_vars)}, + "nonant_start": {("ROOT",i): v._value for i,v in enumerate(nonant_vars)}, + "nonant_names": {("ROOT",i): v.name for i, v in enumerate(nonant_vars)}, + "probability": s._mpisppy_probability, + "sense": pyo.minimize, + "BFs": None + } + # we don't need to attach nonants to s; the agnostic class does it + return gd + + #========= + def scenario_names_creator(self, num_scens,start=None): + return self.model_module.scenario_names_creator(num_scens,start) + + + #========= + def inparser_adder(self, cfg): + self.model_module.inparser_adder(cfg) + + + #========= + def kw_creator(self, cfg): + # creates keywords for scenario creator + return self.model_module.kw_creator(cfg) + + # This is not needed for PH + def sample_tree_scen_creator(self, sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + return self.model_module.sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario, **scenario_creator_kwargs) + + #============================ + def scenario_denouement(self, rank, scenario_name, scenario): + pass + # (the fct in farmer won't work because the Var names don't match) + #self.model_module.scenario_denouement(rank, scenario_name, scenario) + + + ############################################################################ + # begin callouts + # NOTE: the callouts all take the Ag object as their first argument, mainly to see cfg if needed + # the function names correspond to function names in mpisppy + + def attach_Ws_and_prox(self, Ag, sname, scenario): + # Attach W's and prox to the guest scenario. + # Use the nonant index as the index set + gs = scenario._agnostic_dict["scenario"] # guest scenario handle + nonant_idx = list(scenario._agnostic_dict["nonants"].keys()) + assert not hasattr(gs, "W") + gs.W = pyo.Param(nonant_idx, initialize=0.0, mutable=True) + assert not hasattr(gs, "W_on") + gs.W_on = pyo.Param(initialize=0, mutable=True, within=pyo.Binary) + assert not hasattr(gs, "prox_on") + gs.prox_on = pyo.Param(initialize=0, mutable=True, within=pyo.Binary) + assert not hasattr(gs, "rho") + gs.rho = pyo.Param(nonant_idx, mutable=True, default=Ag.cfg.default_rho) + + + def _disable_prox(self, Ag, scenario): + scenario._agnostic_dict["scenario"].prox_on._value = 0 + + + def _disable_W(self, Ag, scenario): + scenario._agnostic_dict["scenario"].W_on._value = 0 + + + def _reenable_prox(self, Ag, scenario): + scenario._agnostic_dict["scenario"].prox_on._value = 1 + + + def _reenable_W(self, Ag, scenario): + scenario._agnostic_dict["scenario"].W_on._value = 1 + + + def attach_PH_to_objective(self, Ag, sname, scenario, add_duals, add_prox): + # TBD: Deal with prox linearization and approximation later, + # i.e., just do the quadratic version + + ### The host has xbars and computes without involving the guest language + ### xbars = scenario._mpisppy_model.xbars + ### but instead, we are going to make guest xbars like other guests + + + gd = scenario._agnostic_dict + gs = gd["scenario"] # guest scenario handle + nonant_idx = list(gd["nonants"].keys()) + # for Pyomo, we can just ask what is the active objective function + # (from some guests, maybe we will have to put the obj function on gd + objfct = sputils.find_active_objective(gs) + ph_term = 0 + gs.xbars = pyo.Param(nonant_idx, mutable=True) + # Dual term (weights W) + if add_duals: + gs.WExpr = pyo.Expression(expr= sum(gs.W[ndn_i] * xvar for ndn_i,xvar in gd["nonants"].items())) + ph_term += gs.W_on * gs.WExpr + + # Prox term (quadratic) + if (add_prox): + prox_expr = 0. + for ndn_i, xvar in gd["nonants"].items(): + # expand (x - xbar)**2 to (x**2 - 2*xbar*x + xbar**2) + # x**2 is the only qradratic term, which might be + # dealt with differently depending on user-set options + if xvar.is_binary(): + xvarsqrd = xvar + else: + xvarsqrd = xvar**2 + prox_expr += (gs.rho[ndn_i] / 2.0) * \ + (xvarsqrd - 2.0 * gs.xbars[ndn_i] * xvar + gs.xbars[ndn_i]**2) + gs.ProxExpr = pyo.Expression(expr=prox_expr) + ph_term += gs.prox_on * gs.ProxExpr + + if gd["sense"] == pyo.minimize: + objfct.expr += ph_term + elif gd["sense"] == pyo.maximize: + objfct.expr -= ph_term + else: + raise RuntimeError(f"Unknown sense {gd['sense'] =}") + + + def solve_one(self, Ag, s, solve_keyword_args, gripe, tee=False, need_solution=True): + # This needs to attach stuff to s (see solve_one in spopt.py) + # Solve the guest language version, then copy values to the host scenario + + # This function needs to update W on the guest right before the solve + + # We need to operate on the guest scenario, not s; however, attach things to s (the host scenario) + # and copy to s. If you are working on a new guest, you should not have to edit the s side of things + + # To acommdate the solve_one call from xhat_eval.py, we need to attach the obj fct value to s + + self._copy_Ws_xbars_rho_from_host(s) + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + + # print(f" in _solve_one {global_rank =}") + if global_rank == 0: + #print(f"{gs.W.pprint() =}") + #print(f"{gs.xbars.pprint() =}") + pass + solver_name = s._solver_plugin.name + solver = pyo.SolverFactory(solver_name) + if 'persistent' in solver_name: + raise RuntimeError("Persistent solvers are not currently supported in the pyomo agnostic example.") + ###solver.set_instance(ef, symbolic_solver_labels=True) + ###solver.solve(tee=True) + else: + solver_exception = None + try: + results = solver.solve(gs, tee=tee, symbolic_solver_labels=True,load_solutions=False) + except Exception as e: + results = None + solver_exception = e + + if (results is None) or (len(results.solution) == 0) or \ + (results.solution(0).status == SolutionStatus.infeasible) or \ + (results.solver.termination_condition == TerminationCondition.infeasible) or \ + (results.solver.termination_condition == TerminationCondition.infeasibleOrUnbounded) or \ + (results.solver.termination_condition == TerminationCondition.unbounded): + + s._mpisppy_data.scenario_feasible = False + + if gripe: + print (f"Solve failed for scenario {s.name} on rank {global_rank}") + if results is not None: + print ("status=", results.solver.status) + print ("TerminationCondition=", + results.solver.termination_condition) + + if solver_exception is not None and need_solution: + raise solver_exception + + else: + s._mpisppy_data.scenario_feasible = True + if gd["sense"] == pyo.minimize: + s._mpisppy_data.outer_bound = results.Problem[0].Lower_bound + else: + s._mpisppy_data.outer_bound = results.Problem[0].Upper_bound + gs.solutions.load_from(results) + # copy the nonant x values from gs to s so mpisppy can use them in s + for ndn_i, gxvar in gd["nonants"].items(): + # courtesy check for staleness on the guest side before the copy + if not gxvar.fixed and gxvar.stale: + try: + float(pyo.value(gxvar)) + except: # noqa + raise RuntimeError( + f"Non-anticipative variable {gxvar.name} on scenario {s.name} " + "reported as stale. This usually means this variable " + "did not appear in any (active) components, and hence " + "was not communicated to the subproblem solver. ") + + s._mpisppy_data.nonant_indices[ndn_i]._value = gxvar._value + + # the next line ignore bundles (other than proper bundles) + s._mpisppy_data._obj_from_agnostic = pyo.value(sputils.get_objs(gs)[0]) + + + # local helper + def _copy_Ws_xbars_rho_from_host(self, s): + # This is an important function because it allows us to capture whatever the host did + # print(f" {s.name =}, {global_rank =}") + gd = s._agnostic_dict + gs = gd["scenario"] # guest scenario handle + for ndn_i, gxvar in gd["nonants"].items(): + assert hasattr(s, "_mpisppy_model"),\ + f"what the heck!! no _mpisppy_model {s.name =} {global_rank =}" + if hasattr(s._mpisppy_model, "W"): + gs.W[ndn_i] = pyo.value(s._mpisppy_model.W[ndn_i]) + gs.rho[ndn_i] = pyo.value(s._mpisppy_model.rho[ndn_i]) + gs.xbars[ndn_i] = pyo.value(s._mpisppy_model.xbars[ndn_i]) + else: + # presumably an xhatter + pass + + + # local helper + def _copy_nonants_from_host(self, s): + # values and fixedness; + gd = s._agnostic_dict + for ndn_i, gxvar in gd["nonants"].items(): + hostVar = s._mpisppy_data.nonant_indices[ndn_i] + guestVar = gd["nonants"][ndn_i] + if guestVar.is_fixed(): + guestVar.fixed = False + if hostVar.is_fixed(): + guestVar.fix(hostVar._value) + else: + guestVar._value = hostVar._value + + + def _restore_nonants(self, Ag, s): + # the host has already restored + self._copy_nonants_from_host(s) + + + def _restore_original_fixedness(self, Ag, s): + self._copy_nonants_from_host(s) + + + def _fix_nonants(self, Ag, s): + # We are assuming the host did the fixing + self._copy_nonants_from_host(s) + + + def _fix_root_nonants(self, Ag, s): + # We are assuming the host did the fixing + self._copy_nonants_from_host(s) + + diff --git a/mpisppy/opt/lshaped.py b/mpisppy/opt/lshaped.py index c79a8a21..faabf050 100644 --- a/mpisppy/opt/lshaped.py +++ b/mpisppy/opt/lshaped.py @@ -238,7 +238,7 @@ def _create_root_with_scenarios(self): if len(ef_scenarios) > 1: def scenario_creator_wrapper(name, **creator_options): scenario = self.scenario_creator(name, **creator_options) - if not hasattr(scenario, '_mpisppy_probability'): + if not hasattr(scenario, '_mpisppy_probability') or scenario._mpisppy_probability == "uniform": scenario._mpisppy_probability = 1./len(self.all_scenario_names) return scenario root = sputils.create_EF( diff --git a/mpisppy/phbase.py b/mpisppy/phbase.py index 27751e62..745e4bc8 100644 --- a/mpisppy/phbase.py +++ b/mpisppy/phbase.py @@ -239,6 +239,8 @@ class PHBase(mpisppy.spopt.SPOpt): Function to set rho values throughout the PH algorithm. variable_probability (callable, optional): Function to set variable specific probabilities. + cfg (config object, optional?) controls (mainly from user) + (Maybe this should move up to spbase) """ def __init__( @@ -275,6 +277,8 @@ def __init__( # Note that options can be manipulated from outside on-the-fly. # self.options (from super) will archive the original options. self.options = options + self.Ag = options.get("Ag", None) # The Agnostic Object + self.options_check() self.ph_converger = ph_converger self.rho_setter = rho_setter @@ -373,7 +377,7 @@ def convergence_diff(self): def _populate_W_cache(self, cache, padding): - """ Copy the W values for noants *for all local scenarios* + """ Copy the W values for nonants *for all local scenarios* Args: cache (np vector) to receive the W's for all local scenarios (for sending) @@ -438,6 +442,8 @@ def _use_rho_setter(self, verbose): def _disable_prox(self): for k, scenario in self.local_scenarios.items(): scenario._mpisppy_model.prox_on = 0 + if self.Ag is not None: + self.Ag.callout_agnostic({"scenario": scenario}) def _disable_W(self): @@ -446,6 +452,8 @@ def _disable_W(self): # probably not mathematically useful for scenario in self.local_scenarios.values(): scenario._mpisppy_model.W_on = 0 + if self.Ag is not None: + self.Ag.callout_agnostic({"scenario": scenario}) def disable_W_and_prox(self): @@ -456,12 +464,16 @@ def disable_W_and_prox(self): def _reenable_prox(self): for k, scenario in self.local_scenarios.items(): scenario._mpisppy_model.prox_on = 1 + if self.Ag is not None: + self.Ag.callout_agnostic({"scenario": scenario}) def _reenable_W(self): # TODO: we should eliminate this method for k, scenario in self.local_scenarios.items(): scenario._mpisppy_model.W_on = 1 + if self.Ag is not None: + self.Ag.callout_agnostic({"scenario": scenario}) def reenable_W_and_prox(self): @@ -637,6 +649,8 @@ def attach_Ws_and_prox(self): scenario._mpisppy_model.rho = pyo.Param(scenario._mpisppy_data.nonant_indices.keys(), mutable=True, default=self.options["defaultPHrho"]) + if self.Ag is not None: + self.Ag.callout_agnostic({"sname":sname, "scenario":scenario}) def attach_smoothing(self): @@ -763,6 +777,11 @@ def attach_PH_to_objective(self, add_duals, add_prox, add_smooth=0): objfct.expr += ph_term else: objfct.expr -= ph_term + + if self.Ag is not None: + self.Ag.callout_agnostic({"sname":sname, "scenario":scenario, + "add_duals": add_duals, "add_prox": add_prox}) + def PH_Prep( diff --git a/mpisppy/scenario_tree.py b/mpisppy/scenario_tree.py index febc2ffd..7898c0f1 100644 --- a/mpisppy/scenario_tree.py +++ b/mpisppy/scenario_tree.py @@ -11,6 +11,7 @@ import logging import pyomo.environ as pyo +import mpisppy.utils.sputils as sputils from pyomo.core.base.indexed_component_slice import IndexedComponent_slice logger = logging.getLogger('mpisppy.scenario_tree') @@ -87,7 +88,7 @@ def __init__(self, name, cond_prob, stage, cost_expression, self.parent_name = parent_name # None for ROOT # now make the vardata lists if self.nonant_list is not None: - self.nonant_vardata_list = build_vardatalist(self, + self.nonant_vardata_list = sputils.build_vardatalist( scen_model, self.nonant_list) else: @@ -96,7 +97,7 @@ def __init__(self, name, cond_prob, stage, cost_expression, self.nonant_vardata_list = [] if self.nonant_ef_suppl_list is not None: - self.nonant_ef_suppl_vardata_list = build_vardatalist(self, + self.nonant_ef_suppl_vardata_list = sputils.build_vardatalist( scen_model, self.nonant_ef_suppl_list) else: diff --git a/mpisppy/spopt.py b/mpisppy/spopt.py index 23d1a036..d17c8824 100644 --- a/mpisppy/spopt.py +++ b/mpisppy/spopt.py @@ -185,66 +185,76 @@ def _vb(msg): # solve_keyword_args["use_signal_handling"] = False pass - try: - results = s._solver_plugin.solve(s, - **solve_keyword_args, - load_solutions=False) - solver_exception = None - except Exception as e: - results = None - solver_exception = e - - pyomo_solve_time = time.time() - solve_start_time - if sputils.not_good_enough_results(results): - s._mpisppy_data.scenario_feasible = False - - if gripe: - name = self.__class__.__name__ - if self.spcomm: - name = self.spcomm.__class__.__name__ - print (f"[{name}] Solve failed for scenario {s.name}") - if results is not None: - print ("status=", results.solver.status) - print ("TerminationCondition=", - results.solver.termination_condition) - else: - print("no results object, so solving agin with tee=True") - solve_keyword_args["tee"] = True - results = s._solver_plugin.solve(s, - **solve_keyword_args, - load_solutions=False) - - if solver_exception is not None: - raise solver_exception - + Ag = getattr(self, "Ag", None) # agnostic + if Ag is not None: + assert not disable_pyomo_signal_handling, "Not thinking about agnostic APH yet" + kws = {"s": s, "solve_keyword_args": solve_keyword_args, "gripe": gripe, "tee": tee, "need_solution": need_solution} + Ag.callout_agnostic(kws) else: try: - if sputils.is_persistent(s._solver_plugin): - s._solver_plugin.load_vars() - else: - s.solutions.load_from(results) - except Exception as e: # catch everything - if need_solution: - raise e - if self.is_minimizing: - s._mpisppy_data.outer_bound = results.Problem[0].Lower_bound - s._mpisppy_data.inner_bound = results.Problem[0].Upper_bound - else: - s._mpisppy_data.outer_bound = results.Problem[0].Upper_bound - s._mpisppy_data.inner_bound = results.Problem[0].Lower_bound - s._mpisppy_data.scenario_feasible = True - # TBD: get this ready for IPopt (e.g., check feas_prob every time) - # propogate down - if self.bundling: # must be a bundle - for sname in s._ef_scenario_names: - self.local_scenarios[sname]._mpisppy_data.scenario_feasible\ - = s._mpisppy_data.scenario_feasible - if s._mpisppy_data.scenario_feasible: - self._check_staleness(self.local_scenarios[sname]) - else: # not a bundle - if s._mpisppy_data.scenario_feasible: - self._check_staleness(s) + results = s._solver_plugin.solve(s, + **solve_keyword_args, + load_solutions=False) + solver_exception = None + except Exception as e: + results = None + solver_exception = e + + if sputils.not_good_enough_results(results): + s._mpisppy_data.scenario_feasible = False + + if gripe: + name = self.__class__.__name__ + if self.spcomm: + name = self.spcomm.__class__.__name__ + print (f"[{name}] Solve failed for scenario {s.name}") + if results is not None: + print ("status=", results.solver.status) + print ("TerminationCondition=", + results.solver.termination_condition) + else: + print("no results object, so solving agin with tee=True") + solve_keyword_args["tee"] = True + results = s._solver_plugin.solve(s, + **solve_keyword_args, + load_solutions=False) + if solver_exception is not None: + raise solver_exception + + else: + try: + if sputils.is_persistent(s._solver_plugin): + s._solver_plugin.load_vars() + else: + s.solutions.load_from(results) + except Exception as e: # catch everything + if need_solution: + raise e + if self.is_minimizing: + s._mpisppy_data.outer_bound = results.Problem[0].Lower_bound + s._mpisppy_data.inner_bound = results.Problem[0].Upper_bound + else: + s._mpisppy_data.outer_bound = results.Problem[0].Upper_bound + s._mpisppy_data.inner_bound = results.Problem[0].Lower_bound + s._mpisppy_data.scenario_feasible = True + # TBD: get this ready for IPopt (e.g., check feas_prob every time) + # propogate down + if self.bundling: # must be a bundle + for sname in s._ef_scenario_names: + self.local_scenarios[sname]._mpisppy_data.scenario_feasible\ + = s._mpisppy_data.scenario_feasible + if s._mpisppy_data.scenario_feasible: + self._check_staleness(self.local_scenarios[sname]) + else: # not a bundle + if s._mpisppy_data.scenario_feasible: + self._check_staleness(s) + + # end of Agnostic bypass + + # Time capture moved down August 2023 + pyomo_solve_time = time.time() - solve_start_time + if self.extensions is not None: results = self.extobject.post_solve(s, results) @@ -364,8 +374,13 @@ def Eobjective(self, verbose=False): """ local_Eobjs = [] for k,s in self.local_scenarios.items(): - objfct = self.saved_objectives[k] - local_Eobjs.append(s._mpisppy_probability * pyo.value(objfct)) + objfct = self.saved_objectives[k] # if bundling? + Ag = getattr(self, "Ag", None) + if Ag is None: + local_Eobjs.append(s._mpisppy_probability * pyo.value(objfct)) + else: + # Agnostic will have attached the objective (and doesn't bundle as of Aug 2023) + local_Eobjs.append(s._mpisppy_probability * s._mpisppy_data._obj_from_agnostic) if verbose: print ("caller", inspect.stack()[1][3]) print ("E_Obj Scenario {}, prob={}, Obj={}, ObjExpr={}"\ @@ -586,6 +601,9 @@ def _restore_original_fixedness(self): for k,s in self.local_scenarios.items(): for ci, _ in enumerate(s._mpisppy_data.nonant_indices): s._mpisppy_data.fixedness_cache[ci] = s._mpisppy_data.original_fixedness[ci] + Ag = getattr(self, "Ag", None) + if Ag is not None: + Ag.callout_agnostic({"s": s}) self._restore_nonants() @@ -628,6 +646,10 @@ def _fix_nonants(self, cache): if persistent_solver is not None: persistent_solver.update_var(this_vardata) + Ag = getattr(self, "Ag", None) + if Ag is not None: + Ag.callout_agnostic({"s": s}) + def _fix_root_nonants(self,root_cache): """ Fix the 1st stage Vars subject to non-anticipativity at given values. Loop over the scenarios to restore, but loop over subproblems @@ -673,7 +695,10 @@ def _fix_root_nonants(self,root_cache): this_vardata.fix() if persistent_solver is not None: persistent_solver.update_var(this_vardata) - + + Ag = getattr(self, "Ag", None) + if Ag is not None: + Ag.callout_agnostic({"s": s}) def _restore_nonants(self, update_persistent=True): @@ -701,6 +726,10 @@ def _restore_nonants(self, update_persistent=True): if persistent_solver is not None: persistent_solver.update_var(vardata) + + Ag = getattr(self, "Ag", None) + if Ag is not None: + Ag.callout_agnostic({"s": s}) def _save_nonants(self): @@ -735,9 +764,6 @@ def _save_original_nonants(self): the variable was fixed is stored in `_PySP_original_fixedness`. """ for k,s in self.local_scenarios.items(): - if hasattr(s,"_PySP_original_fixedness"): - print ("ERROR: Attempt to replace original nonants") - raise if not hasattr(s._mpisppy_data,"nonant_cache"): # uses nonant cache to signal other things have not # been created diff --git a/mpisppy/tests/farmer.mod b/mpisppy/tests/farmer.mod new file mode 100644 index 00000000..1e0c9b3b --- /dev/null +++ b/mpisppy/tests/farmer.mod @@ -0,0 +1,110 @@ +# The farmer's problem in AMPL +# +# Reference: +# John R. Birge and Francois Louveaux. Introduction to Stochastic Programming. +# +# AMPL coding by Victor Zverovich; ## modifed by dlw; now *minization* + +##function expectation; +##function random; + +##suffix stage IN; + +set Crops; + +##set Scen; +##param P{Scen}; # probabilities + +param TotalArea; # acre +param PlantingCost{Crops}; # $/acre +param SellingPrice{Crops}; # $/T +param ExcessSellingPrice; # $/T +param PurchasePrice{Crops}; # $/T +param MinRequirement{Crops}; # T +param BeetsQuota; # T + +# Area in acres devoted to crop c. +var area{c in Crops} >= 0; + +# Tons of crop c sold (at favourable price) under scenario s. +var sell{c in Crops} >= 0, suffix stage 2; + +# Tons of sugar beets sold in excess of the quota under scenario s. +var sell_excess >= 0, suffix stage 2; + +# Tons of crop c bought under scenario s +var buy{c in Crops} >= 0, suffix stage 2; + +# The random variable (parameter) representing the yield of crop c. +##var RandomYield{c in Crops}; +param RandomYield{c in Crops}; + +# Realizations of the yield of crop c. +##param Yield{c in Crops, s in Scen}; # T/acre + +##maximize profit: +## expectation( +## ExcessSellingPrice * sell_excess + +## sum{c in Crops} (SellingPrice[c] * sell[c] - +## PurchasePrice[c] * buy[c])) - +## sum{c in Crops} PlantingCost[c] * area[c]; + +minimize minus_profit: + - ExcessSellingPrice * sell_excess - + sum{c in Crops} (SellingPrice[c] * sell[c] - + PurchasePrice[c] * buy[c]) + + sum{c in Crops} (PlantingCost[c] * area[c]); + +s.t. totalArea: sum {c in Crops} area[c] <= TotalArea; + +s.t. requirement{c in Crops}: + RandomYield[c] * area[c] - sell[c] + buy[c] >= MinRequirement[c]; + +s.t. quota: sell['beets'] <= BeetsQuota; + +s.t. sellBeets: + sell['beets'] + sell_excess <= RandomYield['beets'] * area['beets']; + +##yield: random({c in Crops} (RandomYield[c], {s in Scen} Yield[c, s])); + +data; + +set Crops := wheat corn beets; +#set Scen := below average above; + +param TotalArea := 500; + +##param Yield: +## below average above := +## wheat 2.0 2.5 3.0 +## corn 2.4 3.0 3.6 +## beets 16.0 20.0 24.0; + +param RandomYield := + wheat 2.5 + corn 3.0 + beets 20.0; + +param PlantingCost := + wheat 150 + corn 230 + beets 260; + +param SellingPrice := + wheat 170 + corn 150 + beets 36; + +param ExcessSellingPrice := 10; + +param PurchasePrice := + wheat 238 + corn 210 + beets 100; + +param MinRequirement := + wheat 200 + corn 240 + beets 0; + +param BeetsQuota := 6000; diff --git a/mpisppy/tests/test_agnostic.py b/mpisppy/tests/test_agnostic.py new file mode 100644 index 00000000..c1aea6e9 --- /dev/null +++ b/mpisppy/tests/test_agnostic.py @@ -0,0 +1,260 @@ +############################################################################### +# 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. +############################################################################### +# ruff: noqa: F841 + +import sys +import unittest +import mpisppy.opt.ph +from mpisppy.tests.utils import get_solver +import mpisppy.utils.config as config +import mpisppy.agnostic.agnostic as agnostic +import mpisppy.agnostic.agnostic_cylinders as agnostic_cylinders +import mpisppy.utils.sputils as sputils + +sys.path.insert(0, "../../examples/farmer/agnostic") +import farmer_pyomo_agnostic +try: + import mpisppy.agnostic.gams_guest + have_GAMS = True +except ModuleNotFoundError: + have_GAMS = False +try: + import farmer_ampl_agnostic + have_AMPL = True +except ModuleNotFoundError: + have_AMPL = False +try: + import farmer_gurobipy_agnostic + have_gurobipy = True +except ModuleNotFoundError: + have_gurobipy = False + +__version__ = 0.2 + +solver_available, solver_name, persistent_available, persistent_solver_name = ( + get_solver(persistent_OK=False) +) + +# NOTE Gurobi is hardwired for the AMPL and GAMS tests, so don't install it on github +# (and, if you have gurobi installed the ampl test will fail) + + +def _farmer_cfg(): + cfg = config.Config() + cfg.popular_args() + cfg.ph_args() + cfg.default_rho = 1 + farmer_pyomo_agnostic.inparser_adder(cfg) + return cfg + + +def _get_ph_base_options(): + Baseoptions = {} + Baseoptions["asynchronousPH"] = False + Baseoptions["solver_name"] = solver_name + Baseoptions["PHIterLimit"] = 3 + Baseoptions["defaultPHrho"] = 1 + Baseoptions["convthresh"] = 0.001 + Baseoptions["subsolvedirectives"] = None + Baseoptions["verbose"] = False + Baseoptions["display_timing"] = False + Baseoptions["display_progress"] = False + if "cplex" in solver_name: + Baseoptions["iter0_solver_options"] = {"mip_tolerances_mipgap": 0.001} + Baseoptions["iterk_solver_options"] = {"mip_tolerances_mipgap": 0.00001} + else: + Baseoptions["iter0_solver_options"] = {"mipgap": 0.001} + Baseoptions["iterk_solver_options"] = {"mipgap": 0.00001} + + Baseoptions["display_progress"] = False + + return Baseoptions + + +# ***************************************************************************** + +class Test_Agnostic_pyomo(unittest.TestCase): + + def test_agnostic_pyomo_constructor(self): + cfg = _farmer_cfg() + agnostic.Agnostic(farmer_pyomo_agnostic, cfg) + + + def test_agnostic_pyomo_scenario_creator(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_pyomo_agnostic, cfg) + s0 = Ag.scenario_creator("scen0") + s2 = Ag.scenario_creator("scen2") + + + def test_agnostic_pyomo_PH_constructor(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_pyomo_agnostic, cfg) + s1 = Ag.scenario_creator("scen1") # average case + phoptions = _get_ph_base_options() + ph = mpisppy.opt.ph.PH( + phoptions, + farmer_pyomo_agnostic.scenario_names_creator(num_scens=3), + Ag.scenario_creator, + farmer_pyomo_agnostic.scenario_denouement, + scenario_creator_kwargs=None, # agnostic.py takes care of this + extensions=None + ) + + @unittest.skipIf(not solver_available, + "no solver is available") + def test_agnostic_pyomo_PH(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_pyomo_agnostic, cfg) + s1 = Ag.scenario_creator("scen1") # average case + phoptions = _get_ph_base_options() + phoptions["Ag"] = Ag # this is critical + scennames = farmer_pyomo_agnostic.scenario_names_creator(num_scens=3) + ph = mpisppy.opt.ph.PH( + phoptions, + scennames, + Ag.scenario_creator, + farmer_pyomo_agnostic.scenario_denouement, + scenario_creator_kwargs=None, # agnostic.py takes care of this + extensions=None + ) + conv, obj, tbound = ph.ph_main() + self.assertAlmostEqual(-115405.5555, tbound, places=1) + self.assertAlmostEqual(-110433.4007, obj, places=1) + + +@unittest.skipIf(not have_AMPL, "skipping AMPL") +class Test_Agnostic_AMPL(unittest.TestCase): + def test_agnostic_AMPL_constructor(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_ampl_agnostic, cfg) + + def test_agnostic_AMPL_scenario_creator(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_ampl_agnostic, cfg) + s0 = Ag.scenario_creator("scen0") + s2 = Ag.scenario_creator("scen2") + + @unittest.skipIf(not solver_available, "no solver is available") + def test_agnostic_ampl_PH(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_ampl_agnostic, cfg) + s1 = Ag.scenario_creator("scen1") # average case + phoptions = _get_ph_base_options() + phoptions["Ag"] = Ag # this is critical + phoptions["solver_name"] = "gurobi" # need an ampl solver + scennames = farmer_ampl_agnostic.scenario_names_creator(num_scens=3) + ph = mpisppy.opt.ph.PH( + phoptions, + scennames, + Ag.scenario_creator, + farmer_ampl_agnostic.scenario_denouement, + scenario_creator_kwargs=None, # agnostic.py takes care of this + extensions=None, + ) + conv, obj, tbound = ph.ph_main() + print(f"{obj =}, {tbound}") + print(f"{solver_name=}") + print(f"{tbound=}") + print(f"{obj=}") + message = """ NOTE if you are getting zeros it is because Gurobi is + hardwired for AMPL tests, so don't install it on github (and, if you have + gurobi generally installed on your machine, then the ampl + test will fail on your machine). """ + self.assertAlmostEqual(-115405.5555, tbound, 2, message) + self.assertAlmostEqual(-110433.4007, obj, 2, message) + + def test_agnostic_cylinders_ampl(self): + # just make sure PH runs + print("test_agnostic_cylinders_ampl") + model_fname = "mpisppy.agnostic.examples.farmer_ampl_model" + module = sputils.module_name_to_module(model_fname) + cfg = agnostic_cylinders._setup_args(module) + cfg.module_name = model_fname + cfg.default_rho = 1 + cfg.num_scens = 3 + cfg.solver_name= "gurobi" + cfg.guest_language = "AMPL" + cfg.max_iterations = 5 + cfg.ampl_model_file = "../agnostic/examples/farmer.mod" + agnostic_cylinders.main(model_fname, module, cfg) + + +@unittest.skipIf(not have_GAMS, "skipping GAMS") +class Test_Agnostic_GAMS(unittest.TestCase): + def test_agnostic_cylinders_gams(self): + # just make sure PH runs + print("test_agnostic_cylinders_gams") + model_fname = "mpisppy.agnostic.examples.farmer_gams_model" + module = sputils.module_name_to_module(model_fname) + cfg = agnostic_cylinders._setup_args(module) + cfg.module_name = model_fname + cfg.default_rho = 1 + cfg.num_scens = 3 + cfg.solver_name= "gurobi" + cfg.guest_language = "GAMS" + cfg.max_iterations = 5 + cfg.gams_model_file = "../agnostic/examples/farmer_average.gms" + agnostic_cylinders.main(model_fname, module, cfg) + + +@unittest.skipIf(not have_gurobipy, "skipping gurobipy") +class Test_Agnostic_gurobipy(unittest.TestCase): + def test_agnostic_gurobipy_constructor(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_gurobipy_agnostic, cfg) + + def test_agnostic_gurobipy_scenario_creator(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_gurobipy_agnostic, cfg) + s0 = Ag.scenario_creator("scen0") + s2 = Ag.scenario_creator("scen2") + + def test_agnostic_gurobipy_PH_constructor(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_gurobipy_agnostic, cfg) + s1 = Ag.scenario_creator("scen1") # average case + phoptions = _get_ph_base_options() + ph = mpisppy.opt.ph.PH( + phoptions, + farmer_gurobipy_agnostic.scenario_names_creator(num_scens=3), + Ag.scenario_creator, + farmer_gurobipy_agnostic.scenario_denouement, + scenario_creator_kwargs=None, # agnostic.py takes care of this + extensions=None, + ) + + + @unittest.skipIf(not solver_available, "no solver is available") + # Test Case is failing + def test_agnostic_gurobipy_PH(self): + cfg = _farmer_cfg() + Ag = agnostic.Agnostic(farmer_gurobipy_agnostic, cfg) + s1 = Ag.scenario_creator("scen1") # average case + phoptions = _get_ph_base_options() + phoptions["Ag"] = Ag # this is critical + scennames = farmer_gurobipy_agnostic.scenario_names_creator(num_scens=3) + ph = mpisppy.opt.ph.PH( + phoptions, + scennames, + Ag.scenario_creator, + farmer_gurobipy_agnostic.scenario_denouement, + scenario_creator_kwargs=None, # agnostic.py takes care of this + extensions=None, + ) + conv, obj, tbound = ph.ph_main() + print(f"{solver_name=}") + print(f"{tbound=}") + print(f"{obj=}") + self.assertAlmostEqual(-110433.4007, obj, places=1) + self.assertAlmostEqual(-115405.5555, tbound, places=1) + + +if __name__ == "__main__": + unittest.main() diff --git a/mpisppy/tests/test_headers.py b/mpisppy/tests/test_headers.py index 36c1302a..327f26db 100644 --- a/mpisppy/tests/test_headers.py +++ b/mpisppy/tests/test_headers.py @@ -31,4 +31,5 @@ def test_headers(): if p.stat().st_size == 0: continue nonempty_missing_header.append(p) + print(f"{nonempty_missing_header=}") assert not nonempty_missing_header diff --git a/mpisppy/tests/test_stoch_admmWrapper.py b/mpisppy/tests/test_stoch_admmWrapper.py index 7f9e7618..ccc0ed60 100644 --- a/mpisppy/tests/test_stoch_admmWrapper.py +++ b/mpisppy/tests/test_stoch_admmWrapper.py @@ -119,11 +119,12 @@ def test_values(self): if result.stderr: print("Error output:") print(result.stderr) - + raise RuntimeError("Error encountered as shown above.") # Check the standard output if result.stdout: result_by_line = result.stdout.strip().split('\n') - + else: + raise RuntimeError(f"No results in stdout for {command=}.") target_line = "Iter. Best Bound Best Incumbent Rel. Gap Abs. Gap" precedent_line_target = False i = 0 diff --git a/mpisppy/tests/utils.py b/mpisppy/tests/utils.py index f3b3f102..726f6aae 100644 --- a/mpisppy/tests/utils.py +++ b/mpisppy/tests/utils.py @@ -11,8 +11,10 @@ import pyomo.environ as pyo from math import log10, floor -def get_solver(): - solvers = [n+e for e in ('_persistent', '') for n in ("cplex","gurobi","xpress")] +def get_solver(persistent_OK=True): + solvers = ["cplex","gurobi","xpress"] + if persistent_OK: + solvers = [n+e for e in ('_persistent', '') for n in solvers] for solver_name in solvers: try: diff --git a/mpisppy/utils/sputils.py b/mpisppy/utils/sputils.py index 5118bf65..40558d6f 100644 --- a/mpisppy/utils/sputils.py +++ b/mpisppy/utils/sputils.py @@ -14,6 +14,8 @@ import os import re import numpy as np +import inspect +import importlib import mpisppy.scenario_tree as scenario_tree from pyomo.core import Objective from pyomo.repn import generate_standard_repn @@ -21,11 +23,47 @@ from mpisppy import MPI, haveMPI from pyomo.core.expr.numeric_expr import LinearExpression from pyomo.opt import SolutionStatus, TerminationCondition +from pyomo.core.base.indexed_component_slice import IndexedComponent_slice from mpisppy import tt_timer global_rank = MPI.COMM_WORLD.Get_rank() + +def build_vardatalist(model, varlist=None): + """ + Convert a list of pyomo variables to a list of SimpleVar and _GeneralVarData. If varlist is none, builds a + list of all variables in the model. Written by CD Laird + + Parameters + ---------- + model: ConcreteModel + varlist: None or list of pyo.Var + """ + vardatalist = None + + # if the varlist is None, then assume we want all the active variables + if varlist is None: + raise RuntimeError("varlist is None in scenario_tree.build_vardatalist") + vardatalist = [v for v in model.component_data_objects(pyo.Var, active=True, sort=True)] + elif isinstance(varlist, (pyo.Var, IndexedComponent_slice)): + # user provided a variable, not a list of variables. Let's work with it anyway + varlist = [varlist] + + if vardatalist is None: + # expand any indexed components in the list to their + # component data objects + vardatalist = list() + for v in varlist: + if isinstance(v, IndexedComponent_slice): + vardatalist.extend(v.__iter__()) + elif v.is_indexed(): + vardatalist.extend((v[i] for i in sorted(v.keys()))) + else: + vardatalist.append(v) + return vardatalist + + def not_good_enough_results(results): return (results is None) or (len(results.solution) == 0) or \ (results.solution(0).status == SolutionStatus.infeasible) or \ @@ -106,8 +144,9 @@ def get_objs(scenario_instance, allow_none=False): raise RuntimeError(f"Scenario {scenario_instance.name} has no active " "objective functions.") if (len(scenario_objs) > 1): - print("WARNING: Scenario", scenario_instance.name, "has multiple active " - "objectives. Selecting the first objective.") + print(f"WARNING: Scenario {scenario_instance.name} has multiple active " + "objectives, returning a list.") + return scenario_objs @@ -1037,7 +1076,15 @@ def number_of_nodes(branching_factors): #How many nodes does a tree with a given branching_factors have ? last_node_stage_num = [i-1 for i in branching_factors] return node_idx(last_node_stage_num, branching_factors) - + + +def module_name_to_module(module_name): + if inspect.ismodule(module_name): + module = module_name + else: + module = importlib.import_module(module_name) + return module + if __name__ == "__main__": branching_factors = [2,2,2,3] @@ -1057,3 +1104,4 @@ def number_of_nodes(branching_factors): print(ndn, v) print(f"slices: {slices}") check4losses(numscens, branching_factors, sntr, slices, ranks_per_scenario) + diff --git a/mpisppy/utils/xhat_eval.py b/mpisppy/utils/xhat_eval.py index fd93dd72..cfd378e2 100644 --- a/mpisppy/utils/xhat_eval.py +++ b/mpisppy/utils/xhat_eval.py @@ -60,6 +60,7 @@ def __init__( self.verbose = self.options['verbose'] self._subproblems_solvers_created = False + self.Ag = options.get("Ag", None) def _lazy_create_solvers(self): @@ -88,7 +89,6 @@ def solve_one(self, solver_options, k, s, disable_pyomo_signal_handling=disable_pyomo_signal_handling, update_objective=update_objective) - if compute_val_at_nonant: objfct = self.saved_objectives[k] if self.verbose: @@ -96,6 +96,7 @@ def solve_one(self, solver_options, k, s, print ("E_Obj Scenario {}, prob={}, Obj={}, ObjExpr={}"\ .format(k, s._mpisppy_probability, pyo.value(objfct), objfct.expr)) self.objs_dict[k] = pyo.value(objfct) + return(pyomo_solve_time) @@ -279,7 +280,7 @@ def evaluate(self, nonant_cache, fct=None): ) Eobj = self.Eobjective(self.verbose,fct=fct) - + return Eobj