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

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. 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.

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.

smoothing hydrographic profiles

Abstract.Methods are presented for smoothing hydrographic profiles.

1. Introduction

Smoothing hydrographic profiles with conventional time-series methods is problematic for two reasons: (a) the data are commonly not equi-spaced in depth and (b) the data seldom lack trends in depth. The issues and their solutions are illustrated without much discussion here.

2. Methods

library(oce)
library(signal)
data(ctd)
S <- ctd[['salinity']]
p <- ctd[['pressure']]
## must create equispaced data for filtering to make sense
dp <- median(diff(p))
pp <- seq(min(p), max(p), dp)
S0 <- approx(p, S, pp)$y W <- dp / 2 # critical frequency f1 <- butter(1, W) f2 <- butter(2, W) par(mfrow=c(1, 3)) ## filter raw profile plotProfile(ctd, xtype="salinity", type='l') S0f1 <- filtfilt(f1, S0) S0f2 <- filtfilt(f2, S0) lines(S0f1, pp, col='red') lines(S0f2, pp, col='blue') mtext("(a) ", side=3, adj=1, line=-5/4, cex=3/4) ## filter detrended profile plotProfile(ctd, xtype="salinity", type='l') Sd <- detrend(pp, S0) S1f1 <- filtfilt(f1, Sd$Y) + Sd$a + Sd$b * pp
S1f2 <- filtfilt(f2, Sd$Y) + Sd$a + Sd$b * pp lines(S1f1, pp, col='red') lines(S1f2, pp, col='blue') mtext("(b) ", side=3, adj=1, line=-5/4, cex=3/4) ## smooth-spline raw profile spline <- smooth.spline(pp, S0, df=3/W) # suggestion: try different df values S2 <- predict(spline)$y
plotProfile(ctd, xtype="salinity", type='l')
lines(S2, pp, col='red')
mtext("(c) ", side=3, adj=1, line=-5/4, cex=3/4)


3. Results

The first order filter is less prone to isolated wiggles than the second order one, but both behave poorly near the top and bottom of the profile, unless the data are detrended. With detrending, filtered signals are similar to those calculated with the smoothing spline.

A salinity profile (black line) with various smoothing models. (a) Smoothing by first-order Butterworth filter (red) and second-order Butterworth filter (blue) on raw data. (b) As (a) but using detrended data. (c) Spline smoothing.

4. Discussion and conclusions

With detrending, the filtered results can be made similar to those of smoothing spline. However, splines have some practical advantages: (1) they do not require equispaced data, (2) they do not require detrending, (3) the results are inherently smooth (by construction of the spline) and (4) calculating derivatives is easy with a spline. These advantages explain why splines are used to calculate buoyancy frequency in the Oce package.

inferring halocline depth

Abstract. A method given for inferring halocline depths based on derivatives calculated with a smoothing spline.

1. Introduction.

There are no agreed-upon methods for inferring halocline depth, but a reasonable method might involve locating the depth at which dS/dp is largest, where S is salinity and p is pressure (Kelley 2014 chapter 5).  Calculating the derivative using e.g. diff(S)/diff(p) can be problematic because of sensitivity to noise, especially for data that have not been bin-averaged. Achieving smoothness with conventional filtering has problems at the end-points, which is particularly troublesome for a near-surface halocline (see the next blog entry). A possible solution to such problems is to calculate the derivative with a smoothing spline.

2. Methods.

Pasted below is test code that does this with the ctd dataset in the oce package. The function returns the pressure at which the smoothing spline has highest salinity derivative, and it can also plot the results (which is recommended). The parameter named deltap is used to set the value of df (degrees of freedom) for the spline. One might think of deltap as the thickness (in dbar) of the smoothing interval for each of the sub-components of the spline.

library(oce)

findHalocline <- function(ctd, deltap=5, plot=TRUE)
{
S <- ctd[['salinity']]
p <- ctd[['pressure']]
n <- length(p)
## trim df to be no larger than n/2 and no smaller than 3.
N <- deltap / median(diff(p))
df <- min(n/2, max(3, n / N))
spline <- smooth.spline(S~p, df=df)
SS <- predict(spline, p)
dSSdp <- predict(spline, p, deriv=1)
H <- p[which.max(dSSdp$y)] if (plot) { par(mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0)) plotProfile(ctd, xtype="salinity") lines(SS$y, SS$x, col='red') abline(h=H, col='blue') mtext(sprintf("%.2f m", H), side=4, at=H, cex=3/4, col='blue') mtext(sprintf(" deltap: %.0f, N: %.0f, df: %.0f", deltap, N, df), side=1, line=-1, adj=0, cex=3/4) } return(H) } # Plot two panels to see influence of deltap. par(mfrow=c(1, 2)) data(ctd) findHalocline(ctd) findHalocline(ctd, 1)  3. Results. The graph shows results for a test dataset provided with the oce package, using two values for deltap. 4. Discussion. Lacking further information about the data or a statistical model of the variation, the choice of deltap is somewhat open, and one interpretation of the results would be to state that the halocline depth is likely to be between 4 and 5 m. A spline-based method for halocline detection. Left: with default parameters. Right: setting the smoothing interval to 1 dbar (roughly 1 m). 5. Conclusions. Further work is required to test and calibrate the method with other data. References. 1. Kelley, Dan, 2014. Oceanographic Analysis with R. Springer-Verlag. In preparation. installing the oce package for the R language Several of the blog items have used the oce package. The official version of this can be installed from within R by install.packages("oce")  and more up-to-date versions can be installed using the devtools package written by Hadley Wickham, which is itself installed with install.packages("devtools")  after which installing the latest development version of oce is accomplished with library(devtools) install_github('ocedata', 'dankelley', 'master') install_github('oce', 'dankelley', 'develop')  Note that ocedata does not need to be updated frequently, as it only updated when new datasets are added to oce. The official version of oce is only updated every few months, but the branch named develop (used above) may be updated several times a day, if the author is adding new features or fixing bugs. For more on oce, see the oce website on github. R interface to WebTide (tidal prediction) The previous posting explained how to install WebTide on an OSX machine. This one shows how to hook up to an installed WebTide database, so that R code can get tidal predictions. The following code in the R language will produce a graph in which the top panel mimics the tidal-elevation graph produced by WebTide itself (see previous blog posting for comparison). library(oce) tStart <- as.POSIXct("2013-12-29 14:21:00", tz="UTC") tEnd <- as.POSIXct("2014-01-13 15:21:00", tz="UTC") time<-seq(tStart, tEnd, by=15, units="minutes") prediction <- webtide("predict", lon=-65.06747, lat=45.36544, time=time)  A webtide prediction using the oce package of the R library One of the advantages of accessing the tidal prediction from within oce is to make it easier to undertake further analysis, e.g. a node nearer Halifax has a mixed tide, with the following illustrating in terms of velocity and a so-called progressive vector p <- webtide("predict", node=14569) par(mfrow=c(2,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0)) plot(p$u, p$v, asp=1, type="o") dt <- diff(as.numeric(p$time[1:2]))
x <- dt * cumsum(p$u) y <- dt * cumsum(p$v)
plot(x, y, asp=1, type="o")


Velocity covariance and progressive vector diagram for a mixed tide.

WebTide (tidal prediction) installation (OSX)

Abstract.This blog item explains how to install WebTide, a tidal prediction application, on an OSX machine. The work requires moderate expertise, being carried out in the console, and requiring that a C compiler be found on the host machine.

1. Introduction

WebTide comes with a friendly graphical user interface that makes it easy to get tidal predictions at various locations around the world. Anyone with an interest in oceanography is likely to find it a useful application. Although the interface is slightly quirky (particularly for map navigation), it only takes a few minutes to learn. (An upcoming blog entry will explain how to avoid the WebTide interface, using the oce R package to access the data that underly WebTide.)

2. Installation

1. Download the “linux without Java VM” version from the Bedford Institute of Oceanography website, and store the results in the ~/Downloads directory.

2. Open a terminal or console window, and type the following to enter the WebTide directory and start the process of installing WebTide.

cd ~/Downloads
sh WebTide_unix_0_7.sh


After a moment, a dialog box will pop up, asking where WebTide is to be installed. To install a version just for you, supply /Users/USERNAME/WebTide, where USERNAME is your username. To install a version for all accounts on the computer, use the default, which is /usr/local/WebTide. A second dialog box addresses the issue of symlinks. Click the box indicating that these are not to be used. After a few moments, a final dialog box will appear, stating that the work has been completed. Click the button named “Finish”.

3. Next comes the OSX-specific part of the work. You will need the gcc compiler to do this. It is available for free on the Apple website, as the Xcode development tools. (If you have been using the computer for science, it is quite likely you have already installed a compiler, so just try the steps below, before bothering with a re-install of Xcode.) In a console, type the following to enter the source directory (here assumed to be configured for a single user) and to compile two programs.

cd ~/WebTide/Tidecor_src
gcc -O2 -o tidecor webtidecor2.5.1.c -lm
gcc -O2 -o constituentinterpolator constituentinterpolator1.1.c -lm


In the two gcc lines shown above, it may be necessary to change the names of the ".c" files, to match whatever you find in this tidecor_src directory. You may ignore the warning messages that are likely to appear, but if errors appear, you will need to interrupt this tutorial to find out what is wrong. Possibly leaving a comment on this site would encourage another reader to offer a solution; possibly the scientists at the Bedford Institute of Oceanography (named on the website) could help.

4. Finally, move these newly-compiled programs from the source directory into the bin directory, by typing

mv tidecor ../bin
mv constituentinterpolator ../bin


3. Using WebTide

If you installed WebTide to your home directory, launch WebTide by typing the following in a console window

~/WebTide/WebTide


and use the GUI from that point on. If you installed it to /usr/local/WebTide, type the following instead

/usr/local/WebTide/WebTide


The main WebTide window looks like the following.

A WebTide map, showing a marker in the Bay of Fundy (see also the sea-level elevation diagram)

There is little point in explaining the details of the menu and mouse actions here, since WebTide provides its own documentation. In most cases, the action will be to click the mouse on spot on the map, to set a so-called “marker” for a location at which calculations are to be done. The diagram shown here has a marker in the Bay of Fundy. Once a marker has been set, the user may wish to click a menu item (see at the right-hand side) to get a tidal prediction, as shown in the screenshot below.

A WebTide prediction of sea-level elevation