Hello all
This is my first posting for some years. I am back
using R again and must say I do like the language
(regarding scripting, I also use matlab, perl, and bash).
My question involves plotting a Pareto frontier in
three dimensions. This is strictly a exercise in
visualization, I make no attempt to extract the Pareto
set (aka dominating subset) first.
EXAMPLE PLOTS
For some example plots, please see:
http://www.iet.tu-berlin.de/~morrison/rhelp/pareto-static.png
http://www.iet.tu-berlin.de/~morrison/rhelp/pareto-rotating.gif
MORE SPECIFICALLY
My question is whether these plots duplicate existing
work? And if not, whether my R code (see below) could
or should be placed somewhere more organized for
community usage. The license would be GNU GPLv3
(unless the CC BY-NC-SA or something else is generally
considered more appropriate).
EXISTING R PACKAGES
A search of R packages turned up the following leads:
'mco::paretoFront' : Extract the pareto front or
pareto set from an mco result object.
'TunePareto::plotParetoFronts2D' : Draws a classical
Pareto front plot of 2 objectives of a parameter
tuning
Neither package appears to support three-dimensional
Pareto frontier plots.
BACKGROUND
My PhD project involves simulating a number of energy
system scenarios and then consolidating the costs
across 5-dimensions. At this stage, normally only
3-dimensions are used: financial cost in $, greenhouse
cost in kg CO2e, and depletion cost in J HHV. NOx in
kg and land-use change in m2 are also supported but are
not typically invoked. (CO2e = carbon dioxide
equivalent, HHV = higher heating value, NOx mono-nitrogen oxides.)
These results are collected in a text file and then
read into an R data frame:
leta scenario.name fin ghg nox dep luc
1 + reference energy system 0 0.000 0 0.00 0
2 a indoor temp now 18C was 21C -10 -0.034 0 -0.66 0
3 b far-from-demand ranked-orig-domains 0 0.000 0 0.00 0
4 c ambient air 2C warmer 0 -0.024 0 -0.48 0
5 d uprated building performance 10 -0.089 0 -1.76 0
6 e carbon capture added 190 -1.599 0 20.77 0
7 f nicad storage removed -10 0.000 0 0.00 0
8 g large HV windfarm removed 0 0.003 0 0.06 0
9 h urban domain under 'fin' was 'ghg' 0 0.000 0
0.00 0
10 i zero emissions smelter (experimental) 0 -0.726 0 0.00 0
11 j reduced inflows (dry year) 0 0.159 0 3.17 0
The 'leta' is simply my single character scenario
identifier. The following attributes are then added to
the data frame to facilitate the later annotation of
the plot, in summary:
vec unit label
------------------------
fin $M financial decrease
ghg Tg greenhouse decrease
nox kg NOx decrease
dep PJ depletion decrease
luc m2 landuse decrease
The Pareto frontier captures that subset of scenarios,
dependent on their high-level performance, that are
said to "dominate" and which, in the absence of other
criteria, should alone be examined and traded-off. The
remaining scenarios need not be considered. This
dominating subset is also known as the "Pareto set".
It is possible that one scenario dominates all others
and the Pareto set has only one element.
The problem of identifying the Pareto set is called the
"maximal or maximum vector problem" and is well studied
in computer science.
The Pareto set can also be visualized in 3-space by
plotting the points and then constructing the Pareto
frontier. Moreover is is straightforward to generate
this frontier using the 'rgl' package. And this is
what I do.
LITERATURE
Three-dimensional Pareto frontier visualization is
reported in the literature. For instance (not
paywalled -- thanks guys!):
Madetoja, Elina, Henri Ruotsalainen, Veli-Matti
M?nkk?nen, Jari H?m?l?inen, and Kalyanmoy Deb.
2008. Visualizing multi-dimensional Pareto-optimal
fronts with a 3D virtual reality system. IEEE
Proceedings. ISBN-13 978-83-60810-14-9.
Madetoja etal (2008) also cite earlier papers,
published in 2004 and 2005.
R CODE
<---- R code starts -------------------------------------------->
# Waiver: To the extent possible under law, Robbie
# Morrison <robbie at actrix.co.nz> has waived all
# copyright and related or neighboring rights to this
# R scripting. This work is published from Germany.
# http://creativecommons.org/publicdomain/zero/1.0/
# the various entries in 'opt' are set using command-line arguments
# ---------------------------------
# step 00 : plot preparation
# ---------------------------------
# at the outset define some functions: three for
# displaying shrouds and two for plotting axes
# ---------------------------------
# function : quad
# ---------------------------------
# description : generate quadrilateral from four "flat" vertices
# role : called by function 'rect'
# caution : the 'OpenGL' graphics library uses homogeneous,
# not Euclidean, coordinates
# status : complete
#
# Typical input uses Euclidean coordinates
#
# vertices (3 x 4) : 21.4 0.77 15.3 21.4 0 15.3 21.4 0.77 0 21.4 0 0
# indices (4) : 1 2 4 3
#
# Coordinate systems
#
# This material is not as relevant after adopting "homogeneous = F"
#
# http://www.mail-archive.com/r-help at r-project.org/msg58482.html
#
# One would think that each vertex would have 3
# coordinates (x,y,z) what does the fourth one in
# the definition of the variable vertices stand
# for.
#
# By default qmesh3d uses 4-coordinate homogeneous
# coordinates, because that's the coordinate system
# used by OpenGL. See the ?rgl::matrices help topic
# for a description of how they work.
#
# Report formatting
#
# Note that -1.0 * 0.0 gives -0.0 which 'formatC'
# then outputs as "-0". This is quite normal under
# floating point arithmetic.
#
# ---------------------------------
quad <- function(vertices, # four 3-space vertices
indices) # stated ordering for above
{
# report
if ( opt$debug )
{
cat(" vertices", formatC(unlist(vertices), width = 9))
cat("\n")
}
# active calls
object <- qmesh3d(unlist(vertices),
indices,
homogeneous = F) # 'rgl' quad-mesh function
shade3d(object,
col = quad.color,
alpha = quad.alpha)
}
# ---------------------------------
# function : rect
# ---------------------------------
# description : generate four 3-space coordinates to pass to function
'quad'
# role : utility call
# comment : hardcoded (not the most elegant method)
# status : complete
# ---------------------------------
rect <- function(point, # data point
ofset, # transposed origin
focus) # "x" "y"
"z"
{
switch(focus,
x = {
p1 <- c(point[1], point[2], point[3])
p2 <- c(point[1], ofset[2], point[3])
p3 <- c(point[1], ofset[2], ofset[3])
p4 <- c(point[1], point[2], ofset[3])
},
y = {
p1 <- c(point[1], point[2], point[3])
p2 <- c(ofset[1], point[2], point[3])
p3 <- c(ofset[1], point[2], ofset[3])
p4 <- c(point[1], point[2], ofset[3])
},
z = {
p1 <- c(point[1], point[2], point[3])
p2 <- c(ofset[1], point[2], point[3])
p3 <- c(ofset[1], ofset[2], point[3])
p4 <- c(point[1], ofset[2], point[3])
},
stop("coding error: focus not supported", focus))
list(p1, p2, p3, p4) # positive spin (right-hand
rule) relative to 'focus' axis
}
# ---------------------------------
# function : shroud
# ---------------------------------
# description : create a "shroud" over 'point' relative to
'ofset'
# role : utility call
# status : complete
# ---------------------------------
shroud <- function(point, # 3-space point of interest
ofset) # 3-space transposed origin
{
if ( opt$debug ) message(" point :", format(point, digits = 2,
width
= 6))
if ( opt$debug ) message(" ofset :", format(ofset, digits = 2,
width
= 6))
indices <- c(1, 2, 3, 4) # CAUTION: ordering is
significant (else get "twisted bow-tie" shapes)
quad(rect(point, ofset, "x"), indices)
quad(rect(point, ofset, "y"), indices)
quad(rect(point, ofset, "z"), indices)
}
# ---------------------------------
# function : line
# ---------------------------------
# description : generate line from two vertices
# role : called by function 'rule'
# status : complete
# ---------------------------------
line <- function(vertices) # two 3-space vertices
{
# report
if ( opt$debug )
{
cat(" vertices", formatC(unlist(vertices), width = 9))
cat("\n")
}
# active calls
lines3d(matrix(unlist(vertices), ncol = 3, byrow = T),
col = line.color,
lwd = line.weight,
homogeneous = F) # 'rgl' line function
}
# ---------------------------------
# function : rule
# ---------------------------------
# description : generate two 3-space coordinates to pass to function
'line'
# role : utility call
# status : complete
# ---------------------------------
rule <- function(range, # data range
const, # constant dimensions
focus, # "x" "y"
"z"
add = 0.0)
{
if ( add != 0.0 )
{
len <- abs(range[1] - range[2])
adj <- add * len
range[1] <- range[1] - adj
range[2] <- range[2] + adj
}
switch(focus,
x = {
p1 <- c(range[1], const[1], const[2])
p2 <- c(range[2], const[1], const[2])
},
y = {
p1 <- c(const[2], range[1], const[1])
p2 <- c(const[2], range[2], const[1])
},
z = {
p1 <- c(const[1], const[2], range[1])
p2 <- c(const[1], const[2], range[2])
},
stop("coding error: focus not supported", focus))
list(p1, p2)
}
# ---------------------------------
# step 00 : now plot the data
# ---------------------------------
device.no <- open3d()
bg3d(background.color)
# regarding points: use "size = 10" if 'type' not set
plot3d(data[vecs],
main = main,
xlab = xlab,
ylab = ylab,
zlab = zlab,
col = point.color,
type = "s", # CAUTION: 'type'
not 'pch',
"s" for spheres
size = 0.7,
box = F, # added in below if required
axes = F) # added in below if required
if ( opt$axes == 0 ) axes3d() # add back standard
axes
if ( opt$axes == 1 ) axes3d(c('x--', 'y--', 'z-+')) #
used fixed
specified axes
if ( opt$axes == 2 ) axes3d(c('x+-', 'y--', 'z++')) #
used fixed
specified axes
if ( opt$axes == 3 ) invisible() # do nothing
# extract the offset (looping unavoidable)
offset <- c()
if ( opt$zero )
{
offset <- rep(0.0, 3) # the origin
} else {
for ( vec in vecs ) # usually 'fin'
'ghg' 'dep'
{
raw <- data[, vec]
off <- max(raw)
skirt <- opt$move * (max(raw) - min(raw))
off <- off + skirt
offset <- c(offset, off) # collect
}
}
message("offset (transposed origin) =", format(offset, digits = 2,
width 8))
message()
# add the pareto frontier
if ( ! opt$unpareto )
{
for ( i in 1:nrow(data) )
{
sen <- unFactor(data$leta)[i] # 'sen' is a string,
'unFactor' is a local unfactor function
if ( opt$debug ) message(sprintf(" i = %2d scenario = %s",
i,
sen))
point <- unlist(data[i, vecs]) # each row with selected cols
# active calls
quad.color <- rainbow(20)[i]
shroud(point, offset) # shroud creation call
if ( ! opt$nolabels )
{
text3d(point, # point labeling call
text = sen,
adj = label.adjust,
col = label.color)
}
}
if ( opt$debug ) message()
}
# add an origin marker
if ( opt$origin )
{
message("adding origin marker")
if ( opt$debug ) message()
axes <- c("x", "y", "z")
for ( i in 1:3 )
{
raw <- data[, vecs[i]]
ran <- range(raw)
const <- c(0.0, 0.0)
line(rule(ran, const, axes[i], opt$move))
}
}
# add a transposed origin marker
if ( opt$quip )
{
message("adding transposed origin marker")
if ( opt$debug ) message()
offsets <- c(offset, offset) # useful to be able to roll
around
axes <- c("x", "y", "z")
for ( i in 1:3 )
{
raw <- data[, vecs[i]]
ran <- range(raw)
const <- c(offsets[i + 1], offsets[i + 2])
line(rule(ran, const, axes[i], opt$move))
}
}
if ( opt$origin || opt$quip ) message()
# finally add the bounding box
if ( opt$wire )
{
box3d()
}
<---- R code ends ---------------------------------------------->
all the best
---
Robbie Morrison
PhD student -- policy-oriented energy system simulation
Institute for Energy Engineering (IET)
Technical University of Berlin (TU-Berlin), Germany
University email (redirected) : morrison at iet.tu-berlin.de
Webmail (preferred) : robbie at actrix.co.nz
[from Webmail client]
