1 Testowanie czasu wykonywania kodu

Za pomocą polecenia system.time() możemy zmierzyć czas działania bloku instrukcji.

1.1 Zadanie - czas działania kodu

Sprawdzić działanie poniższych instrukcji (wszystkie robią to samo, jednak w różnym czasie). Wyjaśnić różnice w sposobie implementacji. Uwaga: kod chwilę się wykonuje.

system.time(
  {x=NULL
  for(i in 1:10^5) x=c(x,runif(1))
  }
)

system.time(
  {x=NULL
  for(i in 1:10^5) x[i]=runif(1)
  }
)

system.time(
  {x=numeric(10^5)
  for(i in 1:10^5) x[i]=runif(1)
  }
)

system.time(
  {x=NULL
  x=runif(10^5)
  }
)

2 Microbenchmark

Korzystając z pakietu microbenchmark możemy dostać wyniki z dokładnością do mikrosekundy.

Uwaga: pakiet trzeba zainstalować.

Przykład wykorzystania:

#---------biblioteki
library(microbenchmark)
library(ggplot2) #do graficznej prezentacji wyników

#instrukcje testowe takie jak poprzednio
pomiary=microbenchmark(
  skladanie    = {x=NULL; for(i in 1:1000) x=c(x,runif(1))},
  skladanie2   = {x=NULL; for(i in 1:1000) x[i]=runif(1)},
  inicjacja    = {x=numeric(1000); for(i in 1:1000) x[i]=runif(1)},
  wektoryzacja = {x=NULL; x=runif(1000)}
)

#------ wynik
pomiary

## Unit: microseconds
##          expr    min      lq     mean  median      uq     max neval
##     skladanie 4243.5 4848.00 5614.133 5072.70 5419.30 17789.6   100
##    skladanie2 2383.3 2527.65 3538.424 2780.90 3038.55 19569.4   100
##     inicjacja 2091.3 2263.15 3142.328 2461.25 2616.95 20015.6   100
##  wektoryzacja   24.0   26.85   33.094   28.25   30.65    95.1   100

#------ na obrazku
autoplot(pomiary)

Uwaga: każdy blok instrukcji wykonywany jest 100 razy dla celów testowych, w tabelce i na wykresie mamy wyniki dla tych 100 wykonań. Liczbę powtórzeń można zmienić (argument times w poleceniu microbenchmark).

2.1 Zadania - microbenchmark

Uwaga:

Wyniki w poniższych działaniach należy zilustrować za pomocą wykresów i/lub tabelek.
Otrzymane wyniki można porównać z wynikami kolegów (a to może być ciekawe).
Wyniki mogą zależeć od użytego urządzenia (komputer w pracowni, laptop na baterii, laptop z zasilacza itp.).

2.1.1 Zadanie microbenchmark - norma wektora

Zbudować i przetestować (porównać) za pomocą microbenchmarka kod wykonujący poniższe zadanie w dwóch wersjach: poprzez pętle i w wersji wektorowej.

Zadanie dla kodu:

przygotowanie (przed microbenchmarkiem): wylosować z rozkładu normalnego dwa wektory \(x\) i \(y\) o długości \(N=1000\).
właściwe zadanie: obliczenie normę euklidesową różnicy \(x-y\)

Chcemy porównać 2 wersje:

pętla for - startujemy od sumy równej zero i w pętli dodajemy kolejno kwadraty różnic elementów wektorów. Na koniec liczymy pierwiastek kwadratowy uzyskanej sumy
wektorowo - odejmujemy wektory, robimy “kwadrat” każdego elementu, sumujemy poleceniem sum i pierwiastek kwadratowy na koniec.

Przykładowy wynik

2.1.2 Zadanie microbenchmark - norma wektora - czas działania w zależności od liczby elementów

Ciąg dalszy poprzedniego zadania. Tym razem chcemy sprawdzić iloraz czasów działania wersji iteracyjnej i wektorowej (ich średnich) w zależności od ilości elementów w x (czy y, to jest ta sama długość).

Wyniki należy zilustrować w postaci wykresu zależności długość x vs. iloraz czasu działania (patrz przykład wyniku poniżej).

Zakres długości wektora x do przebadania to od n=1000 do n=20*1000 co 1000.

Może przydać się odpowiednia pętla for z losowaniem wektorów x i y o kolejnych długościach.

Ponieważ microbenchmark zwraca ramka danych jako wynik składającą się z expr oraz time, poszczególne średnie czasów możemy uzyskać poleceniem tapply:

tapply(pomiary_dlugosc$time, pomiary_dlugosc$expr, mean)

##   wektorowo iteracyjnie 
##        5798     3193457

Przykładowy wynik

Z obrazka widać, że iloraz stabilizuje się na wartości ok. 100, czyli można powiedzieć, że wersja wektorowa jest przyzwoicie 100x szybsza od wersji iteracyjnej

2.1.3 Zadanie microbenchmark - `for` vs. `sapply` vs. `vapply` vs. `wektorowo`

Porównamy czas działania pętli for oraz funkcji sapply i vapply w przypadku obliczania 10.000 wywołań wybranej funkcji. Porównany także wersję wektorową.

Do testów zbudujemy funkcję (taką, którą da się zastosować bezpośrednio do wektora), np. użyjemy: \[ f(x) = \sin(2x) + \frac{e^x}{x^2+1}. \] na przedziale \(x \in [0,10]\) podzielonym (seq) na 10.000 części (length.out=10000).

Wynikiem działania ma być wektor y tej samej długości zawierający wartości funkcji w podanych punktach przedziału.

Porównać za pomocą microbenchmarka 4 wersje obliczeń:

poprzez zastosowanie pętli for dla 10.000 elementów
poprzez zastosowanie funkcji vapply
poprzez zastosowanie funkcji sapply
poprzez bezpośrednie wywołanie f(x) dla wektora x

Przykładowy wynik

2.1.4 Zadanie microbenchmark - `runif`

Uwaga: jest to rozwinięcie jednego z wcześniejszych zadań.

Utworzymy wektor losowy x z rozkładu jednostajnego na przedziale [0,1] o długości n=1000 za pomocą polecenia runif. Następnie utworzymy wektor y o tej własności, że y[i]=0 gdy x[i]<0.5 oraz y[i]=1 w przeciwnym przypadku, \(i=1, \ldots, n\).

W każdym przypadku należy zainicjować “pusty = złożony z samych zer” wektor za pomocą polecenia y=numeric(n).

Za pomocą microbenchmarka porównać 4 wersje wykonania tego zadania:

z pętlą for oraz klasycznym, bezpośrednim if... else
z pętlą for oraz if - ale wstawiamy jedynki tylko tam gdzie trzeba, bo “zera” już są z polecenia numeric
za pomocą polecenia ifelse(...) zastosowanym do wektora
za pomocą “logicznej” konstrukcji y[x>=0.5]=1

Moje wyniki dla porównania

2.1.5 Zadanie microbenchmark - paproć Barnsleya lub zbiór Mandelbrota

Do zadania wykorzystamy własny kod generujący paproć Barnsleya lub zbiór Mandelbrota utworzony w Laboratorium 1.

Porównamy 2 wersje:

samo wygenerowanie punktów bez rysunku
wygenerowanie punktów plus rysunek w ggplot2 przypisany do zmiennej (tutaj rysunek najprostszy, samo geom_point wystarczy).

Uwaga:

zamiast wyświetlania obrazka należy przypisać go do zmiennej.
pewne elementy, np. funkcję testową, należy zrobić przed benchmarkiem
liczba iteracji: 10.000

Moje wyniki dla porównania

3 Testowanie składowych bloków - profiler

Do testowania składowych bloku kodu można wykorzystać profiler:

Rprof("profiler.out", interval=0.01)
N=2000
df= as.data.frame(matrix(runif(2*N*N),2*N,N))
tmp =summary(lm(V1~.,data=df))
Rprof(NULL)
summaryRprof("profiler.out")$by.self

##                           self.time self.pct total.time total.pct
## "lm.fit"                      11.32    77.75      11.32     77.75
## "chol2inv"                     2.05    14.08       2.05     14.08
## "runif"                        0.25     1.72       0.25      1.72
## ".External2"                   0.21     1.44       0.40      2.75
## "[.data.frame"                 0.13     0.89       0.14      0.96
## "matrix"                       0.10     0.69       0.35      2.40
## ".External"                    0.10     0.69       0.10      0.69
## "as.vector"                    0.09     0.62       0.09      0.62
## ".subset2"                     0.05     0.34       0.05      0.34
## "makepredictcall.default"      0.04     0.27       0.04      0.27
## ".deparseOpts"                 0.03     0.21       0.05      0.34
## "model.frame.default"          0.02     0.14       0.45      3.09
## "na.omit.data.frame"           0.02     0.14       0.19      1.30
## "as.data.frame.matrix"         0.02     0.14       0.11      0.76
## "is.na"                        0.02     0.14       0.02      0.14
## "pmatch"                       0.02     0.14       0.02      0.14
## "sys.call"                     0.02     0.14       0.02      0.14
## "eval"                         0.01     0.07      14.56    100.00
## "summary.lm"                   0.01     0.07       2.06     14.15
## "[["                           0.01     0.07       0.10      0.69
## "[[.data.frame"                0.01     0.07       0.09      0.62
## "<Anonymous>"                  0.01     0.07       0.06      0.41
## "makepredictcall"              0.01     0.07       0.05      0.34
## "is.factor"                    0.01     0.07       0.01      0.07

summaryRprof zwraca listę 4-elementową, gdzie by.self to informacja, która funkcja jak długo działała, uporządkowana zgodnie z czasem działania wyłącznie danej funkcji bez podfunkcji.

Graficznie: potrzebny dodatkowy pakiet profr

library(profr)
Rprof("profiler.out", interval=0.01)
tmp=table(outer(1:2000,1:2000,"*") %% 10)
tmp=sort(rnorm(10^7))
Rprof(NULL)
wynik=parse_rprof("profiler.out")
plot(wynik)

3.1 Zadanie: Profiler - paproć Barnsleya

Wykorzystać powyższy przykład użycia profilera do zbadania kody paproci Barnsleya (lub zbioru Mandelbrota).

Moje wyniki dla porównania

##                       self.time self.pct total.time total.pct
## "ifs"                      0.02    16.67       0.09     75.00
## "sample.int"               0.02    16.67       0.06     50.00
## "all"                      0.02    16.67       0.02     16.67
## "eval"                     0.01     8.33       0.12    100.00
## "withCallingHandlers"      0.01     8.33       0.12    100.00
## "stopifnot"                0.01     8.33       0.03     25.00
## "c"                        0.01     8.33       0.01      8.33
## "length"                   0.01     8.33       0.01      8.33
## "match.arg"                0.01     8.33       0.01      8.33

4 Formatowanie tabel w pakiecie `gt`

Zaczniemy od załadowania pakietu gt oraz dplyr. Możliwe, że pakiet gt trzeba zainstalować jednorazowo w znany sposób.

#install.packages("gt")
library(dplyr)
library(gt)

Podstawy dotyczące formatowania tabel w pakiecie gt znajdziemy na stronie:

podstawy gt.

W powyższym tutorialu autorzy korzystają z operatora |> mającego takie same działanie jak znany z dplyr operator %>%.

4.1 Zadanie - podstawy `gt()`

Do zadania wykorzystamy pierwsze 10 wierszy wbudowanej w R ramki danych iris, tj.

gt_iris=head(iris,10)
gt(gt_iris) # podstawowa ramka

Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa
4.6	3.4	1.4	0.3	setosa
5.0	3.4	1.5	0.2	setosa
4.4	2.9	1.4	0.2	setosa
4.9	3.1	1.5	0.1	setosa

Następnie, wykorzystując tutorial, należy (nie koniecznie po kolei):

zmienić tytuł i podtytuł tabeli na własny,
dodać źródło,
dodać odnośnik do wybranego elementu tabeli,
pogrupować kolumny działki kielicha i płatka,
przenieść kolumnę z gatunkiem na początek,
nadać kolumną własne, polskie nazwy.

Efekt końcowy powinien przypominać tabelkę wyglądającą jak niżej:

Przykładowy efekt

Gatunek	Działka kielicha		Płatek
Ramka iris
dane wbudowane w R
Gatunek	Długość	Szerokość	Długość	Szerokość
setosa	5.1	3.5	1.4	0.2
setosa	4.9	3.0	1.4	0.2
setosa	4.7	¹ 3.2	1.3	0.2
setosa	4.6	3.1	1.5	0.2
setosa	5.0	3.6	1.4	0.2
setosa	5.4	3.9	1.7	0.4
setosa	4.6	3.4	1.4	0.3
setosa	5.0	3.4	1.5	0.2
setosa	4.4	2.9	1.4	0.2
setosa	4.9	3.1	1.5	0.1
Źródło: R
¹ Mój ulubiony wiersz

4.2 Zadanie - format liczb i kolory

Zaczniemy od pobrania pliku RK4.csv ze strony przedmiotu oraz wczytaniu danych. Plik ten zawiera wyniki eksperymentów numerycznych - metody Rungego - Kutty rzędu 4 zastosowanej do wybranego równania różniczkowego zwyczajnego. Jest to pewne rozwinięcie schematu Eulera omawianego w Laboratorium 1. Kolejne kolumny oznaczają: no. - kolejny numer, N - liczba kroków, h - długość kroku, error - błąd mierzony jako różnica między wynikiem z metody a rozwiązaniem dokładnym, q - iloraz błędów, p \(=\log_2 (q)\) numeryczny rząd zbieżności metody (dobrze żeby był w przybliżeniu równy 4).

Zadanie:

wczytać dane
nadać odpowiedni format liczbowy kolumnom: no. i N - liczby całkowite - czyli pewnie zostawić, jak jest, h i error - format naukowy z 2 cyframi po przecinku, q i p - liczbowy z 2 cyframi po przecinku.
NA - zastąpić znakiem -- (dwa minusy)
wyrównanie do prawej kolumn (o ile byłoby inne)
nadać wybrane kolory tabeli (inwencja własna)
można dodać tytuł i podtytuł tabeli, jak w przykładzie niżej

Przydatne strony (sugeruje przejście od razu do przykładów)

zastępowanie NA
wyrównanie kolumn
format liczbowy - formatowanie liczb. Format naukowy (scientific) oraz inne w sekcji See also na dole strony
wybór kolorów - gotowe “style”

Przykładowy efekt

no.	N	h	error	q	p
Schemat Rungego-Kutty
bez formatowania
1	10	2.000000e-01	6.063362e-05	NA	NA
2	20	1.000000e-01	7.653130e-06	7.9227218	2.9859961
3	40	5.000000e-02	6.691689e-07	11.4367682	3.5156075
4	80	2.500000e-02	5.115456e-08	13.0813146	3.7094356
5	160	1.250000e-02	3.564025e-09	14.3530311	3.8432835
6	320	6.250000e-03	2.355794e-10	15.1287639	3.9192222
7	640	3.125000e-03	1.514799e-11	15.5518539	3.9590147
8	1280	1.562500e-03	9.616752e-13	15.7516740	3.9774333
9	2560	7.812500e-04	6.183942e-14	15.5511670	3.9589509
10	5120	3.906250e-04	8.049117e-15	7.6827586	2.9416244
11	10240	1.953125e-04	4.385381e-15	1.8354430	0.8761283
12	20480	9.765625e-05	1.598721e-14	0.2743056	-1.8661443
13	40960	4.882813e-05	9.547918e-15	1.6744186	0.7436602
14	81920	2.441406e-05	9.825474e-15	0.9717514	-0.0413408
15	163840	1.220703e-05	4.796163e-14	0.2048611	-2.2872820

no.	N	h	error	q	p
Schemat Rungego-Kutty
efekt końcowy
1	10	2.00 × 10⁻¹	6.06 × 10⁻⁵	–	–
2	20	1.00 × 10⁻¹	7.65 × 10⁻⁶	7.92	2.99
3	40	5.00 × 10⁻²	6.69 × 10⁻⁷	11.44	3.52
4	80	2.50 × 10⁻²	5.12 × 10⁻⁸	13.08	3.71
5	160	1.25 × 10⁻²	3.56 × 10⁻⁹	14.35	3.84
6	320	6.25 × 10⁻³	2.36 × 10⁻¹⁰	15.13	3.92
7	640	3.13 × 10⁻³	1.51 × 10⁻¹¹	15.55	3.96
8	1280	1.56 × 10⁻³	9.62 × 10⁻¹³	15.75	3.98
9	2560	7.81 × 10⁻⁴	6.18 × 10⁻¹⁴	15.55	3.96
10	5120	3.91 × 10⁻⁴	8.05 × 10⁻¹⁵	7.68	2.94
11	10240	1.95 × 10⁻⁴	4.39 × 10⁻¹⁵	1.84	0.88
12	20480	9.77 × 10⁻⁵	1.60 × 10⁻¹⁴	0.27	−1.87
13	40960	4.88 × 10⁻⁵	9.55 × 10⁻¹⁵	1.67	0.74
14	81920	2.44 × 10⁻⁵	9.83 × 10⁻¹⁵	0.97	−0.04
15	163840	1.22 × 10⁻⁵	4.80 × 10⁻¹⁴	0.20	−2.29

4.3 Więcej `gt`

Więcej przykładów zastosowania pakietu gt i formatowania tabel można zaleźć w np.:

https://posit.co/blog/rstudio-community-table-gallery/

https://forum.posit.co/c/irl/table-gallery/64

Przetwarzanie i wizualizacja danych w R

Laboratorium 4 - microbenchmark, profiler, formatowanie tabel z pakietem gt

Michał Braś

11 listopad, 2025

1 Testowanie czasu wykonywania kodu

1.1 Zadanie - czas działania kodu

2 Microbenchmark

2.1 Zadania - microbenchmark

2.1.1 Zadanie microbenchmark - norma wektora

2.1.2 Zadanie microbenchmark - norma wektora - czas działania w zależności od liczby elementów

2.1.3 Zadanie microbenchmark - `for` vs. `sapply` vs. `vapply` vs. `wektorowo`

2.1.4 Zadanie microbenchmark - `runif`

2.1.5 Zadanie microbenchmark - paproć Barnsleya lub zbiór Mandelbrota

3 Testowanie składowych bloków - profiler

3.1 Zadanie: Profiler - paproć Barnsleya

4 Formatowanie tabel w pakiecie `gt`

4.1 Zadanie - podstawy `gt()`

4.2 Zadanie - format liczb i kolory

4.3 Więcej `gt`

Przetwarzanie i wizualizacja danych w R

Laboratorium 4 - microbenchmark, profiler, formatowanie tabel z pakietem gt

Michał Braś

11 listopad, 2025

1 Testowanie czasu wykonywania kodu

1.1 Zadanie - czas działania kodu

2 Microbenchmark

2.1 Zadania - microbenchmark

2.1.1 Zadanie microbenchmark - norma wektora

2.1.2 Zadanie microbenchmark - norma wektora - czas działania w zależności od liczby elementów

2.1.3 Zadanie microbenchmark - for vs. sapply vs. vapply vs. wektorowo

2.1.4 Zadanie microbenchmark - runif

2.1.5 Zadanie microbenchmark - paproć Barnsleya lub zbiór Mandelbrota

3 Testowanie składowych bloków - profiler

3.1 Zadanie: Profiler - paproć Barnsleya

4 Formatowanie tabel w pakiecie gt

4.1 Zadanie - podstawy gt()

4.2 Zadanie - format liczb i kolory

4.3 Więcej gt

2.1.3 Zadanie microbenchmark - `for` vs. `sapply` vs. `vapply` vs. `wektorowo`

2.1.4 Zadanie microbenchmark - `runif`

4 Formatowanie tabel w pakiecie `gt`

4.1 Zadanie - podstawy `gt()`

4.3 Więcej `gt`