Department of Statistics, Virginia Tech

*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
```

*Nearly a second to “touch” three numbers isn’t super speedy.*

*Your tasks on this problem are the following.*

*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 `quote`

d 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
```

- Notice how the first evaluation takes a substantial amount of time, but subsequent evaluations are essentially instantaneous.

- The same thing does not happen with consecutive evaluations of the original
`swap`

function.

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

- Each evaluation consistently takes about 8/10 of a second.

*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 `quote`

d 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
```

- Very similar behavior to part a.: slow the first time, but essentially instantaneous afterwards.

*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
```

- So the behavior is quite similar to the original
`swap`

function, perhaps ever-so-slightly faster. - That’s one disadvantage: not nearly as fast as
`swap.eval`

.

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`

- Notice that
`v`

has changed*within*the function.

After:

`v[1:5]`

`## [1] 1 2 3 4 5`

- Back the way it was at the start.

- So the a change has not taken effect in the parent frame. That’s another disadvantage.

Consider the `swap.eval`

alternative:

```
swap.eval()
v[1:5]
```

`## [1] 2 1 3 4 5`

`v`

is indeed changed! That seems more useful.

*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,*

*reduce \(x_l\) and/or increase \(x_r\), and try again.**A reasonable choice is to double the width of the interval, i.e.,*

\[ \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.*

*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

- An
`if`

statement at the top of the`bisection`

fucntion checking if broadening is required, and if so implementing at most`100`

(by default) iterations of the algorithm described above via a new function called`broaden`

. That function has a similar`verb`

osity feature allowing progress to be printed to the screen in real time. - Changes to the S3
`print`

and`summary`

methods allow that progress to be visualized in text form after the fact. Actually, the`summary`

method is unchanged, but since it calls`print`

it benefits from changes therein to report broadening and bracketing iterations separately. - Changes to the
`plot`

S3 method similarly augment visualizations of progress. Broadening iterations are noted by a “+” mark above iteration numbers.

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

`source("bisection2.R")`

*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
```

- Observe how it separately reports on broadening and bracketing iterations.
- Notice that this is the same answer we got with earlier with a less pathological initialization.

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

- Besides the updated print, the output is pretty similar to our earlier version.

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")
```

- Observe how the first plot has indicies with a “+” above them, indicating broadening iterations.

*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
```

- Looks like we found a root of 4.31 after 3 and 30 bracketing iterations.

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")
```

- We can see that substantial broadening was required.
- And also notice how “close” we got to finding a root near \(x=0.5\), but the algorithm wasn’t fooled!.

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 hope your own solution was somewhat more verbose than the above. But I didn’t see any need to regurgitate the relevent help files here.

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")
```

- Ignore the
`##`

above, which comes from`Rmarkdown`

, not from the R script. `render`

creates the`.html`

file and`purl`

creates the`.R`

version.

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.