Problem 2: Immutable objects (30 pts)

Consider the following function which swaps elements i and j of a vector v.

swap <- function(v, i, j) 
  { 
    tmp <- v[i]
    v[i] <- v[j]
    v[j] <- tmp 
  }

A disadvantage of this implementation is that it copies the entire vector, v, in order to work with just two of its elements. Consider the following example.

v <- 1:1000000000
system.time(swap(v, i=1, j=2))
##    user  system elapsed 
##   0.536   0.288   0.826

Your tasks on this problem are the following.

  1. Report on how much time it takes to swap two elements (i=1; j=2) directly on the command line, i.e., without wrapping in a function.

Here is the setup of such an expression, which is quoted below so it can be evaluated for timing purposes later.

i <- 1
j <- 2
expr <- quote({tmp <- v[1]; v[1] <- v[2]; v[2] <- tmp})

Here are the times from five consecutive evaluations.

rbind(system.time(eval(expr)), system.time(eval(expr)), system.time(eval(expr)))
##      user.self sys.self elapsed user.child sys.child
## [1,]     0.496      0.3   0.796          0         0
## [2,]     0.000      0.0   0.000          0         0
## [3,]     0.000      0.0   0.000          0         0
rbind(system.time(swap(v, i=1, j=2)), system.time(swap(v, i=1, j=2)), system.time(swap(v, i=1, j=2)))
##      user.self sys.self elapsed user.child sys.child
## [1,]     0.532    0.280   0.810          0         0
## [2,]     0.504    0.304   0.808          0         0
## [3,]     0.488    0.320   0.808          0         0
  1. Write a new version of the swap function, called swap.eval, which uses quote and eval to perform the calculation just like in part a. but within the function environment and without copying v by working on v in the parent.frame. Although this is a toy example, a similar code might be useful if, say, indicies i and j required substantial pre-calculation within the function before the swap occurred. Demonstrate your swap.eval with i=1; j=2 and report on the time.

The function below takes the same idea behind the quoted function above, but puts it inside a function. The eval call puts the scope of evaluation back in the parent frame so that v need not be copied in at great expense. Note however that all variables must come from the parent frame, including i and j.

swap.eval <- function()
  {
    e <- quote({tmp <- v[i]; v[i] <- v[j]; v[j] <- tmp})
    eval(e, envir=parent.frame())
  }

Now we can try the swap.

rbind(system.time(swap.eval()), system.time(swap.eval()), system.time(swap.eval()))
##      user.self sys.self elapsed user.child sys.child
## [1,]     0.492    0.316   0.809          0         0
## [2,]     0.000    0.000   0.002          0         0
## [3,]     0.000    0.000   0.000          0         0
  1. Write a similar function named swap.do which can be called via do.call that similarly accesses v in the parent frame. Add a print statement at the end of swap.do to show the first five elements of v after the swap occurs. Demonstrate swap.do with i=1; j=2 and report on the time. Are there any disadvantages to swap.do compared to swap.eval?

The function below suites our specifications. However, I have omitted the print statement for now.

swap.do <- function(i, j)
  {
    tmp <- v[i]
    v[i] <- v[j]
    v[j] <- tmp
  }

Lets call it three times and collect timings, just as above.

expr <- quote(do.call("swap.do", args=list(i=1, j=2), quote=TRUE))
rbind(system.time(eval(expr)), system.time(eval(expr)), system.time(eval(expr)))
##      user.self sys.self elapsed user.child sys.child
## [1,]     0.492    0.304   0.795          0         0
## [2,]     0.500    0.296   0.796          0         0
## [3,]     0.488    0.308   0.795          0         0

Now lets add the print statement to explore other behavior.

swap.do <- function(i, j)
  {
    tmp <- v[i]
    v[i] <- v[j]
    v[j] <- tmp
    print(v[1:5])
  }

Consider the state of v before and after the swap.do call …

Before:

print(v[1:5])
## [1] 1 2 3 4 5

During:

do.call("swap.do", args=list(i=1, j=2), quote=TRUE)
## [1] 2 1 3 4 5

After:

v[1:5]
## [1] 1 2 3 4 5

Consider the swap.eval alternative:

swap.eval()
v[1:5]
## [1] 2 1 3 4 5

Problem 3: Bisection broadening (30 pts)

Recall the bisection algorithm and S3 object-oriented implementation from class. The bisection method can be generalized to deal with the case \(f(x_l) f(x_r) > 0\), by broadening the bracket. That is,

\[ \begin{aligned} m &\leftarrow (x_l + x_r)/2, & w &\leftarrow x_r - x_l \\ x_l&\leftarrow m - w, & x_r &\leftarrow m + w. \end{aligned} \]

Your tasks are the following.

  1. Incorporate bracketing into the function we coded. Note that broadening is not guaranteed to find \(x_l\) and \(x_r\) such that \(f(x_l) f(x_r) \leq 0\), so you should include a limit on the number of times broadening is successively tried with a sensible default.

See the revised bisection2.R file, which is a version of the original bisection.R, augmented to include broadening. Some highlights of the modifications include

To load those functions in we have to source the new R file.

source("bisection2.R")
  1. Use your modified function to find a root of the (original) function \(f(x)\) we used in class, but with a different starting interval of \(x_l = 2\) and \(x_r = 3]\), i.e., not containing the root we found in class.

First define the function, as in class.

f <- function(x) log(x) - exp(-x)

Then call the function and demonstrate the print output. Here the default verb=0 is used so as not to clutter up the document, but feel free to try verb=1.

fr <- bisection(f, 2, 3)
fr
## Root of:
## function(x) log(x) - exp(-x)
## <bytecode: 0x5e18630>
## in (2, 3): 1.3098
## found after 2 broading and 29 bracketing iterations
## to a tolerance of 1.490116e-08

The summary output contains more information.

summary(fr)
## Root of:
## function(x) log(x) - exp(-x)
## <bytecode: 0x5e18630>
## in (2, 3): 1.3098
## found after 2 broading and 29 bracketing iterations
## to a tolerance of 1.490116e-08
## 
## Progress is as follows
##          xl       xr            fl           fr
## 1  2.000000 3.000000  5.578119e-01 1.048825e+00
## 2  1.500000 3.500000  1.823349e-01 1.222566e+00
## 3  0.500000 4.500000 -1.299678e+00 1.492968e+00
## 4  0.500000 2.500000 -1.299678e+00 8.342057e-01
## 5  0.500000 1.500000 -1.299678e+00 1.823349e-01
## 6  1.000000 1.500000 -3.678794e-01 1.823349e-01
## 7  1.250000 1.500000 -6.336125e-02 1.823349e-01
## 8  1.250000 1.375000 -6.336125e-02 6.561414e-02
## 9  1.250000 1.312500 -6.336125e-02 2.787367e-03
## 10 1.281250 1.312500 -2.985381e-02 2.787367e-03
## 11 1.296875 1.312500 -1.342726e-02 2.787367e-03
## 12 1.304688 1.312500 -5.293741e-03 2.787367e-03
## 13 1.308594 1.312500 -1.246670e-03 2.787367e-03
## 14 1.308594 1.310547 -1.246670e-03 7.719731e-04
## 15 1.309570 1.310547 -2.369419e-04 7.719731e-04
## 16 1.309570 1.310059 -2.369419e-04 2.676172e-04
## 17 1.309570 1.309814 -2.369419e-04 1.536305e-05
## 18 1.309692 1.309814 -1.107831e-04 1.536305e-05
## 19 1.309753 1.309814 -4.770843e-05 1.536305e-05
## 20 1.309784 1.309814 -1.617229e-05 1.536305e-05
## 21 1.309799 1.309814 -4.045236e-07 1.536305e-05
## 22 1.309799 1.309807 -4.045236e-07 7.479287e-06
## 23 1.309799 1.309803 -4.045236e-07 3.537388e-06
## 24 1.309799 1.309801 -4.045236e-07 1.566434e-06
## 25 1.309799 1.309800 -4.045236e-07 5.809554e-07
## 26 1.309799 1.309800 -4.045236e-07 8.821597e-08
## 27 1.309799 1.309800 -1.581538e-07 8.821597e-08
## 28 1.309800 1.309800 -3.496891e-08 8.821597e-08
## 29 1.309800 1.309800 -3.496891e-08 2.662353e-08
## 30 1.309800 1.309800 -4.172689e-09 2.662353e-08
## 31 1.309800 1.309800 -4.172689e-09 1.122542e-08

The three plot commands below, arranged horizontally, give a nice visualization of progress.

par(mfrow=c(1,3))
plot(fr, main="bisection progress", ylab="f(x)", xlab="x")
plot(fr, main="bisection progress [zoom]", after=20, ylab="f(x)", xlab="x")
plot(fr, main="bisection progress [zoom]", after=25, ylab="f(x)", xlab="x")

  1. Use your modified function find the root of \[ h(x) = (x - 1)^3 - 2x^2 + 10 - \sin(x), \] starting with \(x_l = 1\) and \(x_r = 2\).

First define the new h function.

h <- function(x) (x-1)^3 - 2*x^2 - sin(x)

Here are the results, skipping summary.

hr <- bisection(h, 1, 2)
hr
## Root of:
## function(x) (x-1)^3 - 2*x^2 - sin(x)
## <bytecode: 0x4d9a2c0>
## in (1, 2): 4.307962
## found after 3 broading and 30 bracketing iterations
## to a tolerance of 1.490116e-08

The three plot commands below, arranged horizontally, give a nice visualization of progress.

par(mfrow=c(1,3))
plot(hr, main="bisection progress", ylab="f(x)", xlab="x")
plot(hr, main="bisection progress [zoom]", after=20, ylab="f(x)", xlab="x")
plot(hr, main="bisection progress [zoom]", after=25, ylab="f(x)", xlab="x")

Problem 4: R scripts from the Unix prompt (25 pts)

The two interfaces I was thinking of are Rscript and R CMD BATCH. although R -e offers a more limited yet related capability. The main difference is that BATCH is somewhat more customizable in terms of how it handles input, writes output, etc. To keep it simple here, I will demonstrate Rscript.

I made a file called rendpurl.R which has the following two commands.

lines <- readLines("rendpurl.R")
## Warning in readLines("rendpurl.R"): incomplete final line found on 'rendpurl.R'
cat(paste(lines, collapse="\n"))
## rmarkdown::render("hw2_sol.Rmd")
## knitr::purl("hw2_sol.Rmd")

My build.sh script has the following lines.

lines <- readLines("build.sh")
cat(paste(lines, collapse="\n"))
## #!/bin/bash
## Rscript rendpurl.R

Since the rendpurl.R script works on this very Rmarkdown document, running it here would create an infinite loop, so lets not do that. Simply download the files if you want to try them.