FoRt 2015/16: Lab: Debugging & Profiling
Fabian Scheipl
Aufgabe 1: Rainbow? . . . more like rainbug!
Nachfolgender Code erzeugt eine n_grid×n-Matrix x mit Auswertungen von n zufälligen Funktionen xi (t)
auf je n_grid Gitterpunkten zwischen 0 und 1 (t_grid):
library(rainbow)
set.seed(121212)
n <- 80
n_grid <- 100
t_grid <- seq(0, 1, l = n_grid)
x <- replicate(n,
dbeta(t_grid, shape1 = runif(1, 3, 4), shape2 = runif(1, 2, 7))
rt(n_grid, df = 5)/10
)
+
Das linke Panel in untenstehender Grafik zeigt die Realisierungen der einzelnen xi (t), i = 1, . . . , n.
Das Paket rainbow [pdf] stellt Plotfunktionen für solche funktionalen Daten zur Verfügung, unter anderem
eine Art Boxplot für funktionale Daten mit fboxplot(). Diese Funktion produziert hier zwar den korrekten
Plot (rechtes Panel in der folgenden Grafik), aber auch eine Fehlermeldung:
# fds creates "functional data"-objects, see ?fds and ?fboxplot
x_fds <- fds(x = t_grid, y = x)
layout(t(1:2))
matplot(t_grid, x, lty = 1, col = rgb(0, 0, 0, .2), type = "l", lwd = 1.5)
fboxplot(x_fds)
3
2
1
−2
−1
0
x
−2
−1
0
x
1
2
3
## Error in as.graphicsAnnot(legend): argument "legend" is missing, with no default
0.0
0.2
0.4
0.6
0.8
1.0
0.0
t_grid
0.2
0.4
0.6
0.8
1.0
t_grid
a) Lokalisieren Sie wo im Code der Fehler auftritt, beschreiben Sie die Ursache des Fehlers und ändern Sie
den obigen Code so, dass er nicht mehr auftritt.
b) Wie könnte der Code im rainbow-Paket verbessert werden um den obigen Fehler zu vermeiden? (hier
kein Code gefragt, nur Lösungsideen. . . )
1
Lösung:
a)
Lokalisierung: erst mal mit traceback():
traceback()
##
##
##
##
##
##
##
4: as.graphicsAnnot(legend)
3: legend(legendpos, c(colnames(data$y)[outlier]), col = rainbow(n),
lty = 1, ncol = ncol, ...)
2: fbag(data, factor, xlab = xlab, ylab = ylab, plotlegend = plotlegend,
legendpos = legendpos, ncol = ncol, projmethod = projmethod,
...)
1: fboxplot(x_fds)
Also: fboxplot() ruft fbag() ruft legend() ruft as.graphicsAnnot() auf, aber as.graphicsAnnot()
bekommt von legend() offensichtlich kein legend-Argument übergeben – siehe Fehlermeldung. Das heisst
der Fehler passiert wohl schon in legend(), oder sogar noch früher.
Mit options(error=recover, deparse.max.lines=5) können wir direkt in die Funktionsauswertung springen und uns mit dem interaktiven Debugger dort umsehen:
options(error=recover, deparse.max.lines=5)
fboxplot(x_fds)
##
##
##
##
##
##
##
##
##
Error in as.graphicsAnnot(legend) :
argument "legend" is missing, with no default
Enter a frame number, or 0 to exit
1:
2:
3:
4:
fboxplot(x_fds)
fbag(data, factor, xlab = xlab, ylab = ylab, plotlegend = plotlegend, legendpos = legendpos, ncol
legend(legendpos, c(colnames(data$y)[outlier]), col = rainbow(n), lty = 1, ncol = ncol, ...)
as.graphicsAnnot(legend)
Wir springen in den Aufruf von legend() um nachzusehen ob/wie dort die Variable legend definiert ist die
bei as.graphicsAnnot(legend) zu fehlen scheint:
## Selection: 3
## Called from: legend(legendpos, c(colnames(data$y)[outlier]), col = rainbow(n),
##
lty = 1, ncol = ncol, ...)
Jetzt sind wir im interaktiven Debugger und können uns umschauen:
## Browse[1]> ls.str()
## adj : num [1:2] 0 0.5
## angle : num 45
[...]
## inset : num 0
## legend : <missing>
## lty : num 1
2
Aha, also die legend-Variable ist hier gar nicht definiert: <missing>. Zeit sich die Hilfe für legend()
anzusehen. Dort lesen wir dass legend das dritte Argument von legend() ist, bzw. das zweite Argument
ohne default-Wert. Der Aufruf von legend() hier ist
legend(legendpos, c(colnames(data$y)[outlier]), col = rainbow(n), lty = 1,
ncol = ncol, ...)
das heisst nach den Regeln für das argument matching ist c(colnames(data$y)[outlier]) hier der Wert, der
für das legend-Argument übergeben wird – sieht so aus als würde fboxplot versuchen die farbig markierten
Ausreisser im Datensatz mit den Spaltennamen der übergebenen Daten zu labeln. Springen wir also in
die Funktion fbag() die legend() aufruft und schauen nach was bei der Übergabe des legend-Arguments
schiefläuft:
fboxplot(x_fds)
##
##
##
##
##
##
##
##
##
##
##
##
##
Error in as.graphicsAnnot(legend) :
argument "legend" is missing, with no default
Enter a frame number, or 0 to exit
1: fboxplot(x_fds)
2: fbag(data, factor, xlab = xlab, ylab = ylab, plotlegend = plotlegend, legendpos = legendpos, ncol
3: legend(legendpos, c(colnames(data$y)[outlier]), col = rainbow(n), lty = 1, ncol = ncol, ...)
4: as.graphicsAnnot(legend)
Selection: 2
Called from: fbag(data, factor, xlab = xlab, ylab = ylab, plotlegend = plotlegend,
legendpos = legendpos, ncol = ncol, projmethod = projmethod,
...)
Was wird hier mit c(colnames(data$y)[outlier]) übergeben?
## Browse[1]> c(colnames(data$y)[outlier])
## NULL
AHA! Das ist also das Problem und deswegen ist legend im Aufruf von legend() einfach <missing>:
Das Argument das übergeben wird hat den Wert NULL. Warum ist es NULL? Liegt es an outlier oder an
colnames(data$y)?
##
##
##
##
##
##
##
##
##
##
##
Browse[1]> outlier
[1] 40 53 68
Browse[1]> colnames(data$y)
NULL
Browse[1]> str(data)
List of 4
$ x
: num [1:100] 0 0.0101 0.0202 0.0303 0.0404 ...
$ y
: num [1:100, 1:80] 0.3319 -0.1903 0.2472 0.0283 0.0171 ...
$ xname: chr "t_grid"
$ yname: chr "x"
- attr(*, "class")= chr "fds"
Das ist also das Problem: Das von uns an fboxplot übergebene data-Argument (=x_fds) hat keine
Spaltennamen für die Funktionsauswertungen (die in der $y-Komponente des fds-Objekts abgespeichert sind)
aber fboxplot geht ohne Überprüfung davon aus dass diese Spaltennamen vorhanden sind. . .
3
Um den Fehler zu vermeiden gibt es 2 Möglichkeiten: Wir rufen fboxplot mit plotlegend = FALSE auf um
gar nicht erst zu versuchen die Legende zu zeichnen oder wir definieren Spaltennamen für x und erzeugen
x_fds nochmal neu mit diesem x damit die Legende gezeichnet werden kann:
layout(t(1:2))
fboxplot(x_fds, plotlegend = FALSE)
3
68
1
2
40
53
−2
−1
0
x
−2
−1
0
x
1
2
3
colnames(x) <- 1 : ncol(x)
x_fds <- fds(x = t_grid, y = x)
fboxplot(x_fds)
0.0
0.2
0.4
0.6
0.8
1.0
0.0
t_grid
0.2
0.4
0.6
0.8
1.0
t_grid
b)
Das Problem ist wie so oft mangelhaftes input checking, hier von fboxplot(), und zwar in dem Sinn dass
der Code implizit voraussetzt dass die übergebenen Daten Spaltennamen haben. Korrekturmöglichkeiten:
• Überprüfe ob die Spaltennamen existieren bevor legend() aufgerufen wird, z.B. mit is.null(colnames(data$y)).
Wenn es keine Spaltennamen gibt werden Sie einfach als seq_len(ncol(data$y)) definiert und dann
wird legend() aufgerufen.
• Am Anfang der Funktion wird überprüft ob die Spaltennamen existieren. Wenn nicht wird plotlegend
auf FALSE gesetzt (evtl. mit einer warning) damit legend() gar nicht erst aufgerufen wird.
• Sauberste Lösung: fds() überprüft ob Spaltennamen vorhanden sind und legt sie an falls nicht.
Dann können alle Funktionen die fds-Objekte verarbeiten sich darauf verlassen dass die Spaltennamen
existieren.
Aufgabe 2: Mopsgeschwindigkeit!
Der Code in slow_sim.R implementiert eine (völlig sinnfreie) Simulationsstudie um die Verteilung der
geschätzten Regressionskoeffizienten β̂ in einem Modell y ∼ t(ncp = Xβ, df = 4) mit t(4)-verteilten Fehlern
und linearem Prädiktor Xβ zu bestimmen:
source("slow_sim.R")
#######################
set.seed <- 232323
observations <- 5000
4
covariates <- 10
testdata <- as.data.frame(matrix(rnorm(observations * covariates),
nrow = observations))
system.time(test <- simulate(reps = 100, seed = 20141028, data = testdata))
##
##
user
1.986
system elapsed
0.021
2.008
Die Simulation ist recht ineffizient programmiert. Lesen Sie sich zunächst mal den (bewusst unkommentierten)
Quellcode genau durch, lesen sie die Hilfe für alle dort verwendeten Funktionen die Sie nicht gut kennen.
Benutzen Sie falls nötig den Debugger um zu verstehen was die einzelnen Befehle genau tun.
a) Benutzen Sie die in der Vorlesung kennengelernten Profiling-Methoden um die Stellen zu identifizieren
an denen das Skript in slow_sim.R die meiste Zeit verbringt.
b) Modifizieren Sie den Code in slow_sim.R so, dass er mindestens 4x schneller läuft (ohne dass sich die
Ergebnisse ändern!)
Hinweis: Betrachten Sie zu a) sowohl wo in dem Code von slow_sim.R die meiste Zeit verbraucht wird als
auch welche eingebauten R-Funktionen dort aufgerufen werden und was diese tun und wie. Denken Sie auch
daran dass das line profiling nur für Code der aus einem Skript per source() geladen wurde vernünftig
funktioniert.
Lösung:
a)
Schauen wir zunächt mal nur auf unseren eigen Code, das heisst Rprof mit line.profiling und
summaryRprof() mit lines="show":
Rprof("slow_sim.log", line.profiling=TRUE)
test <- simulate(reps = 100, seed = 20141028, data = testdata)
Rprof(NULL)
(profile_slow_lines <- summaryRprof("slow_sim.log", lines = "show"))
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
$by.self
slow_sim.R#29
slow_sim.R#22
slow_sim.R#24
slow_sim.R#23
slow_sim.R#30
$by.total
slow_sim.R#10
slow_sim.R#18
slow_sim.R#29
slow_sim.R#17
slow_sim.R#22
slow_sim.R#24
slow_sim.R#23
slow_sim.R#30
self.time self.pct total.time total.pct
0.92
56.79
0.92
56.79
0.48
29.63
0.48
29.63
0.14
8.64
0.14
8.64
0.06
3.70
0.06
3.70
0.02
1.23
0.02
1.23
total.time total.pct self.time self.pct
1.62
100.00
0.00
0.00
0.94
58.02
0.00
0.00
0.92
56.79
0.92
56.79
0.68
41.98
0.00
0.00
0.48
29.63
0.48
29.63
0.14
8.64
0.14
8.64
0.06
3.70
0.06
3.70
0.02
1.23
0.02
1.23
5
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
$by.line
slow_sim.R#10
slow_sim.R#17
slow_sim.R#18
slow_sim.R#22
slow_sim.R#23
slow_sim.R#24
slow_sim.R#29
slow_sim.R#30
self.time self.pct total.time total.pct
0.00
0.00
1.62
100.00
0.00
0.00
0.68
41.98
0.00
0.00
0.94
58.02
0.48
29.63
0.48
29.63
0.06
3.70
0.06
3.70
0.14
8.64
0.14
8.64
0.92
56.79
0.92
56.79
0.02
1.23
0.02
1.23
$sample.interval
[1] 0.02
$sampling.time
[1] 1.62
Aha – interessant ist hier die $by.line-Tabelle. Sie zeigt dass unser Simulationscode etwa 30% der Zeit
damit verbringt immer und immmer wieder die selbe Designmatrix X zu erzeugen (Zeile 22 ist design <model.matrix(~., data = data)). Etwa 9% der Zeit gehen damit drauf in jeder Replikation t-verteilte
Fehler zu erzeugen (Zeile 24 ist data$y <- expected + rt(nrow(data), df = df)), da können wir nicht
viel machen.
Wenig überraschend ist zunächst mal auch dass mit Abstand die meiste Zeit (ca. 57%) dafür verbraucht wird
in jeder Replikation das lineare Modell zu schätzen – Zeile 29 ist model <- lm(y ~., data=data). Hier
lohnt es sich also genauer nachzusehen was dort so lange dauert. Dafür ist summaryRprof() mit lines =
"hide" ganz gut geeignet:
profile_slow <- summaryRprof("slow_sim.log", lines = "hide")
profile_slow$by.total
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
"cbind"
"simulate"
"simulate_once"
"estimate_coef"
"lm"
".External2"
"simulate_response"
"model.matrix"
"model.matrix.default"
"model.frame.default"
"na.omit.data.frame"
"na.omit"
"<Anonymous>"
"eval"
"model.frame"
"lm.fit"
"model.response"
"[.data.frame"
"rt"
"["
"as.character"
total.time total.pct self.time self.pct
1.62
100.00
0.00
0.00
1.62
100.00
0.00
0.00
1.62
100.00
0.00
0.00
0.94
58.02
0.00
0.00
0.92
56.79
0.00
0.00
0.84
51.85
0.38
23.46
0.68
41.98
0.00
0.00
0.66
40.74
0.02
1.23
0.64
39.51
0.00
0.00
0.56
34.57
0.02
1.23
0.46
28.40
0.26
16.05
0.46
28.40
0.00
0.00
0.34
20.99
0.00
0.00
0.34
20.99
0.00
0.00
0.22
13.58
0.00
0.00
0.18
11.11
0.18
11.11
0.16
9.88
0.00
0.00
0.14
8.64
0.14
8.64
0.14
8.64
0.14
8.64
0.14
8.64
0.00
0.00
0.12
7.41
0.12
7.41
6
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
"FUN"
"lapply"
"%in%"
"deparse"
"paste"
"sapply"
"match"
".deparseOpts"
"%*%"
"[["
"is.na"
"[[.data.frame"
".getXlevels"
"anyNA"
"as.vector"
"is.data.frame"
"makepredictcall.default"
"pmatch"
"unname"
"makepredictcall"
0.12
0.12
0.10
0.10
0.10
0.10
0.08
0.08
0.06
0.06
0.04
0.04
0.04
0.02
0.02
0.02
0.02
0.02
0.02
0.02
7.41
7.41
6.17
6.17
6.17
6.17
4.94
4.94
3.70
3.70
2.47
2.47
2.47
1.23
1.23
1.23
1.23
1.23
1.23
1.23
0.02
0.00
0.04
0.00
0.00
0.00
0.06
0.00
0.06
0.02
0.04
0.00
0.00
0.02
0.02
0.02
0.02
0.02
0.02
0.00
1.23
0.00
2.47
0.00
0.00
0.00
3.70
0.00
3.70
1.23
2.47
0.00
0.00
1.23
1.23
1.23
1.23
1.23
1.23
0.00
Hier sehen wir jetzt
1. unser Code verbringt die komplette Zeit in cbind, was daran liegt dass die Zeile coefs <- cbind(coefs,
simulate_once(data, true_coef, df)) quasi alle weiteren von uns definierten/benutzten Funktionen
aufruft die nennenswert Rechenzeit benötigen. Hier könnten wir evtl. durch pre-allocation von coefs
ein bisschen Zeit sparen und Duplikate vermeiden.
2. Interessanter ist aber dass wir etwa 35% der Zeit in model.frame.default verbringen (und etwa
41% der Zeit in model.matrix bzw. model.matrix.default, das haben wir aber schon weiter oben
festgestellt. . . ). model.frame wird von lm() aufgerufen und checkt im Datensatz nach NAs, ob alle
Variablen in der Formel auch wirklich im Datensatz vorhanden sind usw., alles Dinge die wir hier nicht
brauchen. Wir können stattdessen direkt lm.fit() aufrufen um das Modell zu schätzen.
b)
In faster_sim.R findet sich Code (hier aufgeführt) der die oben festgestellten Flaschenhäse & Ineffizienzen
beseitigt:
simulate_faster <- function(reps, seed, data, true_coef=0:ncol(data), df=4) {
set.seed(seed)
coefs <- matrix(0, nrow = length(true_coef), ncol = reps)
#change 1: pre-compute design X, expected (X*beta):
design <- model.matrix(~., data = data)
expected <- design %*% true_coef
for(rep in seq_len(reps)){
coefs[, rep] <- simulate_once_faster(expected, design, df)
}
return(structure(coefs,
seed = seed))
}
simulate_once_faster <- function(expected, design, df) {
7
response <- simulate_response_faster(expected, df)
estimate_coef_faster(response, design)
}
simulate_response_faster <- function(expected, df) {
expected + rt(length(expected), df = df)
}
estimate_coef_faster <- function(response, design) {
#change 2: lm.fit instead of lm to avoid unnecessary input checks
model <- lm.fit(y = response, x = design)
coef(model)
}
source("faster_sim.R")
all.equal(simulate(reps = 10, seed = 20141028, data = testdata),
simulate_faster(reps = 10, seed = 20141028, data = testdata))
## [1] TRUE
library(rbenchmark)
benchmark(slow = simulate(reps = 30, seed = 20141028, data = testdata),
faster = simulate_faster(reps = 30, seed = 20141028, data = testdata),
replications = 30, columns = c('test', 'elapsed', 'relative'))
##
test elapsed relative
## 2 faster
2.556
1.000
## 1
slow 21.382
8.365
Das hat sich also durchaus gelohnt. Schauen wir uns an wie jetzt die Rechenzeit verteilt ist:
Rprof("faster_sim.log", line.profiling = TRUE)
test <- simulate_faster(reps = 100, seed = 20141028, data = testdata)
Rprof(NULL)
profile_faster_lines <- summaryRprof("faster_sim.log", lines = "show")
profile_faster_lines$by.line
##
##
##
##
##
##
faster_sim.R#11
faster_sim.R#17
faster_sim.R#18
faster_sim.R#21
faster_sim.R#25
self.time self.pct total.time total.pct
0.00
0.00
0.72
100.00
0.00
0.00
0.30
41.67
0.00
0.00
0.42
58.33
0.30
41.67
0.30
41.67
0.42
58.33
0.42
58.33
Das macht jetzt auch mehr Sinn als vorher: Wir verbringen 58% der Zeit damit das lineare Modell zu schätzen
und 42% der Zeit damit t-verteilte Zufallsvariablen zu generieren. Alle anderen Operationen brauchen keine
nennenswerte Rechenzeit.
Mehr Möglichkeiten für bessere Performance:
1. Parallelisierung: Wie fast alle Simulationsstudien ist das hier natürlich embarrassingly parallel, das
heisst so etwas wie der Code weiter unten wäre eine Möglichkeit weitere Zeit zu sparen wenn man
mehrere Prozesse gleichzeitig laufen lassen kann:
8
simulate_parallel <- function(reps, seed, data, true_coef=0:ncol(data), df=4,
cluster=NULL, cores = NULL) {
set.seed(seed)
design <- model.matrix(~., data = data)
expected <- design %*% true_coef
library(parallel)
if(is.null(cores)){
cores <- detectCores() - 1
}
## parallel apply-function <papply> is platform-dependent:
if(.Platform$OS.type != "windows") {
# simply use forking for unix, mac:
papply <- function(X, FUN) {
mclapply(X = X, FUN = FUN, mc.cores = cores)
}
} else {
# use socket clusters on windows
if(is.null(cluster)) {
cluster <- makePSOCKcluster(rep("localhost", detectCores()-1))
on.exit(stopCluster(cl = cluster))
}
# make sure RNG is set to parallel mode:
RNGkind("L'Ecuyer-CMRG")
clusterSetRNGStream(cluster)
# load needed objects onto cluster:
clusterExport(cl = cluster, varlist = c("simulate_once_faster",
"simulate_response_faster",
"estimate_coef_faster", "expected",
"design", "df"))
papply <- function(X, FUN) {
parLapply(cl = cluster, X = X, fun = FUN)
}
}
coefs <- papply(X = seq_len(reps), FUN = function(i) {
unname(simulate_once_faster(expected, design, df))
})
coefs <- do.call(cbind, coefs)
return(structure(coefs,
seed = seed))
}
all.equal(simulate_faster(reps = 10, seed = 20141028, data = testdata),
simulate_parallel(reps = 10, seed = 20141028, data = testdata))
## [1] "Mean relative difference: 0.00526049"
benchmark(faster = simulate_faster(reps = 100, seed = 20141028, data = testdata),
parallel = simulate_parallel(reps = 100, seed = 20141028,
data = testdata),
replications = 9, columns = c('test', 'elapsed', 'relative'))
9
##
test elapsed relative
## 1
faster
2.411
3.52
## 2 parallel
0.685
1.00
Obacht: Parallelisierte Zufallsgeneratoren generieren andere Zufallszahlen!
Ergebnis hier unter Benutzung von 31 Prozessen.
2. Byte-compilation: bringt hier eher nichts, weil die relevanten Teile des Codes großteils sowieso schon
optimierte high-level Funktionen aus base und stats aufrufen. . .
3. Spezialisierte Pakete: RcppArmadillo hat eine Funktion fastLm() bzw. fastLmPure(), mal schauen
was die kann:
simulate_rcpp <- function(reps, seed, data, true_coef=0:ncol(data), df=4) {
set.seed(seed)
coefs <- matrix(0, nrow=length(true_coef), ncol=reps)
design <- model.matrix(~., data = data)
expected <- design %*% true_coef
for(rep in seq_len(reps)){
coefs[, rep] <- simulate_once_rcpp(expected, design, df)
}
return(structure(coefs,
seed = seed))
}
simulate_once_rcpp <- function(expected, design, df) {
response <- simulate_response_faster(expected, df)
estimate_coef_rcpp(response, design)
}
estimate_coef_rcpp <- function(response, design) {
stopifnot(require(RcppArmadillo))
model <- fastLmPure(y = response, X = design)
model$coefficients
}
all.equal(simulate_faster(reps = 10, seed = 20141028, data = testdata),
simulate_rcpp(reps = 10, seed = 20141028, data = testdata))
## Loading required package: RcppArmadillo
## [1] TRUE
benchmark(faster = simulate_faster(reps = 30, seed = 20141028, data = testdata),
rcpp = simulate_rcpp(reps = 30, seed = 20141028, data = testdata),
replications = 30, columns = c('test', 'elapsed', 'relative'))
##
test elapsed relative
## 1 faster
3.221
1.38
## 2
rcpp
2.334
1.00
Also nochmal bissle schneller durch modernste C++-Magie!
10
© Copyright 2024 ExpyDoc
