# R in Action (2nd ed): Chapter 20
# Advanced R programming
# requires packages ggplot2, reshape2, foreach, doParallel
# install.packages(c("ggplot2", "reshap2e", "foreach", "doParallel"))
# Atomic vectors
passed <- c(TRUE, TRUE, FALSE, TRUE)
ages <- c(15, 18, 25, 14, 19)
cmplxNums <- c(1+2i, 0+1i, 39+3i, 12+2i)
names <- c("Bob", "Ted", "Carol", "Alice")
# Matrices
x <- c(1,2,3,4,5,6,7,8)
class(x)
## [1] "numeric"
print(x)
## [1] 1 2 3 4 5 6 7 8
attr(x, "dim") <- c(2,4)
print(x)
## [,1] [,2] [,3] [,4]
## [1,] 1 3 5 7
## [2,] 2 4 6 8
class(x)
## [1] "matrix"
attributes(x)
## $dim
## [1] 2 4
attr(x, "dimnames") <- list(c("A1", "A2"),
c("B1", "B2", "B3", "B4"))
print(x)
## B1 B2 B3 B4
## A1 1 3 5 7
## A2 2 4 6 8
attr(x, "dim") <- NULL
class(x)
## [1] "numeric"
print(x)
## [1] 1 2 3 4 5 6 7 8
# Generic vectors (lists)
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
unclass(iris)
## $Sepal.Length
## [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4
## [18] 5.1 5.7 5.1 5.4 5.1 4.6 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5
## [35] 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8 5.1 4.6 5.3 5.0 7.0
## [52] 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8
## [69] 6.2 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4
## [86] 6.0 6.7 6.3 5.6 5.5 5.5 6.1 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8
## [103] 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8 6.4 6.5 7.7 7.7
## [120] 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7
## [137] 6.3 6.4 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9
##
## $Sepal.Width
## [1] 3.5 3.0 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 3.7 3.4 3.0 3.0 4.0 4.4 3.9
## [18] 3.5 3.8 3.8 3.4 3.7 3.6 3.3 3.4 3.0 3.4 3.5 3.4 3.2 3.1 3.4 4.1 4.2
## [35] 3.1 3.2 3.5 3.6 3.0 3.4 3.5 2.3 3.2 3.5 3.8 3.0 3.8 3.2 3.7 3.3 3.2
## [52] 3.2 3.1 2.3 2.8 2.8 3.3 2.4 2.9 2.7 2.0 3.0 2.2 2.9 2.9 3.1 3.0 2.7
## [69] 2.2 2.5 3.2 2.8 2.5 2.8 2.9 3.0 2.8 3.0 2.9 2.6 2.4 2.4 2.7 2.7 3.0
## [86] 3.4 3.1 2.3 3.0 2.5 2.6 3.0 2.6 2.3 2.7 3.0 2.9 2.9 2.5 2.8 3.3 2.7
## [103] 3.0 2.9 3.0 3.0 2.5 2.9 2.5 3.6 3.2 2.7 3.0 2.5 2.8 3.2 3.0 3.8 2.6
## [120] 2.2 3.2 2.8 2.8 2.7 3.3 3.2 2.8 3.0 2.8 3.0 2.8 3.8 2.8 2.8 2.6 3.0
## [137] 3.4 3.1 3.0 3.1 3.1 3.1 2.7 3.2 3.3 3.0 2.5 3.0 3.4 3.0
##
## $Petal.Length
## [1] 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 1.5 1.6 1.4 1.1 1.2 1.5 1.3
## [18] 1.4 1.7 1.5 1.7 1.5 1.0 1.7 1.9 1.6 1.6 1.5 1.4 1.6 1.6 1.5 1.5 1.4
## [35] 1.5 1.2 1.3 1.4 1.3 1.5 1.3 1.3 1.3 1.6 1.9 1.4 1.6 1.4 1.5 1.4 4.7
## [52] 4.5 4.9 4.0 4.6 4.5 4.7 3.3 4.6 3.9 3.5 4.2 4.0 4.7 3.6 4.4 4.5 4.1
## [69] 4.5 3.9 4.8 4.0 4.9 4.7 4.3 4.4 4.8 5.0 4.5 3.5 3.8 3.7 3.9 5.1 4.5
## [86] 4.5 4.7 4.4 4.1 4.0 4.4 4.6 4.0 3.3 4.2 4.2 4.2 4.3 3.0 4.1 6.0 5.1
## [103] 5.9 5.6 5.8 6.6 4.5 6.3 5.8 6.1 5.1 5.3 5.5 5.0 5.1 5.3 5.5 6.7 6.9
## [120] 5.0 5.7 4.9 6.7 4.9 5.7 6.0 4.8 4.9 5.6 5.8 6.1 6.4 5.6 5.1 5.6 6.1
## [137] 5.6 5.5 4.8 5.4 5.6 5.1 5.1 5.9 5.7 5.2 5.0 5.2 5.4 5.1
##
## $Petal.Width
## [1] 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 0.2 0.2 0.1 0.1 0.2 0.4 0.4
## [18] 0.3 0.3 0.3 0.2 0.4 0.2 0.5 0.2 0.2 0.4 0.2 0.2 0.2 0.2 0.4 0.1 0.2
## [35] 0.2 0.2 0.2 0.1 0.2 0.2 0.3 0.3 0.2 0.6 0.4 0.3 0.2 0.2 0.2 0.2 1.4
## [52] 1.5 1.5 1.3 1.5 1.3 1.6 1.0 1.3 1.4 1.0 1.5 1.0 1.4 1.3 1.4 1.5 1.0
## [69] 1.5 1.1 1.8 1.3 1.5 1.2 1.3 1.4 1.4 1.7 1.5 1.0 1.1 1.0 1.2 1.6 1.5
## [86] 1.6 1.5 1.3 1.3 1.3 1.2 1.4 1.2 1.0 1.3 1.2 1.3 1.3 1.1 1.3 2.5 1.9
## [103] 2.1 1.8 2.2 2.1 1.7 1.8 1.8 2.5 2.0 1.9 2.1 2.0 2.4 2.3 1.8 2.2 2.3
## [120] 1.5 2.3 2.0 2.0 1.8 2.1 1.8 1.8 1.8 2.1 1.6 1.9 2.0 2.2 1.5 1.4 2.3
## [137] 2.4 1.8 1.8 2.1 2.4 2.3 1.9 2.3 2.5 2.3 1.9 2.0 2.3 1.8
##
## $Species
## [1] setosa setosa setosa setosa setosa setosa
## [7] setosa setosa setosa setosa setosa setosa
## [13] setosa setosa setosa setosa setosa setosa
## [19] setosa setosa setosa setosa setosa setosa
## [25] setosa setosa setosa setosa setosa setosa
## [31] setosa setosa setosa setosa setosa setosa
## [37] setosa setosa setosa setosa setosa setosa
## [43] setosa setosa setosa setosa setosa setosa
## [49] setosa setosa versicolor versicolor versicolor versicolor
## [55] versicolor versicolor versicolor versicolor versicolor versicolor
## [61] versicolor versicolor versicolor versicolor versicolor versicolor
## [67] versicolor versicolor versicolor versicolor versicolor versicolor
## [73] versicolor versicolor versicolor versicolor versicolor versicolor
## [79] versicolor versicolor versicolor versicolor versicolor versicolor
## [85] versicolor versicolor versicolor versicolor versicolor versicolor
## [91] versicolor versicolor versicolor versicolor versicolor versicolor
## [97] versicolor versicolor versicolor versicolor virginica virginica
## [103] virginica virginica virginica virginica virginica virginica
## [109] virginica virginica virginica virginica virginica virginica
## [115] virginica virginica virginica virginica virginica virginica
## [121] virginica virginica virginica virginica virginica virginica
## [127] virginica virginica virginica virginica virginica virginica
## [133] virginica virginica virginica virginica virginica virginica
## [139] virginica virginica virginica virginica virginica virginica
## [145] virginica virginica virginica virginica virginica virginica
## Levels: setosa versicolor virginica
##
## attr(,"row.names")
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [18] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
## [35] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
## [52] 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
## [69] 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
## [86] 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
## [103] 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
## [120] 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
## [137] 137 138 139 140 141 142 143 144 145 146 147 148 149 150
attributes(iris)
## $names
## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## [5] "Species"
##
## $class
## [1] "data.frame"
##
## $row.names
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [18] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
## [35] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
## [52] 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
## [69] 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
## [86] 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
## [103] 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
## [120] 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
## [137] 137 138 139 140 141 142 143 144 145 146 147 148 149 150
set.seed(1234)
fit <- kmeans(iris[1:4], 3)
names(fit)
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
unclass(fit)
## $cluster
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [71] 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 3 3 3
## [106] 3 2 3 3 3 3 3 3 2 2 3 3 3 3 2 3 2 3 2 3 3 2 2 3 3 3 3 3 2 3 3 3 3 2 3
## [141] 3 3 2 3 3 3 2 3 3 2
##
## $centers
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053
##
## $totss
## [1] 681.3706
##
## $withinss
## [1] 15.15100 39.82097 23.87947
##
## $tot.withinss
## [1] 78.85144
##
## $betweenss
## [1] 602.5192
##
## $size
## [1] 50 62 38
##
## $iter
## [1] 2
##
## $ifault
## [1] 0
sapply(fit, class)
## cluster centers totss withinss tot.withinss
## "integer" "matrix" "numeric" "numeric" "numeric"
## betweenss size iter ifault
## "numeric" "integer" "integer" "integer"
# Indexing atomic vectors
x <- c(20, 30, 40)
x[3]
## [1] 40
x[c(2,3)]
## [1] 30 40
x <- c(A=20, B=30, C=40)
x[c(2,3)]
## B C
## 30 40
x[c("B", "C")]
## B C
## 30 40
# Indexing lists
fit[c(2,7)]
## $centers
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053
##
## $size
## [1] 50 62 38
fit[2]
## $centers
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053
fit[[2]]
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053
fit$centers
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053
fit[[2]][1,]
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.006 3.428 1.462 0.246
# fit$centers$Petal.Width # should give an error
# Listing 20.1 - Plotting the centroides from a k-mean cluster analysis
fit <- kmeans(iris[1:4], 3)
means <- fit$centers
library(reshape2)
dfm <- melt(means)
names(dfm) <- c("Cluster", "Measurement", "Centimeters")
dfm$Cluster <- factor(dfm$Cluster)
head(dfm)
## Cluster Measurement Centimeters
## 1 1 Sepal.Length 5.006000
## 2 2 Sepal.Length 6.850000
## 3 3 Sepal.Length 5.901613
## 4 1 Sepal.Width 3.428000
## 5 2 Sepal.Width 3.073684
## 6 3 Sepal.Width 2.748387
library(ggplot2)
ggplot(data=dfm,
aes(x=Measurement, y=Centimeters, group=Cluster)) +
geom_point(size=3, aes(shape=Cluster, color=Cluster)) +
geom_line(size=1, aes(color=Cluster)) +
ggtitle("Profiles for Iris Clusters")
# for loops
for(i in 1:5) print(1:i)
## [1] 1
## [1] 1 2
## [1] 1 2 3
## [1] 1 2 3 4
## [1] 1 2 3 4 5
for(i in 5:1)print(1:i)
## [1] 1 2 3 4 5
## [1] 1 2 3 4
## [1] 1 2 3
## [1] 1 2
## [1] 1
# ifelse
pvalues <- c(.0867, .0018, .0054, .1572, .0183, .5386)
results <- ifelse(pvalues <.05, "Significant", "Not Significant")
results
## [1] "Not Significant" "Significant" "Significant" "Not Significant"
## [5] "Significant" "Not Significant"
pvalues <- c(.0867, .0018, .0054, .1572, .0183, .5386)
results <- vector(mode="character", length=length(pvalues))
for(i in 1:length(pvalues)){
if (pvalues[i] < .05) results[i] <- "Significant"
else results[i] <- "Not Significant"
}
results
## [1] "Not Significant" "Significant" "Significant" "Not Significant"
## [5] "Significant" "Not Significant"
# Creating functions
f <- function(x, y, z=1){
result <- x + (2*y) + (3*z)
return(result)
}
f(2,3,4)
## [1] 20
f(2,3)
## [1] 11
f(x=2, y=3)
## [1] 11
f(z=4, y=2, 3)
## [1] 19
args(f)
## function (x, y, z = 1)
## NULL
# object scope
x <- 2
y <- 3
z <- 4
f <- function(w){
z <- 2
x <- w*y*z
return(x)
}
f(x)
## [1] 12
x
## [1] 2
y
## [1] 3
z
## [1] 4
# Working with environments
x <- 5
myenv <- new.env()
assign("x", "Homer", env=myenv)
ls()
## [1] "ages" "cmplxNums" "dfm" "f" "fit"
## [6] "i" "means" "myenv" "names" "passed"
## [11] "pvalues" "results" "x" "y" "z"
ls(myenv)
## [1] "x"
x
## [1] 5
get("x", env=myenv)
## [1] "Homer"
myenv <- new.env()
myenv$x <- "Homer"
myenv$x
## [1] "Homer"
parent.env(myenv)
## <environment: R_GlobalEnv>
# function closures
trim <- function(p){
trimit <- function(x){
n <- length(x)
lo <- floor(n*p) + 1
hi <- n + 1 - lo
x <- sort.int(x, partial = unique(c(lo, hi)))[lo:hi]
}
trimit
}
x <- 1:10
trim10pct <- trim(.1)
y <- trim10pct(x)
y
## [1] 2 3 4 5 6 7 8 9
trim20pct <- trim(.2)
y <- trim20pct(x)
y
## [1] 3 4 5 6 7 8
ls(environment(trim10pct))
## [1] "p" "trimit"
get("p", env=environment(trim10pct))
## [1] 0.1
makeFunction <- function(k){
f <- function(x){
print(x + k)
}
}
g <- makeFunction(10)
g (4)
## [1] 14
k <- 2
g (5)
## [1] 15
ls(environment(g))
## [1] "f" "k"
environment(g)$k
## [1] 10
# Generic functions
summary(women)
## height weight
## Min. :58.0 Min. :115.0
## 1st Qu.:61.5 1st Qu.:124.5
## Median :65.0 Median :135.0
## Mean :65.0 Mean :136.7
## 3rd Qu.:68.5 3rd Qu.:148.0
## Max. :72.0 Max. :164.0
fit <- lm(weight ~ height, data=women)
summary(fit)
##
## Call:
## lm(formula = weight ~ height, data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7333 -1.1333 -0.3833 0.7417 3.1167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
## height 3.45000 0.09114 37.85 1.09e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
## F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14
class(women)
## [1] "data.frame"
class(fit)
## [1] "lm"
methods(summary)
## [1] summary.aov summary.aovlist*
## [3] summary.aspell* summary.check_packages_in_dir*
## [5] summary.connection summary.data.frame
## [7] summary.Date summary.default
## [9] summary.ecdf* summary.factor
## [11] summary.ggplot* summary.glm
## [13] summary.infl* summary.lm
## [15] summary.loess* summary.manova
## [17] summary.matrix summary.mlm*
## [19] summary.nls* summary.packageStatus*
## [21] summary.PDF_Dictionary* summary.PDF_Stream*
## [23] summary.POSIXct summary.POSIXlt
## [25] summary.ppr* summary.prcomp*
## [27] summary.princomp* summary.proc_time
## [29] summary.srcfile summary.srcref
## [31] summary.stepfun summary.stl*
## [33] summary.table summary.tukeysmooth*
## [35] summary.warnings
## see '?methods' for accessing help and source code
# Listing 20.2 - An example of a generic function
mymethod <- function(x, ...) UseMethod("mymethod")
mymethod.a <- function(x) print("Using A")
mymethod.b <- function(x) print("Using B")
mymethod.default <- function(x) print("Using Default")
x <- 1:5
y <- 6:10
z <- 10:15
class(x) <- "a"
class(y) <- "b"
mymethod(x)
## [1] "Using A"
mymethod(y)
## [1] "Using B"
mymethod(z)
## [1] "Using Default"
class(z) <- c("a", "b")
mymethod(z)
## [1] "Using A"
class(z) <- c("c", "a", "b")
mymethod(z)
## [1] "Using A"
# Vectorization and efficient code
set.seed(1234)
mymatrix <- matrix(rnorm(10000000), ncol=10)
accum <- function(x){
sums <- numeric(ncol(x))
for (i in 1:ncol(x)){
for(j in 1:nrow(x)){
sums[i] <- sums[i] + x[j,i]
}
}
}
system.time(accum(mymatrix)) # using loops
## user system elapsed
## 1.407 0.010 1.432
system.time(colSums(mymatrix)) # using vectorization
## user system elapsed
## 0.014 0.001 0.014
# Correctly size objects
set.seed(1234)
k <- 100000
x <- rnorm(k)
y <- 0
system.time(for (i in 1:length(x)) y[i] <- x[i]^2)
## user system elapsed
## 0.043 0.007 0.051
y <- numeric(k)
system.time(for (i in 1:k) y[i] <- x[i]^2)
## user system elapsed
## 0.014 0.000 0.014
y <- numeric(k)
system.time(y <- x^2)
## user system elapsed
## 0.001 0.000 0.000
# Listing 20.3 - Parallelization with foreach and doParallel
library(foreach)
library(doParallel)
## Loading required package: iterators
## Loading required package: parallel
registerDoParallel(cores=4)
eig <- function(n, p){
x <- matrix(rnorm(100000), ncol=100)
r <- cor(x)
eigen(r)$values
}
n <- 1000000
p <- 100
k <- 500
system.time(
x <- foreach(i=1:k, .combine=rbind) %do% eig(n, p)
)
## user system elapsed
## 10.049 0.516 10.839
system.time(
x <- foreach(i=1:k, .combine=rbind) %dopar% eig(n, p)
)
## user system elapsed
## 15.102 1.048 5.082
# Listing 20.4 - A sample debugging session
args(mad)
## function (x, center = median(x), constant = 1.4826, na.rm = FALSE,
## low = FALSE, high = FALSE)
## NULL
debug(mad)
mad(1:10)
## debugging in: mad(1:10)
## debug: {
## if (na.rm)
## x <- x[!is.na(x)]
## n <- length(x)
## constant * if ((low || high) && n%%2 == 0) {
## if (low && high)
## stop("'low' and 'high' cannot be both TRUE")
## n2 <- n%/%2 + as.integer(high)
## sort(abs(x - center), partial = n2)[n2]
## }
## else median(abs(x - center))
## }
## debug: if (na.rm) x <- x[!is.na(x)]
## debug: n <- length(x)
## debug: constant * if ((low || high) && n%%2 == 0) {
## if (low && high)
## stop("'low' and 'high' cannot be both TRUE")
## n2 <- n%/%2 + as.integer(high)
## sort(abs(x - center), partial = n2)[n2]
## } else median(abs(x - center))
## debug: median(abs(x - center))
## exiting from: mad(1:10)
## [1] 3.7065
# enters debugging mode
# Q to quit - see text
undebug(mad)
# Listing 20.5 - Sample debugging session with recover()
f <- function(x, y){
z <- x + y
g(z)
}
g <- function(x){
z <- round(x)
h(z)
}
h <- function(x){
set.seed(1234)
z <- rnorm(x)
print(z)
}
options(error=recover)
f(2,3)
## [1] -1.2070657 0.2774292 1.0844412 -2.3456977 0.4291247
# f(2, -3) # enters debugging mode at this point