| I l@ve RuBoard |
Using the Examples and DemosI want to briefly describe how to use the book's examples here. In general, though, please see the following text files in the examples distribution directory for more details:
Of these, the README-PP2E.txt file is the most informative, and the PP2E\Config directory contains all configuration file examples. I give an overview here, but the files listed give a complete description. The Short StoryIf you want to see some Python examples right away, do this:
The Launch_* scripts start Python programs portably[4] and require only that Python be installed -- you don't need to set environment variables first or tweak the included PP2E\Config setup files to run them. LaunchBrowser will work if it can find a web browser on your machine, even if you don't have an Internet link (though some Internet examples won't work completely without a live link).
If installing Python isn't an option, you can still run a few Python web demos by visiting http://starship.python.net/~lutz/PyInternetDemos.html with your browser. Because these examples execute scripts on a server, they tend to work best when run live from this site, rather than from the book's CD. The DetailsTo help organize the new examples, I've provided a demo launcher program, PyDemos.pyw, in the top-level PP2E directory of the examples distribution. Figure P-1 shows PyDemos in action on Windows after pressing a few buttons. The launcher bar appears on the left of the screen; with it, you can run most of the major graphical examples in the book with a mouse click. The demo launcher bar can also be used to start major Internet book examples if a browser can be located on your machine (see the following launcher description). Besides launching demos, the PyDemos source code provides pointers to major examples in the distribution; see its source code for details. You'll also find Linux automated build scripts for the Python/C integration examples in the top-level examples directory, which serve as indexes to major C examples. Figure P-1. The PyDemos launcher with pop-ups and demos (Guido's photo reprinted with permission from Dr. Dobb's Journal)I've also included a top-level program called PyGadgets.py, and its relative PyGadgets_bar.pyw, to launch some of the more useful GUI book examples for real use instead of demonstration (mostly, the programs I use often; configure as desired). Figure P-2 shows what PyGadgets_bar looks like on Windows, along with a few of the utilities that its buttons can launch. All of the programs are presented in this book and included in the examples distribution. Most gadgets require a Python with Tkinter support, but that is the default configuration for the standard Windows port on the book's CD (see http://examples.oreilly.com/python2). Figure P-2. The PyGadgets utilities launcher bar, with gadgetsTo run the files listed in the prior paragraph directly, you'll need to set up your Python module search path (see the top-level PP2E/Config/setup* files for hints). But if you want to run a collection of Python demos from the book and don't want to bother with setting up your environment first, simply run the self-launching utility scripts in the PP2E directory instead: Launch_PyDemos.pyw, Launch_PyGadgets.py, and Launch_PyGadgets_bar.pyw. These Python-coded launcher scripts assume Python has already been installed, but will automatically find your Python executable and the book examples distribution, and set up your Python module and system search paths as needed to run the demos. You can probably run these launch scripts by simply clicking on their names in a file explorer, and you should also be able to run them directly from the book's CD-ROM (see http://examples.oreilly.com/python2). See the comments at the top of Launcher.py for more details (or read about these scripts in Chapter 4 ). Many of the browser-based Internet examples from the book can also be found online at http://starship.python.net/~lutz/PyInternetDemos.html, where you can test-drive a few Python Internet scripting examples. Because these examples run in your web browser, they can be tested even if you haven't installed Python (or Python's Tk support) on your machine. The PyDemos program also attempts to launch a web browser on the major example web pages by starting the LaunchBrowser.py script in the examples root directory. That script tries to find a usable browser on your machine, with generally good results; see the script for more details if it fails. Provided LaunchBrowser can find a browser on your machine, some demo buttons will pop up web pages automatically, whether you have a live Internet connection or not (if not, you'll see local files in your browser). Figure P-3 shows what the PyInternetDemos page looks like under Internet Explorer on Windows. Figure P-3. The PyInternetDemos web pageOf special interest, the getfile.html link on this page allows you to view the source code of any other file on the book's site -- HTML code, Python CGI scripts, and so on; see Chapter 12 for details. To summarize, here is what you'll find in the top-level PP2E directory of the book's examples distribution:
There are also subdirectories for examples from each major topic area of the book. In addition, the top-level PP2E\PyTools directory contains Python-coded command-line utilities for converting line-feeds in all example text files to DOS or Unix format (useful if they look odd in your text editor), making all example files writable (useful if you drag-and-drop off the book's CD (see http://examples.oreilly.com/python2), deleting old .pyc bytecode files in the tree, and more. Again, see the distribution's README-PP2E.txt file for more details on all example issues. Where It's AtThe book examples distribution can be found on the CD-ROM that accompanies this book. See the CD's top-level README file for usage details, or browse the CD's example root directory in your favorite file explorer for a quick tour. In addition to the book examples, the CD also contains various Python-related packages, including a full Windows self-installer program with Python and Tk support (double-click and say "yes" at all prompts to install), the full Python source code distribution (unpack and compile on your machine), and Python's standard documentation set in HTML form (click to view in your web browser). Extra open source packages such as the latest releases (at the time of publication) of the SWIG code generator and JPython are also included, but you can always find up-to-date releases of Python and other packages at Python's web site, http://www.python.org. |
| I l@ve RuBoard |


This year's episodes
- : An update
- : Deconstructing Formula One
- : Pioneer 10: The Lost Traveller
- : Bizarre history of birth control
- : Beer Talk: Why learn history?
- : History of Timekeeping
- : Coffee rules the world
- : Serial Killer Files Part 4: DC Sniper
- : Why do the continents move
- : Beer Talk: Facebook, Google and online privacy
- : The Race to the South Pole
- : Blockchain Primer
- : The Match Fixing Scandal of 2000
- : Life as a comic with Rupen Paul
- : Submarines
- : Restaurant Roundup - 02
- : Serial Killer Files Part 3: The Green River Killer
- : To Moon and Back - Part 5
- : History of Motion Pictures: Part 1
- : Evolution of Languages: Dravidian - Part 1
- : History of the Toothbrush
2017
- : Serial Killer Files Part 2: Ed Gein - The Real Psycho
- : The Genographic Project (Human Migration)
- : What is the Law of Attraction?
- : To Moon and Back - Part 4
- : Freelancing and Startups
- : The War Stories: World War I & II
- : Serial Killer Files Part 1: The Night Stalker
- : Restaurant Roundup - 01
- : Tools for writing
- : Our musical evolution
- : The Neanderthals - Where did they go?
- : Podcasting Infrastructure, Hardware and Software
- : Intro to podcasting
- : Mars One - too ambitious?
- : One Mr. Sherlock Holmes
- : To Moon and Back - Part 3
- : Mangalyaan: To Mars on budget
- : 3D Printing for Practical Purposes
- : AI - evolution or revolution?
- : Mass Hysteria feat. Monkey Man
- : To Moon and Back - Part 2
- : To Moon and Back - Part 1
- : Phantom Time Hypothesis
- : Losing weight
- : Social media gives you wings
- : Does formal education do justice?
- : Welcome to the show


SERVICES
COMPANY
PORTFOLIO
CAREERS
BLOG
CONTACT
Your Digital
Workshop
Building Business Solutions & Bright Ideas


Adobe Analytics and Launch 연습을 위한 페이지입니다.
YHK | Data Consulting
Our templates are updated regularly so they don't break.
You can use this theme as is, or you can make changes!
We update dependencies to keep things fresh.
You have to make your websites with love these days!
Ready to start your next project with us? That's great! Give us a call or send us an email and we will get back to you as soon as possible!
123-456-6789


“We chose Kubernetes because of its flexibility, ease of management and the way it improves our engineers’ productivity.”
Read about Pearson“With Kubernetes, we’re simplifying our environment and making it easier for developers to build the tools that make wikis run better.”
Read about WikimediaInside eBay’s shift to Kubernetes and containers atop OpenStack
Read about eBayKubernetes at Box: Microservices at Maximum Velocity
Read about Box


--- title: 'Using Pyomo from R through the magic of Reticulate' author: 'Notesofdabbler' date: '2020-07-01' slug: 'rpyomo' comments: true categories: ["R"] tags: ["R"] ---
Pyomo is a python based open-source package for modeling optimization problems. It makes it easy to represent optimization problems and can send it to different solvers (both open-source and commercial) to solve the problem and return the results in python. The advantage of pyomo compared to commercial software such as GAMS and AMPL is the ability to code using standard python syntax (with some modifications for pyomo constructs). Another open source package for modeling optimization problems is JuMP in Julia language.
My goal in this blog is to see how far I can get in terms of using Pyomo from R using the reticulate package. The simplest option would be to develop the model in pyomo and call it from R using reticulate. However, it still requires writing the pyomo model in python. I want to use reticulate to write the pyomo model using R. In this blog post, I describe two examples in detail where I developed the pyomo model in R and discuss my learnings. I first discuss set-up in terms of packages needed followed by discussion of the two examples.
# load libraries
library(reticulate)
library(dplyr)
library(tidyr)
library(ggplot2)
For using pyomo package through reticulate, it is necessary to have python and pyomo already installed. We can import pyomo as follows:
# Import pyomo
pyo = import("pyomo.environ", convert = FALSE)
bi = import_builtins()
When trying different examples, I saw some issues when not using the convert = FALSE
option and so started using this option when importing pyomo.
I also ran into some signal handling errors when using pyomo solvers. I added the following two commands to turn off signal handling (based on this thread). It seems to have worked but I don’t know yet if it has other side effects.
pyulib = import("pyutilib.subprocess.GlobalData")
pyulib$DEFINE_SIGNAL_HANDLERS_DEFAULT = FALSE
The first example is a unconstrained nonlinear minimization of rosenbrock function \[ \min \;\; (x - 1)^2 + 100(y - x^2)^2 \] Since arithmetic operators in pyomo are overloaded to generate pyomo expressions, I wrote functions for arithmetic operators based on this discussion thread and sourced them here.
source("operators.R")
The pyomo model using python for this problem is in this location. Reticulate makes it easy to write code close to what it is in python (with certain modifications)
# Create model object
M = pyo$ConcreteModel()
# Define and initialize the variables*
M$x = pyo$Var(initialize = 1.5)
M$y = pyo$Var(initialize = 1.5)
# Define the objective function
M$obj = pyo$Objective(expr = (1 - M$x)**2 + 100 * (M$y - M$x**2)**2)
We can check the model that has been created with the pprint
command. While the R command M$pprint()
works and shows output in the console, I couldn’t get it to show in the R markdown output. So I am using the python chunk to show the output.
r.M.pprint()
## 2 Var Declarations
## x : Size=1, Index=None
## Key : Lower : Value : Upper : Fixed : Stale : Domain
## None : None : 1.5 : None : False : False : Reals
## y : Size=1, Index=None
## Key : Lower : Value : Upper : Fixed : Stale : Domain
## None : None : 1.5 : None : False : False : Reals
##
## 1 Objective Declarations
## obj : Size=1, Index=None, Active=True
## Key : Active : Sense : Expression
## None : True : minimize : (1.0 - x)**2.0 + 100.0*(y - x**2.0)**2.0
##
## 3 Declarations: x y obj
Pyomo can interface with different open-source and commercial solvers depending on the type of optimization problem. Since this is a nonlinear problem, we use the open-source solver IPOPT to solve the problem. Usually solvers have to be separately installed to be used with pyomo.
solver = pyo$SolverFactory("ipopt")
res = solver$solve(M, logfile = 'tmp.log')
Since the console output from solver doesn’t show in R markdown output, I wrote it to a temporary file and read it back again.
x = paste(readLines("tmp.log", encoding = "UTF-8"))
writeLines(x[2:63]) # lines 1 and 64 cause encoding issues in blog rendering
##
## Ipopt 3.12.12:
##
## ******************************************************************************
## This program contains Ipopt, a library for large-scale nonlinear optimization.
## Ipopt is released as open source code under the Eclipse Public License (EPL).
## For more information visit http://projects.coin-or.org/Ipopt
## ******************************************************************************
##
## This is Ipopt version 3.12.12, running with linear solver mumps.
## NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
##
## Number of nonzeros in equality constraint Jacobian...: 0
## Number of nonzeros in inequality constraint Jacobian.: 0
## Number of nonzeros in Lagrangian Hessian.............: 3
##
## Total number of variables............................: 2
## variables with only lower bounds: 0
## variables with lower and upper bounds: 0
## variables with only upper bounds: 0
## Total number of equality constraints.................: 0
## Total number of inequality constraints...............: 0
## inequality constraints with only lower bounds: 0
## inequality constraints with lower and upper bounds: 0
## inequality constraints with only upper bounds: 0
##
## iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
## 0 5.6500000e+01 0.00e+00 1.00e+02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
## 1 2.4669972e-01 0.00e+00 2.22e-01 -1.0 7.40e-01 - 1.00e+00 1.00e+00f 1
## 2 1.6256267e-01 0.00e+00 2.04e+00 -1.7 1.48e+00 - 1.00e+00 2.50e-01f 3
## 3 8.6119444e-02 0.00e+00 1.08e+00 -1.7 2.36e-01 - 1.00e+00 1.00e+00f 1
## 4 4.3223836e-02 0.00e+00 1.23e+00 -1.7 2.61e-01 - 1.00e+00 1.00e+00f 1
## 5 1.5610508e-02 0.00e+00 3.54e-01 -1.7 1.18e-01 - 1.00e+00 1.00e+00f 1
## 6 5.3544798e-03 0.00e+00 5.51e-01 -1.7 1.67e-01 - 1.00e+00 1.00e+00f 1
## 7 6.1281576e-04 0.00e+00 5.19e-02 -1.7 3.87e-02 - 1.00e+00 1.00e+00f 1
## 8 2.8893076e-05 0.00e+00 4.52e-02 -2.5 4.53e-02 - 1.00e+00 1.00e+00f 1
## 9 3.4591761e-08 0.00e+00 3.80e-04 -2.5 3.18e-03 - 1.00e+00 1.00e+00f 1
## iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
## 10 1.2680803e-13 0.00e+00 3.02e-06 -5.7 3.62e-04 - 1.00e+00 1.00e+00f 1
## 11 7.0136460e-25 0.00e+00 1.72e-12 -8.6 2.13e-07 - 1.00e+00 1.00e+00f 1
##
## Number of Iterations....: 11
##
## (scaled) (unscaled)
## Objective...............: 1.5551321399859192e-25 7.0136459513364959e-25
## Dual infeasibility......: 1.7239720368203862e-12 7.7751138860599418e-12
## Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
## Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
## Overall NLP error.......: 1.7239720368203862e-12 7.7751138860599418e-12
##
##
## Number of objective function evaluations = 18
## Number of objective gradient evaluations = 12
## Number of equality constraint evaluations = 0
## Number of inequality constraint evaluations = 0
## Number of equality constraint Jacobian evaluations = 0
## Number of inequality constraint Jacobian evaluations = 0
## Number of Lagrangian Hessian evaluations = 11
## Total CPU secs in IPOPT (w/o function evaluations) = 0.037
## Total CPU secs in NLP function evaluations = 0.000
##
## EXIT: Optimal Solution Found.
The solution of the model can be access using the defined variables of the model
# Display model solution
M$x(); M$y()
## 1.0000000000008233
## 1.0000000000016314
Given a set of markets with demands and plants that can supply the demands, the goal in the transportation problem is to determine minimum cost shipping scenario that meets market demands while satisfying supply capacity of plants. The python pyomo version of this example is in this location. The decision variables and parameters of the problem are:
The mathematical formulation of the problem is: \[ \mbox{(Minimize Cost)}\;\;\min \;\; \sum_{i=1}^P\sum_{j=1}^Mc_{ij}x_{ij} \\ \mbox{(Supply from a plant below capacity)} \;\;\sum_{j=1}^Mx_{ij} \leq Cap_i, \;\; i=1,2,\ldots,P \\ \mbox{(Supply to a market meets demand)}\;\;\sum_{i=1}^Px_{ij} \geq Dem_j, \;\; j=1,2,\ldots,M \\ \mbox{(shipments are non-negative)}\;\;x_{ij}\geq 0, \;\; i=1,2,\ldots,P, \;\; j=1,2,\ldots,M \]
Cost is a dictionary where keys are tuples of the form (plant, market). I couldn’t figure out how to specify a R list that would translate to such a dictionary in python. So I used a round about way by converting a tibble to a pandas dataframe to a dictionary.
# shipping cost between plants and markets
costdf = tribble(
~plants, ~mkts, ~dtab,
"seattle", "new-york", 2.5,
"seattle", "chicago", 1.7,
"seattle", "topeka", 1.8,
"san-diego", "new-york", 2.5,
"san-diego", "chicago", 1.8,
"san-diego", "topeka", 1.4
)
costdf = costdf %>% mutate(cost = dtab * 90)
dict_from_df = function(df) {
names(df) = c("x", "y", "z")
dfpy = r_to_py(df)$set_index(c("x", "y"))
dfdict = dfpy$to_dict()["z"]
return(dfdict)
}
cost = dict_from_df(costdf %>% select(plants, mkts, cost))
cost
## {('seattle', 'new-york'): 225.0, ('seattle', 'chicago'): 153.0, ('seattle', 'topeka'): 162.0, ('san-diego', 'new-york'): 225.0, ('san-diego', 'chicago'): 162.0, ('san-diego', 'topeka'):<PHONE_NUMBER>9999999}
Next, we create a model object and define the parameters and decision variables
# create the model object
M = pyo$ConcreteModel()
# set the model parameters
M$plants = pyo$Set(initialize = c("seattle", "san-diego"))
M$mkts = pyo$Set(initialize = c("new-york", "chicago", "topeka"))
M$cap = pyo$Param(M$plants, initialize = list("seattle" = 350, "san-diego" = 600))
M$dem = pyo$Param(M$mkts, initialize = list("new-york" = 325, "chicago" = 300, "topeka" = 275))
M$cost = pyo$Param(M$plants, M$mkts, initialize = cost)
# define the model decision variables (shipment quantities between a plant and market)
# bounds specify that x>=0
M$x = pyo$Var(M$plants, M$mkts, bounds = tuple(0, NULL))
Next, we will set the constraint on supply from plants being below their capacities. In pyomo, the constraints can be described with a function that describes the rule to construct the constraint. In python pyomo, the supply constraint will be listed as follows:
def supply_rule(m, i):
expr = sum(m.x[i, j] for j in mkts) <= cap[i]
return expr
M.supply_cons = pyo.Constraint(plants, rule = supply_rule)
The function supply_rule
is a function that tells what the constraint expression is for a particular plant \(i\). Then M.supply_cons
line ensures that the constraint is generated for each of the plants. Python shines here since this is represented succintly using list comprehensions. Since I am not sure how to do an equivalent representation using R, I used loops in the function to construct the constraint. I also wasn’t sure if using R function to represent the rule will work when passing the R function to create the supply constraint. But it worked fine and I guess that’s the magic of reticulate. Below is the supply constraint constructed in R.
# supply from plant <= plant capacity
supply_rule = function(m, i) {
supply = 0
for (j in bi$list(m$mkts)) {
supply = supply + m$x[tuple(i, j)]
}
expr = supply <= m$cap[[i]]
return(expr)
}
M$supply_cons = pyo$Constraint(M$plants, rule = supply_rule)
We can check if the supply constraint is constructed correctly
r.M.supply_cons.pprint()
## supply_cons : Size=2, Index=plants, Active=True
## Key : Lower : Body : Upper : Active
## san-diego : -Inf : x[san-diego,new-york] + x[san-diego,chicago] + x[san-diego,topeka] : 600.0 : True
## seattle : -Inf : x[seattle,new-york] + x[seattle,chicago] + x[seattle,topeka] : 350.0 : True
Similarly, we can construct the demand constraint
# supply to a market meets its demand
demand_rule = function(m, j) {
demand = 0
for (i in bi$list(m$plants)) {
demand = demand + m$x[tuple(i, j)]
}
expr = demand >= m$dem[[j]]
return(expr)
}
M$demand_cons = pyo$Constraint(M$mkts, rule = demand_rule)
r.M.demand_cons.pprint()
## demand_cons : Size=3, Index=mkts, Active=True
## Key : Lower : Body : Upper : Active
## chicago : 300.0 : x[seattle,chicago] + x[san-diego,chicago] : +Inf : True
## new-york : 325.0 : x[seattle,new-york] + x[san-diego,new-york] : +Inf : True
## topeka : 275.0 : x[seattle,topeka] + x[san-diego,topeka] : +Inf : True
Next, we construct the objective of total shipment cost and set it to be minimized.
# objective to minimize shipping cost
objective_rule = function(m) {
totcost = 0
for (i in bi$list(m$plants)) {
for (j in bi$list(m$mkts)) {
totcost = totcost + m$cost[tuple(i, j)] * m$x[tuple(i, j)]
}
}
return(totcost)
}
M$objective = pyo$Objective(rule = objective_rule, sense = pyo$minimize)
Since this a linear program, we use the open-source solver glpk to solve this.
# solve using solver glpk
opt = pyo$SolverFactory('glpk')
results = opt$solve(M, logfile = "tmp.log")
Optimal shipping cost based on solving this problem is:
r.M.objective()
## 153675.0
Next we display the shipment quantities. It is a bit more verbose since we worked with the option of convert = FALSE
resulting in explicitly extracting desired values.
res_L = list()
k = 1
for(i in bi$list(M$plants)) {
for(j in bi$list(M$mkts)) {
res_L[[k]] = list(plant = i, mkt = j, qty = py_to_r(M$x[tuple(i, j)]()))
k = k + 1
}
}
res_df = bind_rows(res_L)
res_df
## # A tibble: 6 x 3
## plant mkt qty
## <chr> <chr> <dbl>
## 1 seattle new-york 0
## 2 seattle chicago 300
## 3 seattle topeka 0
## 4 san-diego new-york 325
## 5 san-diego chicago 0
## 6 san-diego topeka 275
There are 3 more examples where I developed the pyomo model in R. I have listed them below along with location of the code.
Here I covered two examples to show how to develop a pyomo model from R using the reticulate package. While it might still be easier to develop the pyomo model in python (since it was meant to be that way), I found that it is possible to develop pyomo models in R also fairly easily albeit with some modifications (some maybe less elegant compred to the python counterpart). It may still be better to develop more involved pyomo models in python but reticulate offers a way to develop simple to intermediate levels models directly in R. I am summarizing key learnings:
convert = FALSE
to retain pyomo objects as python objects potentially avoid issues that are hard to troubleshoot.