Author Archives: dankelleyns

one-dimensional optimization

Abstract. The method of one-dimensional optimization is explained in the context of capillary-gravity waves.

1. Introduction.

R provides functions for both one-dimensional and multi-dimensional optimization. The second topic is much more complicated than the former (see e.g. Nocedal 1999) and will be left for another day.

A convenient function for 1D optimization is optimize(), also known as optimise(). Its first argument is a function whose minimum (or maximum) is sought, and the second is a two-element vector giving the range of values of the independent variable to be searched. (See ?optimize for more.)

2. Application.

As an example, consider the phase speed of deep gravity-capillary waves, which is given by \omega/k where \omega is the frequency and k is the wavenumber, and the two are bound together with the dispersion relationship \omega^2=gk+\sigma k^3/\rho, where g is the acceleration due to gravity, \sigma is the surface tension parameter (0.074 N/m for an air-water interface) and \rho is the water density (1000 kg/m^3 for fresh water). This yields wave speed given by the following R function.

phaseSpeed <- function(k) {
    g <- 9.8
    sigma <- 0.074  # water-air
    rho <- 1000  # fresh water
    omega2 <- g * k + sigma * k^3/rho
    sqrt(omega2)/k
}

Readers with a background in the topic of waves may know that there is a minimum phase speed at wavelengths \lambda of about 0.02m, or a wavenumber k=2\pi/\lambda of roughly 300. It always makes sense to plot a function to be optimized, if only to check that it has been coded correctly, so that's the next step. We’ll use a range of half an order of magnitude (factor of 3) smaller and larger k.

k <- seq(100, 1000, length.out = 100)
par(mar = c(3, 3, 1, 1), mgp = c(2, 0.7, 0))
plot(k, phaseSpeed(k), type = "l",
     xlab = "Wavenumber", ylab = "Phase speed")
Wave speed for capillary-gravity waves in infinitely deep water

Wave speed for capillary-gravity waves in infinitely deep water

The results suggest that the range of k illustrated in the diagram contains the minimum, so we provide that to optimize().

o <- optimize(phaseSpeed, range(k))
phaseSpeed(o$minimum)
## [1] 0.2321

This speed is not especially fast; it would take about a heartbeat to move past your hand. (Or, about as fast as the waves that pulse your heart … presumably with different physics, though!)

3. Exercises

  1. Use str(o) to learn about the contents of the optimized solution.

  2. Use abline() to indicate the wavenumber at the speed minimum.

  3. Try other functions that are of interest to you.

  4. Use the multi-dimensional optimizer named optim() on this problem.

References

Jorge Nocedal and Stephen J. Wright, 1999. Numerical optimization. Springer series in operations research. New York, NY, USA.

using the plyr package in R

Abstract. The merits of the plyr package are illustrated.

1. Introduction.

The base R system provides lapply() and related functions, and the package plyr provides alternatives that are worth considering. It will be assumed that readers are familiar with lapply() and are willing to spend a few moments reading the plyr documentation, to see why the illustration here will use the ldply() function.

The test task will be extraction of latitude (and then both latitude and longitude) from the section dataset in the oce package. (Users of that package may be aware that there is a built-in accessor for doing this, so results can easily be checked.)

2. Methods

First, load the data

library(oce)
data(section)

Next, find latitudes using lapply

lat <- unlist(lapply(section[["station"]], function(x) x[["latitude"]]))

Next, find latitudes with ldply

library(plyr)
lat <- ldply(section[["station"]], function(x) x[["latitude"]])

3. Results

The reader can check that the results match, although ldply() returns a data frame, not a vector as in the first method. Tests of speed indicate such a difference too small to be of much interest, at least in this case

library(microbenchmark)
microbenchmark(ldply(section[["station"]], function(x) x[["latitude"]])$V1, 
    unlist(lapply(section[["station"]], function(x) x[["latitude"]])))
## Unit: milliseconds
##                                                               expr   min
##        ldply(section[[&quot;station&quot;]], function(x) x[[&quot;latitude&quot;]])$V1 18.99
##  unlist(lapply(section[[&quot;station&quot;]], function(x) x[[&quot;latitude&quot;]])) 18.36
##     lq median    uq   max neval
##  20.26  20.56 21.02 36.05   100
##  19.71  19.93 20.64 63.18   100

.

4. Discussion

Since ldply() returns a data frame, it is more flexible than unlist(), which returns a vector. For example, the following creates a data frame with columns for lat and lon:

latlon <- ldply(section[["station"]], function(x) c(x[["latitude"]], x[["longitude"]]))

A station plot is produced as follows.

if (!interactive()) png("plyr.png", unit = "in", 
    res = 150, width = 7, height = 7, pointsize = 10)
mapPlot(coastlineWorld, projection = "orthographic",
     orientation = c(20, -40, 0))
mapPoints(latlon$V2, latlon$V1, pch = "+", 
    cex = 1/2, col = "red")
if (!interactive()) dev.off()
Map of station locations, produced by extracting locations using a function from the plyr package

Map of station locations, produced by extracting locations using a function from the plyr package

5. Conclusions

The effort of learning how to use the plyr package is likely to pay off in more flexible code, particularly because of the use of data frames in that package. On this theme, note that the author of plyr is developing a similar package called dplry, which centres more closely on data frames and offers many new features; see http://blog.rstudio.org/2014/01/17/introducing-dplyr/ for a blog item introducing dplyr.

using R markdown on wordpress

Abstract. An explanation is given for using the vim editor to compose wordpress blog items in normal-markdown or R-markdown format.

1. The gist of R markdown format.

R markdown is an extension of normal-markdown that permits the embedding of R code and output (including figures). WordPress seems not to handle the figure-inclusion mechanism, but there is a workaround: create the figures locally then use the ‘insert media’ button on wordpress to upload and include them in the text. For more details of R-markdown format, see e.g. the rstudio documentation.

2. Setting up wordpress to accept markdown format.

As noted in a comment on this posting, go to the “Dashboard” and then select “Settings”, then select “Writing” and then click the box next to “Use markdown for posts and pages”.

3. Guidelines for vim editor.

WordPress will insert line breaks whenever they are found in the text, so it is important not to use auto-fill in Vim. This can be set up by inserting the line

autocmd! BufNewFile,BufReadPre,FileReadPre *.Rmd so ~/.vim/rmd.vim

in the ~/.vimrc file, and then creating a new file in ~/.vim named rmd.vim that contains

" wrap long lines in the editor window ...
set wrap
" ... and prevent automatic end-of-line insertion
set textwidth=0 wrapmargin=0

4. Guidelines for emacs editor.

(Something may be added here if someone tells me what to add…)

5. Converting R markdown to normal markdown

R-markdown may contain R code chunks (with the token {r} after the triple-backtick). When these are run through the knit() function provided in the knitr package, the braces are removed, which is the standard for syntax-enabled markdown as used on wordpress.

Conversion from R markdown to normal markdown can be done easily inside R, with the following (in this case for the file used for the present posting).

library(knitr)
knit("editing_vim.Rmd")

colourizing along a trajectory

Abstract. Demonstrate how to colourize an (x,y) curve according to z.

1. Introduction.

In Oceanography it can be useful to use colour to display z values along an (x,y) trajectory. For example, CTD data might be displayed in this way, with x being distance along track, y being depth, and z being temperature. This post shows how one might do this generally with fake data.

2. Methods.

The R code given below demonstrates this with fake data. The core idea is to use segments(), here with head() and tail() to chop up the trajectory.

library(oce)
x <- 1:50
y <- x * 2/10 + x^2/10
z <- seq(0, 5, length.out = length(x))
zlim <- range(z)
npalette <- 10

palette <- oceColorsJet(npalette)
drawPalette(zlim = zlim, col = palette)
plot(x, y, type = "l")
segments(head(x, -1), head(y, -1),
  tail(x, -1), tail(y, -1),
  col = palette[findInterval(z, 
    seq(zlim[1], zlim[2], length.out = npalette))],
  lwd = 10)
points(x, y, pch = 20)

3. Results.

The graph shows that this works reasonably well. The dots will probably not be useful in an actual application, and are used here just to indicate colourization in groups.

Demonstration of colourizing an (x,y) trajectory according to the value of a third variable, z.

Demonstration of colourizing an (x,y) trajectory according to the value of a third variable, z.


4. Conclusions.

The method works well, and is flexible in terms of colour schemes.

Appendix 1: exercises find a way to blend colour between the points, perhaps by defining a euclidian distance in (x,y) space (which will of course require scaling if x and y have different units or ranges) and then using approx().

overshoot in butterworth low-pass filters

Abstract. Demonstrate overshoot for butterworth filters in R.

1. Introduction.

Butterworth filters with order other than 1 have an overshoot phenomenon that can be problematic in some cases. For example, if smoothing is used on an estimate of kinetic energy, overshoots might yield negative values that are nonphysical. This post simply illustrates this with made-up data that the reader can experiment with.

2. Methods.

The code given below creates the graph shown in the next section

library(signal)
n <- 100
x <- 1:n
y <- ifelse(0.3*n < x & x < 0.7*n, 1, 0)
if (!interactive()) png("butter.png", width=7, height=5, unit="in", res=150, pointsize=9)
plot(x, y, type='o', pch=20, ylim=c(-0.1, 1.1))
b1 <- butter(1, 0.1)
y1 <- filtfilt(b1, y)
points(x, y1, type='o', pch=20, col='red')

b2 <- butter(2, 0.1)
y2 <- filtfilt(b2, y)
points(x, y2, type='o', pch=20, col='blue')

b3 <- butter(3, 0.1)
y3 <- filtfilt(b3, y)
points(x, y3, type='o', pch=20, col='forestgreen')
if (!interactive()) dev.off()

3. Results

The graph illustrates; for the colour coding, see the code above.

Demonstration of first, second, and third order butterworth low-pass filters.

4. Conclusions

Caution is recommended in using butterworth filters of order 2 and higher. The advantage of the rapid cutoff in frequency space may be outweighed by the disadvantage of overshooting in time space.

cabelling

Abstract R code is provided in aide of laboratory demonstration of
cabelling.

1. Introduction

Setting up a cabelling experiment requires creating two watermasses of equal density, and if only S and T can be measured, that means calculating densities. Using a TS diagram and graphical interpolation is one approach to that task, but another is to use R to do the calculation.

2. Methods

The code given below will do the calculation for specified Sa, Ta and Sb, where the second letter indicates the watermass. The code uses uniroot() to find the temperature Tb that yields equal densities for watermasses “a” and “b”.

# Alter next three lines as desired; a and b are watermasses.
Sa <- 30
Ta <- 10
Sb <- 40
library(oce)
# Should not need to edit below this line
rho0 <- swRho(Sa, Ta, 0)
Tb <- uniroot(function(T) rho0-swRho(Sb,T,0), lower=0, upper=100)$root
Sc <- (Sa + Sb) /2
Tc <- (Ta + Tb) /2
## density change, and equiv temp change
drho <- swRho(Sc, Tc, 0) - rho0
dT <- drho / rho0 / swAlpha(Sc, Tc, 0)
if (!interactive()) png("cabelling.png", width=7, height=7,
                        unit="in", res=200, pointsize=12)
plotTS(as.ctd(c(Sa, Sb, Sc), c(Ta, Tb, Tc), 0), pch=20, cex=2)
drawIsopycnals(levels=rho0, col="red", cex=0)
segments(Sa, Ta, Sb, Tb, col="blue")
text(Sb, Tb, "b", pos=4)
text(Sa, Ta, "a", pos=4)
text(Sc, Tc, "c", pos=4)
legend("topleft",
       legend=sprintf("Sa=%.1f, Ta=%.1f, Sb=%.1f  ->  Tb=%.1f, drho=%.2f, dT=%.2f",
                      Sa, Ta, Sb, Tb, drho, dT),
       bg="white")
if (!interactive()) dev.off()

If run non-interactively, the code will produce a PNG file like that given below.

3. Results

The legend summarizes the results, indicating also the density change and the temperature change that would be equivalent to that density change (at the midpoint, “c”).

TS diagram for the setup of a cabelling experiment.

TS diagram for the setup of a cabelling experiment.

4. Conclusions

If the design goal is that the density mismatch between watermasses “a” and “b” should be, say, 10 percent of the density difference of the mixture watermass (“c”), then in the illustrated case the temperature would have to be controlled to within a quarter of a degree Celcius — a task that is challenging enough to argue against this as an informal classroom demonstration.

Exercises.

  1. Alter the R code to calculate Sb in terms of Tb.
  2. Consider (calculate or measure) the convection associated with sidewall heat conduction.