Функции system.time(), benchmark(), microbenchmark() позволяют оценить общее время выполнения выражения и нивелировать возможные вариации за счет множества попыток, однако для более полной оценки и определения стратегии оптимизации кода необходимо также выявлять и «узкие» места в выполняемом коде. К подобным «узким» местам можно отнести те вызовы, которые занимают наибольшее количество времени или расходуют больше всего памяти (при профилирования расхода памяти). Профилирование кода позволяет собрать статистику работы кода в целом, а также всех выполняемых вызовов.

Существуют специальные пакеты, предназначенные для организации профилирования или для более наглядной демонстрации результатов, однако все они основываются на использовании функции Rprof() из пакета base.

Пакет base

Принцип работы функции Rprof() таков: с заданным интервалом (аргумент interval) делаются «снимки» вызовов и записываются в лог-файл (аргумент filename). Следует помнить о том, что Rprof() делает «снимки» каждые 0.02 секунды (значение по умолчанию), поэтому если выражение выполняется быстрее, то Rprof() не сможет его отследить. В таком случае периодичность «снимков» можно задать вручную с помощью аргумента interval. Также можно увеличить объём данных, чтобы код выполнялся дольше.

Ещё одной специфической особенностью профилирования является то, что при профилирования записываются абсолютно все вызовы. Например, если вы используется функцию replicate() для организации повторного выполнения выражения, то вызовы replicate() также будут записаны в лог, что может затруднить его анализ. Аналогичная ситуация будет происходить при вызове скрипта через функцию source()

Обратим внимание, что при профилировании мы можем использовать многократное выполнение анализируемой функции для сбора более надёжной статистики вызовов.

Схематически профилирование кода со значениями аргументов по умолчанию можно представить следующим образом:

Rprof() # включаем профилирование
# наш код, скрипт и т.д.
Rprof(NULL) # отключаем проифлирование
summaryRprof() # выводим статистику

В качестве примера, ниже представлены шаги процесса профилирования функции, которая рассчитывает первичные описательные статистики (среднее, медиану, стандартное отклонение, минимум и максимум) по сгенерированным \(10^6\) нормально распределённым значениям.

Итак, для начала сгенерируем вектор нормально распределённых данных и добавим немного пропущенных значений:

x <- rnorm(10^6)
x[sample.int(length(x), size = 10^3)] <- NA

Объявим данную функцию, которая будет расчитывать описательные статистики:

desc <- function(x) {
    x <- na.omit(x)
    n <- length(x)
    mean <- mean(x)
    median <- median(x)
    sd <- sd(x)
    min <- min(x)
    max <- max(x)
    return(c(n = n,  mean = mean, median = median, sd = sd, minx = min, max = max))
}

Теперь мы можем приступить непосредственно к профилированию.

Rprof() # Включаем профилирование
# Выполняем код несколько раз для получения более надёжных результатов
for (i in seq_len(100)) invisible(desc(x))
Rprof(NULL) # отключаем профилирование

Рассмотрим результаты профилирования функции desc():

summaryRprof()
#> $by.self
#>                   self.time self.pct total.time total.pct
#> "sort.int"             4.52    41.77       4.98     46.03
#> "na.omit.default"      3.60    33.27       4.04     37.34
#> "is.na"                0.92     8.50       0.92      8.50
#> ".Call"                0.58     5.36       0.58      5.36
#> "any"                  0.26     2.40       0.26      2.40
#> "mean.default"         0.24     2.22       0.24      2.22
#> "min"                  0.24     2.22       0.24      2.22
#> "max"                  0.20     1.85       0.20      1.85
#> "as.double"            0.16     1.48       0.16      1.48
#> "seq_along"            0.10     0.92       0.10      0.92
#> 
#> $by.total
#>                   total.time total.pct self.time self.pct
#> "desc"                 10.82    100.00      0.00     0.00
#> "median"                5.36     49.54      0.00     0.00
#> "median.default"        5.36     49.54      0.00     0.00
#> "mean"                  5.22     48.24      0.00     0.00
#> "sort.int"              4.98     46.03      4.52    41.77
#> "sort"                  4.98     46.03      0.00     0.00
#> "sort.default"          4.98     46.03      0.00     0.00
#> "na.omit.default"       4.04     37.34      3.60    33.27
#> "na.omit"               4.04     37.34      0.00     0.00
#> "is.na"                 0.92      8.50      0.92     8.50
#> "sd"                    0.74      6.84      0.00     0.00
#> "var"                   0.74      6.84      0.00     0.00
#> ".Call"                 0.58      5.36      0.58     5.36
#> "any"                   0.26      2.40      0.26     2.40
#> "mean.default"          0.24      2.22      0.24     2.22
#> "min"                   0.24      2.22      0.24     2.22
#> "max"                   0.20      1.85      0.20     1.85
#> "as.double"             0.16      1.48      0.16     1.48
#> "is.data.frame"         0.16      1.48      0.00     0.00
#> "seq_along"             0.10      0.92      0.10     0.92
#> 
#> $sample.interval
#> [1] 0.02
#> 
#> $sampling.time
#> [1] 10.82

Результатом подобных операций является 2 таблицы, отсортированные по разным основаниям: по времени выполнения функции (self.time) и общему времени выполнения функции (total.time). В столбцах представлены абсолютные и относительные (в процентах) показатели времени выполнения.

При профилировании скриптов включение опции line.profiling проставляет номера строк, в которых производился вызов, что позволяет проанализировать производительность всего кода:

Rprof(line.profiling = TRUE)
source("/path/scritpt.R")
Rprof(NULL)
summaryRprof()

Пакет proftools

Данный пакет представляет ряд функций для более наглядного представления результатов профилирования. Для использования графических возможностей пакета, нам необходимо установить пакет graph, Rgraphviz из Bioconductor. Сделать это можно следующим образом:

source("http://bioconductor.org/biocLite.R")
biocLite("graph")
biocLite("Rgraphviz")

Функция flatProfile() выводит таблицу, сходную с таблицей «by.total» возвращаемую функцией summaryRprof(), в одну и сортирует их по общему времени выполнения в процентах или по времени выполнения функции в процентах, если указать byTotal = FALSE. Перед использование функции flatProfile() нам необходимо загрузить и обработать лог, созданный в результате профилирования, с помощью функции readProfileData():

library(proftools)
prof_data <- readProfileData("Rprof.out")
flatProfile(prof_data)
#>                 total.pct total.time self.pct self.time
#> desc               100.00      10.82     0.00      0.00
#> median              49.54       5.36     0.00      0.00
#> median.default      49.54       5.36     0.00      0.00
#> mean                48.24       5.22     0.00      0.00
#> sort                46.03       4.98     0.00      0.00
#> sort.default        46.03       4.98     0.00      0.00
#> sort.int            46.03       4.98    41.77      4.52
#> na.omit             37.34       4.04     0.00      0.00
#> na.omit.default     37.34       4.04    33.27      3.60
#> is.na                8.50       0.92     8.50      0.92
#> sd                   6.84       0.74     0.00      0.00
#> var                  6.84       0.74     0.00      0.00
#> .Call                5.36       0.58     5.36      0.58
#> any                  2.40       0.26     2.40      0.26
#> mean.default         2.22       0.24     2.22      0.24
#> min                  2.22       0.24     2.22      0.24
#> max                  1.85       0.20     1.85      0.20
#> as.double            1.48       0.16     1.48      0.16
#> is.data.frame        1.48       0.16     0.00      0.00
#> seq_along            0.92       0.10     0.92      0.10

Функция интересна тем, что предоставляет более компактный вывод при той же информативности, что и функция summaryRprof().

Пакет proftools содержит функцию printProfileCallGraph(), которая позволяет наглядн представить всю иерархию дочерних вызовов.

printProfileCallGraph(prof_data)
#> Call graph
#> 
#> index    % time     % self   % children     name
#> 
#> [1]      100.00       0.00     100.00     desc [1]    
#>                       1.85       0.00         max [17] 
#>                       0.00       2.22         mean [4] 
#>                       0.00      49.54         median [2] 
#>                       2.22       0.00         min [16] 
#>                       0.00      37.34         na.omit [8] 
#>                       0.00       6.84         sd [11] 
#> -----------------------------------------------
#>                       0.00      49.54         desc [1] 
#> [2]       49.54       0.00      49.54     median [2]    
#>                       0.00      49.54         median.default [3] 
#> -----------------------------------------------
#>                       0.00      49.54         median [2] 
#> [3]       49.54       0.00      49.54     median.default [3]    
#>                       0.92       0.00         any [14] 
#>                       2.59       0.00         is.na [10] 
#>                       0.00      46.03         mean [4] 
#> -----------------------------------------------
#>                       0.00      46.03         median.default [3] 
#>                       0.00       2.22         desc [1] 
#> [4]       48.24       0.00      48.24     mean [4]    
#>                       2.22       0.00         mean.default [15] 
#>                       0.00      46.03         sort [5] 
#> -----------------------------------------------
#>                       0.00      46.03         mean [4] 
#> [5]       46.03       0.00      46.03     sort [5]    
#>                       0.00      46.03         sort.default [6] 
#> -----------------------------------------------
#>                       0.00      46.03         sort [5] 
#> [6]       46.03       0.00      46.03     sort.default [6]    
#>                      41.77       4.25         sort.int [7] 
#> -----------------------------------------------
#>                      41.77       4.25         sort.default [6] 
#> [7]       46.03      41.77       4.25     sort.int [7]    
#>                       1.48       0.00         any [14] 
#>                       2.77       0.00         is.na [10] 
#> -----------------------------------------------
#>                       0.00      37.34         desc [1] 
#> [8]       37.34       0.00      37.34     na.omit [8]    
#>                      33.27       4.07         na.omit.default [9] 
#> -----------------------------------------------
#>                      33.27       4.07         na.omit [8] 
#> [9]       37.34      33.27       4.07     na.omit.default [9]    
#>                       3.14       0.00         is.na [10] 
#>                       0.92       0.00         seq_along [20] 
#> -----------------------------------------------
#>                       2.77       0.00         sort.int [7] 
#>                       3.14       0.00         na.omit.default [9] 
#>                       2.59       0.00         median.default [3] 
#> [10]       8.50       8.50       0.00     is.na [10]   
#> -----------------------------------------------
#>                       0.00       6.84         desc [1] 
#> [11]       6.84       0.00       6.84     sd [11]   
#>                       0.00       6.84         var [12] 
#> -----------------------------------------------
#>                       0.00       6.84         sd [11] 
#> [12]       6.84       0.00       6.84     var [12]   
#>                       5.36       0.00         .Call [13] 
#>                       0.00       1.48         is.data.frame [19] 
#> -----------------------------------------------
#>                       5.36       0.00         var [12] 
#> [13]       5.36       5.36       0.00     .Call [13]   
#> -----------------------------------------------
#>                       1.48       0.00         sort.int [7] 
#>                       0.92       0.00         median.default [3] 
#> [14]       2.40       2.40       0.00     any [14]   
#> -----------------------------------------------
#>                       2.22       0.00         mean [4] 
#> [15]       2.22       2.22       0.00     mean.default [15]   
#> -----------------------------------------------
#>                       2.22       0.00         desc [1] 
#> [16]       2.22       2.22       0.00     min [16]   
#> -----------------------------------------------
#>                       1.85       0.00         desc [1] 
#> [17]       1.85       1.85       0.00     max [17]   
#> -----------------------------------------------
#>                       1.48       0.00         is.data.frame [19] 
#> [18]       1.48       1.48       0.00     as.double [18]   
#> -----------------------------------------------
#>                       0.00       1.48         var [12] 
#> [19]       1.48       0.00       1.48     is.data.frame [19]   
#>                       1.48       0.00         as.double [18] 
#> -----------------------------------------------
#>                       0.92       0.00         na.omit.default [9] 
#> [20]       0.92       0.92       0.00     seq_along [20]   
#> -----------------------------------------------

Каждый вызов здесь пронумерован (в квадратных скобках), что позволяет детально проанализировать каждый вызов и его дочерние вызовы, а также время их выполнения.

Также пакет proftools позволяет представить иерархию вызовов в виде графика, что существенно упрощает понимание структуры вызовов:

plotProfileCallGraph(prof_data)

Стоит отметить, что по графику мы не можем оценить время, которое занял тот или иной вызов.

Функция proftable()

Рассмотрим ещё один способ представления результатов профилирования — функция proftable(), написанная Noam Ross. Исходный код данной функции доступен под лицензией GNU GPL v2 и размещён в открытом доступе на github.

Я переработал данную функцию, удалив зависимость от пакета plyr и существенно ускорив её работу, при этом сохранив возвращаемый результат. Код модифицированного варианта доступен в git-репозитории

. Импортировать данный скрипт можно с помощью функции devtools::source_gist():

devtools::source_gist("ce84b9f1e1c0f3af6d64")

Данная функция работает с лог-файлом, полученным в результате профилирования и используется как дополнение к функции summaryRprof().

Для более наглядного представления результатов профилирование необходимо провести с включённой опцией line.profiling. Для удобства последующего анализа результатов можно сохранить функцию и её вызов в файл скрипта и вызывать его с помощью функции source(). В нашем примере скрипт для профилирования имеет следующее содержимое:

desc <- function(x) {
    x <- na.omit(x)
    n <- length(x)
    mean <- mean(x)
    median <- median(x)
    sd <- sd(x)
    min <- min(x)
    max <- max(x)
    return(c(n = n,  mean = mean, median = median, sd = sd, minx = min, max = max))
}

for (i in seq_len(100)) desc(x)

Данный файл можно сформировать прямо из консоли R:

script_file <- tempfile(fileext = ".R")
dump("desc", file = script_file)
write("x <- rnorm(10^6)", file = script_file, append = TRUE)
write("for (i in seq_len(100)) desc(x)", file = script_file, append = TRUE)

Процедура профилирования осуществляется описанным выше способом:

Rprof(line.profiling = TRUE)
source(script_file)
Rprof(NULL)

Результаты работы функции proftable() представлены ниже.

proftable("Rprof.out")
#> Calls:
#>  RealTime PctTime Call                                                                          
#>  3.56     44.61   desc > median > median.default > mean > sort > sort.default > sort.int        
#>  1.56     19.55   desc > na.omit > na.omit.default                                              
#>  0.56      7.02   desc > sd > var > .Call                                                       
#>  0.34      4.26   desc > na.omit > na.omit.default > is.na                                      
#>  0.34      4.26   desc > median > median.default > mean > sort > sort.default > sort.int > is.na
#>  0.30      3.76   desc > min                                                                    
#>  0.24      3.01   desc > na.omit > na.omit.default > seq_along                                  
#>  0.22      2.76   desc > median > median.default > any                                          
#>  0.22      2.76   desc > median > median.default > is.na                                        
#>  0.22      2.76   desc > max                                                                    
#> 
#> Files: None 
#> 
#> Parent Call: source > withVisible > eval > eval 
#> 
#> Total Time: 7.98 seconds
#> 
#> Percent of run time represented: 94.7%

Разберём вывод функции proftable(). В первом столбце отображается время выполнения того или иного вызова, во втором — название вызываемой функции. Далее по цепочке приведена иерархия вызовов. Весь вывод отсортирован по времени выполнения вызовов по убыванию. Вывод наглядно демонстрирует наиболее медленные вызовы вплоть до самого низкого уровня иерархии.

Обратите внимание, что функция proftable() имеет второй аргумент — lines, который отвечает за количество строк в выводе (значение по умолчанию — 10).

Пакет profr

Пакет profr предоставляет несколько функций для упрощения процесса профилирования, а также позволяет графически представить результаты профилирования. Функция profr() использует рассмотренную ранее функцию Rprof() для организации профилирования. Результаты работы функции profr() представляют собой таблицу, с хронологическим перечислением вызовов. Поэтому при использовании данной функции необходимо использовать однократное выполнение выражения. Вышеприведённый листинг по профилированию работы функции desc() будет выглядеть следующим образом:

library(profr)
prof_res <- profr(desc(x))
print(prof_res)
#>    level g_id t_id               f start  end n  leaf time     source
#> 24     1    1    1            desc  0.00 0.10 1 FALSE 0.10 .GlobalEnv
#> 25     2    1    1         na.omit  0.00 0.02 1 FALSE 0.02  S4Vectors
#> 26     2    2    1          median  0.02 0.08 1 FALSE 0.06    IRanges
#> 27     2    3    1             min  0.08 0.10 1  TRUE 0.02       base
#> 28     3    1    1 standardGeneric  0.00 0.02 1 FALSE 0.02       base
#> 29     3    2    1 standardGeneric  0.02 0.08 1 FALSE 0.06       base
#> 30     4    1    1         na.omit  0.00 0.02 1 FALSE 0.02  S4Vectors
#> 31     4    2    1          median  0.02 0.08 1 FALSE 0.06    IRanges
#> 32     5    1    1 na.omit.default  0.00 0.02 1  TRUE 0.02       <NA>
#> 33     5    2    1  median.default  0.02 0.08 1 FALSE 0.06      stats
#> 34     6    1    1            mean  0.02 0.08 1 FALSE 0.06    IRanges
#> 35     7    1    1            sort  0.02 0.08 1 FALSE 0.06    IRanges
#> 36     8    1    1    sort.default  0.02 0.08 1 FALSE 0.06       base
#> 37     9    1    1        sort.int  0.02 0.08 1 FALSE 0.06       base
#> 38    10    1    1           is.na  0.02 0.04 1  TRUE 0.02       base

Содержание столбцов: * «f» — название функции; * «level» — уровень в иерархии вызовов; * «time» — общее время выполнения функции; * «start» — время начала выполнения функции; * «end» — время окончания выполнения функции; * «leaf» — TRUE, если функция вызывается другими функциями; * «source» — название пакеты, содержащего данную функцию.

Графическое представление результатов профилирования осуществляется с помощью стандартной функции plot() или ggplot():

plot(prof_res)

«Стаки» (вложенные вызовы) одного вызова на графике расположены вертикально. Ширина столбцов пропорциональна относительному времени выполнения данного вызова.

Пакет timeit

Пакет timeit также является «обёрткой» для функций Rprof() и summaryRprof(), но при этом имеет встроенную поддержку многократного повторения кода для случаев, когда код выполняется слишком быстро. Для профилирования используется функция timeit(). Есть возможность задать временной интервалам между «снимками» (аргумент interval) и количество повторений (аргумент replications).

library(timeit)
prof_res <- timeit(desc(x), replications = 100)
#> Running iteration 1 of 10 ...
#> Running iteration 2 of 10 ...
#> Running iteration 3 of 10 ...
#> Running iteration 4 of 10 ...
#> Running iteration 5 of 10 ...
#> Running iteration 6 of 10 ...
#> Running iteration 7 of 10 ...
#> Running iteration 8 of 10 ...
#> Running iteration 9 of 10 ...
#> Running iteration 10 of 10 ...
head(prof_res)
#>                   self.time self.pct total.time total.pct mem.total replications iteration              func
#> "na.omit.default"    0.0519   38.388     0.0548    16.723    13.557          100         1 "na.omit.default"
#> "sort.int"           0.0499   36.908     0.0542    16.540     9.657          100         1        "sort.int"
#> "is.na"              0.0103    7.618     0.0103     3.143     4.014          100         1           "is.na"
#> "as.double"          0.0072    5.325     0.0072     2.197     0.307          100         1       "as.double"
#> ".Call"              0.0052    3.846     0.0052     1.587     2.065          100         1           ".Call"
#> "any"                0.0036    2.663     0.0036     1.099     1.950          100         1             "any"

Для получения усреднённой сводной таблицы результатов профилирования используется функция mean():

mean(prof_res)
#>                 func self.time  self.pct total.time total.pct mem.total replications iteration iterations
#> 22        "sort.int"   0.05011 37.644866    0.05401 18.272330    9.2948         1000       5.5         10
#> 21 "na.omit.default"   0.04932 37.051271    0.05228 17.742028   13.4959         1000       5.5         10
#> 20           "is.na"   0.01110  8.339944    0.01110  3.742493    4.5291         1000       5.5         10
#> 19       "as.double"   0.00736  5.530378    0.00736  2.498459    0.4138         1000       5.5         10
#> 18           ".Call"   0.00465  3.492723    0.00465  1.574199    1.6902         1000       5.5         10
#> 17             "min"   0.00287  2.156964    0.00287  0.975188    1.0939         1000       5.5         10
#> 16             "max"   0.00266  1.998433    0.00266  0.903927    0.3991         1000       5.5         10
#> 15             "any"   0.00262  1.967254    0.00262  0.902661    1.2125         1000       5.5         10
#> 14    "mean.default"   0.00162  1.218017    0.00162  0.551193    1.0656         1000       5.5         10
#> 13       "seq_along"   0.00049  0.367852    0.00049  0.171103    0.1111         1000       5.5         10
#> 6             "mean"   0.00008  0.059655    0.03379  9.966799    6.1880          600       3.1          6
#> 7           "median"   0.00004  0.029957    0.01838  6.302348    3.5124          300       1.9          3
#> 9             "sort"   0.00004  0.030026    0.02184  6.347843    3.5652          400       1.5          4
#> 10 "standardGeneric"   0.00004  0.030034    0.05113 14.489284   11.5146          400       1.4          4
#> 8          "na.omit"   0.00003  0.022540    0.01546  5.033082    4.0519          300       2.2          3
#> 1             "desc"   0.00002  0.015061    0.02656  5.866367    5.9240          200       1.2          2
#> 2             "eval"   0.00001  0.007508    0.01332  4.053561    3.0708          100       0.7          1
#> 3    "is.data.frame"   0.00001  0.007435    0.00086  0.267413    0.0383          100       0.4          1
#> 4        "is.factor"   0.00001  0.007513    0.00001  0.004978    0.0076          100       0.5          1
#> 5        "is.finite"   0.00001  0.007508    0.00001  0.003043    0.0077          100       0.7          1
#> 11       "stopifnot"   0.00001  0.007564    0.00001  0.004013    0.0076          100       0.6          1
#> 12             "var"   0.00001  0.007496    0.00149  0.327689    0.1532          100       1.0          1

Вывод функции mean.timeit() сходен с выводом функции summaryRprof() и даёт представление об общем и относительном врмени выполнении каждого вызова. Также пакет timeit содержит функцию для графического представления результатов, но график получается не очень информативным. Судите сами:

plot(prof_res)

Выводы

Из примера выше мы можем сделать вывод, что наиболее медленной функцией из состава функции desc() является функция median(), которая через некоторую последовательность вызовов вызывает функцию sort.int(). Именно функция sort.int() является наиболее «узким» местом всех расчётов.

Иногда, чтобы выяснить какой функции принадлежит тот или иной вызов, приходится анализировать содержимое функций. Чтобы увидеть содержимое функции, можно воспользоваться функцией getAnywhere(). Например:

getAnywhere(median.default)
#> A single object matching 'median.default' was found
#> It was found in the following places
#>   package:stats
#>   registered S3 method for median from namespace stats
#>   namespace:stats
#> with value
#> 
#> function (x, na.rm = FALSE) 
#> {
#>     if (is.factor(x) || is.data.frame(x)) 
#>         stop("need numeric data")
#>     if (length(names(x))) 
#>         names(x) <- NULL
#>     if (na.rm) 
#>         x <- x[!is.na(x)]
#>     else if (any(is.na(x))) 
#>         return(x[FALSE][NA])
#>     n <- length(x)
#>     if (n == 0L) 
#>         return(x[FALSE][NA])
#>     half <- (n + 1L)%/%2L
#>     if (n%%2L == 1L) 
#>         sort(x, partial = half)[half]
#>     else mean(sort(x, partial = half + 0L:1L)[half + 0L:1L])
#> }
#> <bytecode: 0x4bc02a8>
#> <environment: namespace:stats>

Если искомая функция принадлежит какому-то пакету, то этот пакет должен быть предварительно загружен.

С учётом выявленных «узких» мест нашей функции и проанализировав исходный код функций, мы можем оптимизировать нашу исходную функцию для достижения максимальной производительности. Наш вариант оптимизированной функции будет выглядеть следюущим образом:

desc <- function(x) {
    x <- x[!is.na(x)]
    n <- length(x)
    mean <- sum(x) / n
    sd <- sqrt(sum((x - mean)^2) / (n - 1))
    min <- min(x)
    max <- max(x)
    half <- (n + 1L) %/% 2L
    if (n %% 2L == 1L)
        median <- .Internal(psort(x, half))[half]
    else
        median <- sum(.Internal(psort(x, c(half, half + 1)))[c(half, half + 1)]) / 2L
    return(c(n, mean, median, sd, min, max))
}

Вместо функций mean() и sd() мы использовали общеизвестные формулы. А для расчёта медианы была взята часть оригинальной функции median.default, в которой вызов функций sort() был заменён на более быстрый вызов .internal(psort). Также мы заменили медленную функцию na.omit() на более быстрый вариант.