Документация Engee

Советы по производительности

В следующих разделах мы вкратце рассмотрим ряд приемов, которые могут помочь максимально ускорить выполнение кода на Julia.

Код, для которого важна производительность, следует помещать в функции

Любой код, для которого важна производительность, следует помещать в функции. Код внутри функций, как правило, выполняется намного быстрее, чем код верхнего уровня, из-за особенностей работы компилятора Julia.

Функции важны не только с точки зрения производительности: они лучше поддаются повторному использованию и тестированию, а также упрощают понимание того, какие действия выполняются и каковы входные и выходные данные. В руководстве по стилю Julia также рекомендуется писать функции, а не просто скрипты.

Функции должны принимать аргументы, а не работать напрямую с глобальными переменными (см. следующий пункт).

Избегайте нетипизированных глобальных переменных

Значение нетипизированной глобальной переменной может измениться в любой момент, что может привести к изменению ее типа. Поэтому компилятору трудно оптимизировать код, в котором применяются глобальные переменные. Это также относится к переменным, значениями которых являются типы, то есть к псевдонимам типов на глобальном уровне. Переменные должны быть локальными или передаваться в функции как аргументы, где это возможно.

Глобальные имена часто являются константами, и если объявлять их таким образом, можно существенно повысить производительность:

const DEFAULT_VAL = 0

Если известно, что тип глобальной переменной меняться не будет, его следует аннотировать.

Использование нетипизированных глобальных переменных можно оптимизировать, аннотировав их типы в момент использования:

global x = rand(1000)

function loop_over_global()
    s = 0.0
    for i in x::Vector{Float64}
        s += i
    end
    return s
end

С точки зрения стиля лучше передавать аргументы в функции. Это улучшает возможность повторного использования кода и позволяет лучше понять входные и выходные данные.

Весь код в REPL вычисляется в глобальной области, поэтому переменная, определенная и заданная на верхнем уровне, будет глобальной. Переменные, определенные в области верхнего уровня внутри модулей, также являются глобальными.

В сеансе REPL данный код:

julia> x = 1.0

эквивалентен следующему коду:

julia> global x = 1.0

поэтому применимы все описанные выше особенности, касающиеся производительности.

Измерение производительности с помощью @time и контроль выделения памяти

Для измерения производительности полезным будет макрос @time. Здесь мы повторим приведенный выше пример с глобальной переменной, но на этот раз удалим аннотацию типа:

julia> x = rand(1000);

julia> function sum_global()
           s = 0.0
           for i in x
               s += i
           end
           return s
       end;

julia> @time sum_global()
  0.011539 seconds (9.08 k allocations: 373.386 KiB, 98.69% compilation time)
523.0007221951678

julia> @time sum_global()
  0.000091 seconds (3.49 k allocations: 70.156 KiB)
523.0007221951678

При первом вызове (@time sum_global()) функция компилируется. (Если в рамках сеанса макрос @time еще не использовался, будут также скомпилированы функции, необходимые для отсчета времени.) Воспринимать всерьез результаты этого выполнения не стоит. Обратите внимание, что для второго выполнения, помимо времени, также сообщается о выделении значительного объема памяти. В данном случае просто вычисляется сумма всех элементов в векторе 64-битных чисел с плавающей запятой, поэтому выделять память (в куче) не требуется.

Следует пояснить, что макрос @time сообщает о выделении памяти именно в куче. Обычно такое выделение требуется либо для изменяемых объектов, либо для создания или увеличения контейнеров переменного размера (например, Array, Dict, строк или нестабильных по типу объектов, тип которых известен только во время выполнения). Выделение (или отмена выделения) таких блоков памяти может потребовать ресурсоемкого системного вызова (например, malloc в C) и отслеживания сборщиком мусора. Напротив, хранение неизменяемых значений, таких как числа (кроме bignum), кортежи и неизменяемые структуры (struct), требует гораздо меньше ресурсов (например, они могут храниться в стеке или памяти регистров ЦП), поэтому при выделении памяти для них не приходится беспокоиться о производительности.

Непредвиденное выделение памяти — это почти всегда признак какой-то проблемы с кодом, обычно связанной со стабильностью типов или созданием множества небольших временных массивов. Поэтому, даже если не брать во внимание расход памяти, высока вероятность того, что код, сгенерированный для функции, далек от оптимального. Обращайте внимание на такие симптомы и следуйте приведенному ниже совету.

В данном конкретном случае память выделяется вследствие использования нестабильной по типу глобальной переменной x. Если вместо этого передать x в функцию как аргумент, память больше выделяться не будет (некоторый объем выделяемой памяти ниже связан с тем, что макрос @time выполняется в глобальной области). После первого вызова такой вариант будет гораздо быстрее:

julia> x = rand(1000);

julia> function sum_arg(x)
           s = 0.0
           for i in x
               s += i
           end
           return s
       end;

julia> @time sum_arg(x)
  0.007551 seconds (3.98 k allocations: 200.548 KiB, 99.77% compilation time)
523.0007221951678

julia> @time sum_arg(x)
  0.000006 seconds (1 allocation: 16 bytes)
523.0007221951678

Выделение 1 в выходных данных обусловлено тем, что сам макрос @time выполняется в глобальной области. Если отсчет времени вести в самой функции, то память не будет выделяться вообще:

julia> time_sum(x) = @time sum_arg(x);

julia> time_sum(x)
  0.000002 seconds
523.0007221951678

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

Для более основательного тестирования производительности можно воспользоваться пакетом BenchmarkTools.jl, который, помимо прочего, оценивает функцию несколько раз, чтобы уменьшить влияние случайных факторов.

Инструменты

Язык Julia и его экосистема пакетов предлагают инструменты, которые могут помочь в диагностике проблем и повышении производительности кода:

  • Профилирование позволяет измерять производительность выполняемого кода и выявлять строки, которые являются узким местом. В сложных проектах пакет ProfileView может помочь визуализировать результаты профилирования.

  • Пакет Traceur может помочь в выявлении распространенных проблем с производительностью кода.

  • Неожиданно большой объем выделяемой памяти согласно выходным данным макроса @time, @allocated или профилировщика (при вызове подпрограмм сборки мусора) намекает на возможные проблемы с кодом. Если иной причины для такого избыточного выделения не обнаружено, в первую очередь можно предположить проблему с типами. Вы также можете запустить Julia с параметром --track-allocation=user и изучить полученные файлы *.mem, чтобы узнать, где происходит выделение памяти. См. раздел Анализ выделения памяти.

  • Макрос @code_warntype создает представление кода, которое может помочь найти выражения, приводящие к неоднозначности типов. См. раздел, посвященный макросу @code_warntype.

Избегайте контейнеров с параметрами абстрактных типов

При работе с параметризованными типами, включая массивы, лучше по возможности избегать параметризации с абстрактными типами.

Рассмотрим следующий код:

julia> a = Real[]
Real[]

julia> push!(a, 1); push!(a, 2.0); push!(a, π)
3-element Vector{Real}:
 1
 2.0
 π = 3.1415926535897...

Так как a — это массив абстрактного типа Real, в нем должна быть возможность хранить любое значение Real. Так как объекты Real могут быть произвольного размера и структуры, переменная a должна быть представлена массивом указателей на отдельные объекты Real в памяти. Однако, если бы в a было разрешено хранить числа только одного типа, например Float64, эффективность хранения повысилась бы:

julia> a = Float64[]
Float64[]

julia> push!(a, 1); push!(a, 2.0); push!(a,  π)
3-element Vector{Float64}:
 1.0
 2.0
 3.141592653589793

При присваивании чисел переменной a теперь они будут преобразовываться в тип Float64, и a будет храниться как непрерывный блок 64-битных значений с плавающей запятой, с которыми можно эффективно производить операции.

Если не получается избежать использования контейнеров с абстрактными типами значений, иногда лучше параметризировать их с типом Any, чтобы во время выполнения типы не проверялись. Например, IdDict{Any, Any} будет работать быстрее, чем IdDict{Type, Vector}.

См. также обсуждение в разделе Параметрические типы.

Объявления типов

Во многих языках с необязательными объявлениями типов добавление объявлений — основной способ повысить быстродействие кода. В Julia это не так. В Julia компилятору, как правило, известны типы всех аргументов функций, локальных переменных и выражений. Однако есть ряд особых ситуаций, в которых объявления могут быть полезны.

Избегайте полей с абстрактными типами

При объявлении полей можно не указывать их типы:

julia> struct MyAmbiguousType
           a
       end

В результате a может быть любого типа. Это часто полезно, но есть и оборотная сторона: для объектов типа MyAmbiguousType компилятор не сможет сгенерировать высокопроизводительный код. Причина в том, что при определении способа сборки кода компилятор отталкивается от типов объектов, а не от их значений. К сожалению, из типа MyAmbiguousType выяснить можно немногое:

julia> b = MyAmbiguousType("Hello")
MyAmbiguousType("Hello")

julia> c = MyAmbiguousType(17)
MyAmbiguousType(17)

julia> typeof(b)
MyAmbiguousType

julia> typeof(c)
MyAmbiguousType

Значения b и c относятся к одному типу, однако их базовое представление в памяти сильно различается. Даже если бы в поле a хранились только числовые значения, из того, что представление типов UInt8 и Float64 в памяти различается, также следует, что они должны обрабатываться процессором с использованием разных инструкций. Так как необходимой информации в типе не содержится, решение должно приниматься во время выполнения. Это снижает производительность.

Чтобы улучшить ситуацию, можно объявить тип a. Здесь мы рассматриваем ситуацию, когда a может иметь любой из нескольких типов. В таком случае напрашивающимся решением будет использование параметров. Пример:

julia> mutable struct MyType{T<:AbstractFloat}
           a::T
       end

Такой вариант лучше, чем

julia> mutable struct MyStillAmbiguousType
           a::AbstractFloat
       end

потому что в первом случае тип a определяется на основе типа объекта-оболочки. Пример:

julia> m = MyType(3.2)
MyType{Float64}(3.2)

julia> t = MyStillAmbiguousType(3.2)
MyStillAmbiguousType(3.2)

julia> typeof(m)
MyType{Float64}

julia> typeof(t)
MyStillAmbiguousType

Тип поля a можно легко определить на основании типа m, но не типа t. Ведь действительно, в t тип поля a может измениться:

julia> typeof(t.a)
Float64

julia> t.a = 4.5f0
4.5f0

julia> typeof(t.a)
Float32

Напротив, после создания m тип m.a изменить нельзя:

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float64

То, что тип m.a выводится из m, в совокупности с тем фактом, что этот тип не может измениться внутри функции, позволяет компилятору генерировать оптимизированный код для таких объектов, как m, но не t.

Безусловно, все это справедливо только в том случае, если m создается с конкретным типом. Если создать его явным образом с указанием абстрактного типа, ситуация изменится:

julia> m = MyType{AbstractFloat}(3.2)
MyType{AbstractFloat}(3.2)

julia> typeof(m.a)
Float64

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float32

На практике такие объекты всегда ведут себя так же, как объекты типа MyStillAmbiguousType.

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

func(m::MyType) = m.a+1

с помощью

code_llvm(func, Tuple{MyType{Float64}})
code_llvm(func, Tuple{MyType{AbstractFloat}})

По причине большой длины результаты здесь не приводятся, но вы можете попробовать сами. Так как в первом случае тип указан полностью, компилятору не нужно создавать код для разрешения типа во время выполнения. Это делает код лаконичнее и быстрее.

Кроме того, следует учитывать, что не полностью параметризированные типы ведут себя как абстрактные. Например, хотя полностью указанный тип Array{T,n} является конкретным, сам по себе тип Array без параметров не является таковым:

julia> !isconcretetype(Array), !isabstracttype(Array), isstructtype(Array), !isconcretetype(Array{Int}), isconcretetype(Array{Int,1})
(true, true, true, true, true)

В данном случае лучше было бы не объявлять MyType с полем a::Array, а вместо этого объявить поле как a::Array{T,N} или как a::A, где {T,N} или A — это параметры MyType.

Избегайте полей с абстрактными контейнерами

Те же рекомендации относятся и к типам контейнеров:

julia> struct MySimpleContainer{A<:AbstractVector}
           a::A
       end

julia> struct MyAmbiguousContainer{T}
           a::AbstractVector{T}
       end

julia> struct MyAlsoAmbiguousContainer
           a::Array
       end

Пример:

julia> c = MySimpleContainer(1:3);

julia> typeof(c)
MySimpleContainer{UnitRange{Int64}}

julia> c = MySimpleContainer([1:3;]);

julia> typeof(c)
MySimpleContainer{Vector{Int64}}

julia> b = MyAmbiguousContainer(1:3);

julia> typeof(b)
MyAmbiguousContainer{Int64}

julia> b = MyAmbiguousContainer([1:3;]);

julia> typeof(b)
MyAmbiguousContainer{Int64}

julia> d = MyAlsoAmbiguousContainer(1:3);

julia> typeof(d), typeof(d.a)
(MyAlsoAmbiguousContainer, Vector{Int64})

julia> d = MyAlsoAmbiguousContainer(1:1.0:3);

julia> typeof(d), typeof(d.a)
(MyAlsoAmbiguousContainer, Vector{Float64})

Для MySimpleContainer объект полностью определяется его типом и параметрами, поэтому компилятор может генерировать оптимизированные функции. В большинстве случаев этого, скорее всего, будет достаточно.

Хотя теперь компилятор отлично справляется со своей задачей, бывают случаи, когда вам может потребоваться реализовать разные операции в коде в зависимости от типа элементов a. Обычно для этого лучше всего заключить определенную операцию (в данном случае foo) в отдельную функцию:

julia> function sumfoo(c::MySimpleContainer)
           s = 0
           for x in c.a
               s += foo(x)
           end
           s
       end
sumfoo (generic function with 1 method)

julia> foo(x::Integer) = x
foo (generic function with 1 method)

julia> foo(x::AbstractFloat) = round(x)
foo (generic function with 2 methods)

Это упрощает код и позволяет компилятору всегда генерировать оптимизированный код.

Однако иногда вам может потребоваться объявить разные версии внешней функции для разных типов элементов или типов AbstractVector поля a в MySimpleContainer. Это можно сделать так:

julia> function myfunc(c::MySimpleContainer{<:AbstractArray{<:Integer}})
           return c.a[1]+1
       end
myfunc (generic function with 1 method)

julia> function myfunc(c::MySimpleContainer{<:AbstractArray{<:AbstractFloat}})
           return c.a[1]+2
       end
myfunc (generic function with 2 methods)

julia> function myfunc(c::MySimpleContainer{Vector{T}}) where T <: Integer
           return c.a[1]+3
       end
myfunc (generic function with 3 methods)
julia> myfunc(MySimpleContainer(1:3))
2

julia> myfunc(MySimpleContainer(1.0:3))
3.0

julia> myfunc(MySimpleContainer([1:3;]))
4

Аннотируйте значения, поступающие из нетипизированных объектов

Зачастую бывает удобно работать со структурами данных, которые могут содержать значения любых типов (массивами типа Array{Any}). Однако если при использовании такой структуры вам точно известен тип элемента, будет полезно сообщить его компилятору:

function foo(a::Array{Any,1})
    x = a[1]::Int32
    b = x+1
    ...
end

Здесь мы знаем, что первый элемент массива a будет иметь тип Int32. Дополнительное преимущество такой аннотации заключается в том, что если значение будет иметь отличный от требуемого тип, во время выполнения произойдет ошибка. Возможно, это позволит выявить другие ошибки ранее в коде.

Если тип элемента a[1] точно не известен, x можно объявить так: x = convert(Int32, a[1])::Int32. Благодаря использованию функции convert элемент a[1] может быть любым объектом, преобразуемым в Int32 (например, UInt8). Это повышает универсальность кода благодаря менее строгим требованиям к типам. Обратите внимание: в этом контексте функция convert сама требует аннотации типа для достижения стабильности типов. Причина в том, что компилятор не может вывести тип возвращаемого значения функции, даже convert, если неизвестны типы каких-либо ее аргументов.

Аннотация типа не повысит (и может даже снизить) производительность, если тип является абстрактным или конструируется во время выполнения. Это связано с тем, что компилятор не может использовать аннотацию для специализации последующего кода, а проверка типов сама по себе занимает время. Например, в следующем коде:

function nr(a, prec)
    ctype = prec == 32 ? Float32 : Float64
    b = Complex{ctype}(a)
    c = (b + 1.0f0)::Complex{ctype}
    abs(c)
end

аннотация c негативно сказывается на производительности. Для написания высокопроизводительного кода с типами, конструируемыми во время выполнения, используйте описываемый ниже прием барьеров между функциями и включите конструируемый тип в типы аргументов функции ядра, чтобы выполняемые ядром операции специализировались компилятором надлежащим образом. Например, в приведенном ниже фрагменте кода, как только создается переменная b, ее можно передать в другую функцию k, то есть ядро. Например, если в функции k аргумент b объявлен как имеющий тип Complex{T}, где T — это параметр типа, то аннотация типа в операторе присваивания внутри k следующего вида:

c = (b + 1.0f0)::Complex{T}

не снижает производительность (но и не повышает ее), так как компилятор может определить тип c во время компиляции функции k.

Учитывайте то, когда специализация в Julia не производится автоматически

Как правило, параметры типов аргументов не специализируются автоматически в Julia в трех конкретных случаях: Type, Function и Vararg. Специализация всегда производится в Julia, когда аргумент используется внутри метода, а не просто передается далее в другую функцию. Обычно это не влияет на производительность во время выполнения и повышает производительность компилятора. Если в вашем случае оказывается, что влияние на производительность есть, вы можете инициировать специализацию, добавив параметр типа в объявление метода. Вот ряд примеров.

В данном случае специализация не производится:

function f_type(t)  # или t::Type
    x = ones(t, 10)
    return sum(map(sin, x))
end

А в данном производится:

function g_type(t::Type{T}) where T
    x = ones(T, 10)
    return sum(map(sin, x))
end

В данных случаях специализация не производится:

f_func(f, num) = ntuple(f, div(num, 2))
g_func(g::Function, num) = ntuple(g, div(num, 2))

А в данном производится:

h_func(h::H, num) where {H} = ntuple(h, div(num, 2))

В данном случае специализация не производится:

f_vararg(x::Int...) = tuple(x...)

А в данном производится:

g_vararg(x::Vararg{Int, N}) where {N} = tuple(x...)

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

h_vararg(x::Vararg{Any, N}) where {N} = tuple(x...)

Обратите внимание, что макрос @code_typed и другие подобные инструменты всегда показывают специализированный код, даже если в Julia данный вызов метода обычно не специализируется. Если вы хотите узнать, создаются ли специализации при изменении типов аргументов, то есть содержит ли (@which f(...)).specializations специализации для соответствующего аргумента, нужно изучить внутреннее устройство метода.

Разделение функций на несколько определений

Если функция реализована в виде нескольких небольших определений, компилятор может напрямую вызвать наиболее подходящий код или даже встроить его.

Вот пример «составной функции», которую на самом деле следует реализовать в виде нескольких определений:

using LinearAlgebra

function mynorm(A)
    if isa(A, Vector)
        return sqrt(real(dot(A,A)))
    elseif isa(A, Matrix)
        return maximum(svdvals(A))
    else
        error("mynorm: invalid argument")
    end
end

Более лаконичным и эффективным будет такой вариант:

norm(x::Vector) = sqrt(real(dot(x, x)))
norm(A::Matrix) = maximum(svdvals(A))

Однако следует отметить, что компилятор весьма эффективно справляется с оптимизацией тупиковых ветвей в коде, как в примере mynorm.

Пишите функции, стабильные с точки зрения типов

Желательно, чтобы функция по возможности всегда возвращала значение одного и того же типа. Рассмотрим следующее определение:

pos(x) = x < 0 ? 0 : x

Хотя оно кажется вполне безобидным, проблема в том, что 0 — это целое число (типа Int), а x может быть любого типа. Таким образом, в зависимости от значения x, эта функция может возвращать значение одного из двух типов. Такое поведение допустимо, а в некоторых случаях может быть желательно. Однако ситуацию можно легко исправить следующим образом:

pos(x) = x < 0 ? zero(x) : x

Существуют также функция oneunit и более общая функция oftype(x, y), которая возвращает аргумент y, преобразованный в тип x.

Старайтесь не изменять тип переменной

Аналогичная проблема «стабильности типов» актуальна и для переменных, используемых повторно внутри функции:

function foo()
    x = 1
    for i = 1:10
        x /= rand()
    end
    return x
end

Локальная переменная x сначала является целочисленной, но после одной итерации цикла становится числом с плавающей запятой (в результате применения оператора /). Поэтому компилятору труднее оптимизировать тело цикла. Решить эту проблему можно несколькими способами.

  • Инициализируйте x со значением x = 1.0.

  • Явным образом объявите тип x как x::Float64 = 1.

  • Используйте явное преобразование: x = oneunit(Float64).

  • Выполните цикл for i = 2:10 с инициализацией x = 1 / rand() в первой итерации.

Отделяйте функции ядра (то есть используйте барьеры между функциями)

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

julia> function strange_twos(n)
           a = Vector{rand(Bool) ? Int64 : Float64}(undef, n)
           for i = 1:n
               a[i] = 2
           end
           return a
       end;

julia> strange_twos(3)
3-element Vector{Int64}:
 2
 2
 2

Правильно было бы написать ее так:

julia> function fill_twos!(a)
           for i = eachindex(a)
               a[i] = 2
           end
       end;

julia> function strange_twos(n)
           a = Vector{rand(Bool) ? Int64 : Float64}(undef, n)
           fill_twos!(a)
           return a
       end;

julia> strange_twos(3)
3-element Vector{Int64}:
 2
 2
 2

Компилятор Julia специализирует код с учетом типов аргументов на границах функций, поэтому в первоначальной реализации тип a во время выполнения цикла неизвестен (так как выбирается случайным образом). Вторая версия, скорее всего, будет работать быстрее, так как внутренний цикл можно перекомпилировать в рамках fill_twos! для разных типов a.

Вторая форма также является более правильной с точки зрения стиля и облегчает повторное использование кода.

Такой шаблон неоднократно применяется в модуле Base Julia. Например, можно ознакомиться с функциями vcat и hcat в abstractarray.jl, а также с функцией fill!, которую можно было бы использовать вместо написания собственной реализации fill_twos!.

Такие функции, как strange_twos, находят применение при работе с данными неопределенного типа, например, загружаемыми из входного файла, который может содержать целые числа, числа с плавающей запятой, строки или что-то еще.

Типы со значениями в качестве параметров

Допустим, вы хотите создать N-мерный массив с размером 3 по каждой оси. Такой массив можно создать следующим образом:

julia> A = fill(5.0, (3, 3))
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

Данный подход очень эффективен: компилятор может определить, что A — это Array{Float64,2}, так как ему известны тип значения для инициализации (5.0::Float64) и размерность ((3, 3)::NTuple{2,Int}). Из этого следует, что компилятор может сгенерировать очень эффективный код для дальнейшего использования массива A в той же функции.

Но давайте теперь предположим, что нужно написать функцию, которая создает массив 3×3×…​ произвольной размерности. Напрашивается следующая реализация:

julia> function array3(fillval, N)
           fill(fillval, ntuple(d->3, N))
       end
array3 (generic function with 1 method)

julia> array3(5.0, 2)
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

Она работает, но (как вы можете убедиться сами с помощью @code_warntype array3(5.0, 2)) проблема в том, что выходной тип вывести невозможно: аргумент N является значением типа Int, и механизм вывода типов не может предугадать его значение. Это означает, что код, в котором используется результат этой функции, должен быть консервативным: тип должен проверяться при каждом обращении к A. Такой код будет работать очень медленно.

Один из очень эффективных способов решения подобных проблем — использование приема барьеров между функциями. Однако в некоторых случаях может потребоваться полностью устранить нестабильность типов. Одним из подходов может быть передача размерности в качестве параметра, например посредством Val{T}() (см. раздел Типы значения):

julia> function array3(fillval, ::Val{N}) where N
           fill(fillval, ntuple(d->3, Val(N)))
       end
array3 (generic function with 1 method)

julia> array3(5.0, Val(2))
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

В Julia есть специальная версия ntuple, которая принимает экземпляр Val{::Int} в качестве второго параметра; передавая N в качестве параметра-типа, вы сообщаете его «значение» компилятору. Следовательно, эта версия array3 позволяет компилятору предугадывать возвращаемый тип.

Однако использование этого приема может скрывать в себе неожиданные тонкости. Например, вызов array3 из функции наподобие следующей будет бесполезен:

function call_array3(fillval, n)
    A = array3(fillval, Val(n))
end

Здесь возникает все та же проблема: компилятор не может определить n, поэтому ему неизвестен тип Val(n). Неправильное использование Val зачастую может привести к снижению производительности. (Применение приведенного выше шаблона кода целесообразно только при грамотном сочетании Val с приемом барьера между функциями, делающем функцию ядра эффективнее.)

Вот пример правильного использования Val.

function filter3(A::AbstractArray{T,N}) where {T,N}
    kernel = array3(1, Val(N))
    filter(A, kernel)
end

В этом примере N передается как параметр, поэтому его «значение» известно компилятору. По сути, Val(T) работает только в том случае, если T жестко закодировано, является литералом (Val(3)) или уже указано в области типа.

Риски некорректного использования множественной диспетчеризации (продолжение темы типов со значениями в качестве параметров)

Оценив преимущества множественной диспетчеризации, программист может увлечься и начать применять ее во всех ситуациях подряд. Например, он может использовать следующую структуру для хранения информации:

struct Car{Make, Model}
    year::Int
    ...more fields...
end

а затем производить диспетчеризацию объектов: Car{:Honda,:Accord}(year, args...).

Это может быть целесообразно при соблюдении одного из следующих условий.

  • Каждый объект Car требует сложной обработки, которая будет гораздо эффективнее, если Make и Model известны во время компиляции, а общее число разных Make или Model будет не слишком велико.

  • Имеются однородные списки объектов Car одного типа, которые можно хранить в массиве Array{Car{:Honda,:Accord},N}.

Если верно последнее, функцию, обрабатывающую такой однородный массив, можно эффективно специализировать: Тип каждого элемента известен Julia заранее (у всех объектов в контейнере один и тот же конкретный тип), поэтому при компиляции функции Julia может определять подходящие вызовы методов (так что этого не нужно будет делать во время выполнения), тем самым генерируя эффективный код для обработки всего списка.

Если эти условия не соблюдаются, скорее всего, пользы не будет. Более того, обусловленный этим подходом «комбинаторный взрыв типов» будет иметь обратный эффект. Если тип items[i+1] отличается от типа item[i], среде Julia придется определять тип во время выполнения, искать подходящий метод в таблицах методов, определять соответствующий метод (посредством пересечения типов), устанавливать, была ли уже произведена его JIT-компиляция (и производить ее, если это не так), а затем выполнять вызов. По сути, вы даете указание системе типов и механизму JIT-компиляции выполнить в коде действие, аналогичное оператору ветвления или поиску по словарю.

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

Возможно, эффект во время компиляции еще серьезнее эффекта во время выполнения: Julia компилирует специализированные функции для каждого отдельного объекта Car{Make, Model}. Если типов сотни или тысячи, то для каждой функции, принимающей такой объект в качестве параметра (собственной пользовательской функции get_year, универсальной функции push! из модуля Base Julia или любой иной), будут компилироваться сотни или тысячи вариантов. Каждый из них увеличивает размер кэша компилируемого кода, длину внутренних списков методов и т. д. Излишний энтузиазм в отношении значений в качестве параметров может легко привести к исчерпанию ресурсов.

Доступ к массивам в порядке их размещения в памяти (по столбцам)

Многомерные массивы Julia хранятся в памяти по столбцам, то есть размещаются столбец за столбцом. Это можно проверить с помощью функции vec или синтаксической конструкции [:], как показано ниже (обратите внимание, что элементы массива следуют в порядке [1 3 2 4], а не [1 2 3 4]).

julia> x = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> x[:]
4-element Vector{Int64}:
 1
 3
 2
 4

Такое соглашение об упорядочении массивов принято во многих языках, включая Фортран, Matlab, R и другие. Альтернативой размещению в памяти по столбцам является размещение по строкам; такое соглашение принято, помимо прочего, в C и Python (numpy). Порядок размещения массивов в памяти может серьезно сказываться на производительности при их переборе в циклах. Помните простое правило: если массив хранится в памяти по столбцам, быстрее всего меняется первый индекс. Это означает, что перебор будет выполняться быстрее, если индекс, используемый в самом глубоко вложенном цикле, является первым в выражении среза. Имейте в виду, что при индексировании массива с помощью : выполняется неявный цикл с поочередным обращением к каждому элементу в определенном измерении. Например, столбцы извлекаются быстрее, чем строки.

Рассмотрим следующий искусственный пример. Предположим, нам нужно написать функцию, которая принимает объект Vector и возвращает квадратную матрицу Matrix, строки или столбцы которой заполнены копиями входного вектора. Допустим, нам неважно, будут ли заполняться этими копиями строки или столбцы (возможно, соответствующее поведение можно настроить в другой части кода). Реализовать это можно по крайней мере четырьмя способами (помимо рекомендуемого вызова встроенной функции repeat):

function copy_cols(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[:, i] = x
    end
    return out
end

function copy_rows(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[i, :] = x
    end
    return out
end

function copy_col_row(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for col = inds, row = inds
        out[row, col] = x[row]
    end
    return out
end

function copy_row_col(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for row = inds, col = inds
        out[row, col] = x[col]
    end
    return out
end

Теперь мы измерим время выполнения каждой из этих функций, подав на вход вектор размером 10000 на 1:

julia> x = randn(10000);

julia> fmt(f) = println(rpad(string(f)*": ", 14, ' '), @elapsed f(x))

julia> map(fmt, [copy_cols, copy_rows, copy_col_row, copy_row_col]);
copy_cols:    0.331706323
copy_rows:    1.799009911
copy_col_row: 0.415630047
copy_row_col: 1.721531501

Обратите внимание, что функция copy_cols работает гораздо быстрее, чем copy_rows. Это вполне ожидаемо, потому что функция copy_cols учитывает порядок размещения объекта Matrix в памяти по столбцам и заполняет его по одному столбцу. Кроме того, функция copy_col_row работает гораздо быстрее, чем copy_row_col, потому что она следует тому правилу, что первый элемент в выражении среза должен использоваться в самом глубоко вложенном цикле.

Предварительное выделение памяти для выходных данных

Если функция возвращает Array или объект другого сложного типа, ей может потребоваться выделять память. К сожалению, зачастую выделение памяти и обратная операция — сборка мусора — являются серьезными узкими местами.

Иногда избежать необходимости выделять память при каждом вызове функции можно путем предварительного выделения памяти для выходных данных. Рассмотрим простейший пример.

julia> function xinc(x)
           return [x, x+1, x+2]
       end;

julia> function loopinc()
           y = 0
           for i = 1:10^7
               ret = xinc(i)
               y += ret[2]
           end
           return y
       end;

с

julia> function xinc!(ret::AbstractVector{T}, x::T) where T
           ret[1] = x
           ret[2] = x+1
           ret[3] = x+2
           nothing
       end;

julia> function loopinc_prealloc()
           ret = Vector{Int}(undef, 3)
           y = 0
           for i = 1:10^7
               xinc!(ret, i)
               y += ret[2]
           end
           return y
       end;

Результаты замера времени:

julia> @time loopinc()
  0.529894 seconds (40.00 M allocations: 1.490 GiB, 12.14% gc time)
50000015000000

julia> @time loopinc_prealloc()
  0.030850 seconds (6 allocations: 288 bytes)
50000015000000

У предварительного выделения памяти есть и другие преимущества. Например, вызывающая сторона может контролировать выходной тип алгоритма. В приведенном выше примере при желании можно было бы передать SubArray вместо Array.

Злоупотребление предварительным выделением памяти может сделать код менее элегантным, поэтому может потребоваться оценить целесообразность его использования и замерить производительность. Однако для «векторизированных» функций (выполняемых поэлементно) можно использовать удобный синтаксис x .= f.(y) для операций на месте с объединенными циклами без временных массивов (см. раздел, посвященный точечному синтаксису для векторизации функций).

Еще больше точек: Объединение векторизированных операций

В Julia есть специальный точечный синтаксис, который преобразует любую скалярную функцию в «векторизированный» вызов функции, а любой оператор — в «векторизированный» оператор. При этом он обладает особым свойством: вложенные «точечные вызовы» объединяются на синтаксическом уровне в единый цикл без выделения памяти для временных массивов. При использовании .= и аналогичных операторов присваивания результат также может храниться на месте в предварительно размещенном в памяти массиве (см. выше).

В контексте линейной алгебры это означает, что, хотя такие операции, как vector + vector и vector * scalar, определены, вместо них может быть лучше использовать vector .+ vector и vector .* scalar, так как получающиеся циклы можно объединять с соседними вычислениями. Например, рассмотрим две функции:

julia> f(x) = 3x.^2 + 4x + 7x.^3;

julia> fdot(x) = @. 3x^2 + 4x + 7x^3; # эквивалентно 3 .* x.^2 .+ 4 .* x .+ 7 .* x.^3

Функции f и fdot производят одни и те же вычисления. Однако функция fdot (определенная с помощью макроса @.) выполняется применительно к массиву гораздо быстрее:

julia> x = rand(10^6);

julia> @time f(x);
  0.019049 seconds (16 allocations: 45.777 MiB, 18.59% gc time)

julia> @time fdot(x);
  0.002790 seconds (6 allocations: 7.630 MiB)

julia> @time f.(x);
  0.002626 seconds (8 allocations: 7.630 MiB)

То есть fdot(x) работает в десять раз быстрее и выделяет в шесть раз меньше памяти, чем f(x), потому что каждая операция * и + в f(x) выделяет память для нового временного массива и выполняется в отдельном цикле. В данном примере быстродействие функции f.(x) такое же, как у функции fdot(x), но во многих ситуациях удобнее просто добавить немного точек в выражения, чем определять отдельную функцию для каждой векторизированной операции.

Используйте представления для срезов

В Julia выражение «среза» массива, например array[1:5, :], создает копию данных (но не в том случае, если оно используется в левой части присваивания: array[1:5, :] = ... выполняет присваивание на месте соответствующей части массива array). Если со срезом производится множество операций, с точки зрения производительности это может быть выгодно, так как эффективнее работать с небольшой непрерывной копией, чем обращаться по индексам к исходному массиву. С другой стороны, если со срезом производится небольшое количество простых операций, затраты на выделение памяти и копирование могут быть существенными.

Альтернативой является создание «представления» массива, то есть объекта-массива (SubArray), который ссылается на данные исходного массива на месте без создания копии. (При записи в представление также изменяются данные в исходном массиве.) Это можно делать для отдельных срезов, вызывая функцию view, или, что еще проще, для всего выражения или блока кода, указывая @views перед выражением. Пример:

julia> fcopy(x) = sum(x[2:end-1]);

julia> @views fview(x) = sum(x[2:end-1]);

julia> x = rand(10^6);

julia> @time fcopy(x);
  0.003051 seconds (3 allocations: 7.629 MB)

julia> @time fview(x);
  0.001020 seconds (1 allocation: 16 bytes)

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

Копирование данных — это не всегда плохо

Массивы хранятся в памяти в виде непрерывных блоков, что делает возможным векторизацию ЦП и сокращение количества обращений к памяти за счет кэширования. По этим же причинам к массивам рекомендуется обращаться по столбцам (см. выше). Случайные схемы доступа и несвязные представления могут серьезно замедлить операции с массивами из-за непоследовательного доступа к памяти.

Копирование данных, обращение к которым происходит случайным образом, в непрерывный массив перед повторным доступом может привести к существенному приросту скорости, как в примере ниже. В этом случае происходит обращение к элементам матрицы по случайным индексам, а затем их умножение. Копирование в обычные массивы ускоряет умножение, несмотря на расходы на операции копирования и выделения памяти.

julia> using Random

julia> A = randn(3000, 3000);

julia> x = randn(2000);

julia> inds = shuffle(1:3000)[1:2000];

julia> function iterated_neural_network(A, x, depth)
           for _ in 1:depth
               x .= max.(0, A * x)
           end
           argmax(x)
       end

julia> @time iterated_neural_network(view(A, inds, inds), x, 10)
  0.324903 seconds (12 allocations: 157.562 KiB)
1569

julia> @time iterated_neural_network(A[inds, inds], x, 10)
  0.054576 seconds (13 allocations: 30.671 MiB, 13.33% gc time)
1569

При условии, что памяти достаточно, затраты на копирование представления в массив с лихвой компенсируются приростом скорости за счет повторного выполнения умножения в непрерывном массиве.

Используйте StaticArrays.jl для небольших операций с векторами и матрицами фиксированного размера

Если в вычислениях задействовано множество небольших (< 100 элементов) массивов фиксированного размера (который известен перед выполнением), то можно попробовать воспользоваться пакетом StaticArrays.jl. Он позволяет представить такие массивы так, что не потребуется лишних выделений памяти в куче, а компилятор сможет специализировать код с учетом размера массива, например путем полного развертывания векторных операций (с устранением циклов) и хранения элементов в регистрах ЦП.

Например, при двумерных геометрических вычислениях может быть много операций с двухкомпонентными векторами. Благодаря использованию типа SVector из пакета StaticArrays.jl становятся доступны удобная векторная нотация и такие операции, как norm(3v - w) применительно к векторам v и w, а компилятор может развернуть код в минимальное вычисление, эквивалентное @inbounds hypot(3v[1]-w[1], 3v[2]-w[2]).

Избегайте интерполяции строк для ввода-вывода

При записи данных в файл (или на другое устройство ввода-вывода) построение промежуточных строк влечет дополнительные затраты. Вместо:

println(file, "$a $b")

используйте:

println(file, a, " ", b)

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

println(file, "$(f(a))$(f(b))")

с таким вариантом:

println(file, f(a), f(b))

Оптимизируйте сетевой ввод-вывод во время параллельного выполнения

При выполнении удаленной функции в параллельном режиме такой код:

using Distributed

responses = Vector{Any}(undef, nworkers())
@sync begin
    for (idx, pid) in enumerate(workers())
        @async responses[idx] = remotecall_fetch(foo, pid, args...)
    end
end

будет быстрее, чем такой:

using Distributed

refs = Vector{Any}(undef, nworkers())
for (idx, pid) in enumerate(workers())
    refs[idx] = @spawnat pid foo(args...)
end
responses = [fetch(r) for r in refs]

Первый вариант предполагает цикл обращения по сети к каждому рабочему процессу, а во втором происходят два сетевых вызова: первый выполняется макросом @spawnat, а второй — методом fetch (или даже wait). Кроме того, методы fetch и wait выполняются последовательно, что в целом снижает производительность.

Устранение предупреждений о выводе из использования

Устаревшая функция выполняет внутренний поиск для однократного вывода соответствующего предупреждения. Эта дополнительная операция поиска может существенно замедлить работу, поэтому все случаи использования устаревших функций следует скорректировать согласно указаниям в предупреждениях.

Тонкости оптимизации

Есть еще ряд менее важных моментов, которые следует учитывать при работе со сплошными внутренними циклами.

  • Избегайте массивов там, где они не нужны. Например, вместо sum([x,y,z]) используйте x+y+z.

  • Используйте abs2(z) вместо abs(z)^2 для комплексного числа z. В целом старайтесь писать код так, чтобы для комплексных аргументов вместо функции abs использовалась функция abs2.

  • Используйте div(x,y) для деления целых чисел с усечением вместо trunc(x/y), fld(x,y) вместо floor(x/y) и cld(x,y) вместо ceil(x/y).

Аннотации производительности

Иногда можно улучшить оптимизацию, заявив определенные свойства программы.

  • Используйте макрос @inbounds, чтобы отключить проверку границ массивов в выражениях. Но это следует делать, только если вы уверены в результате. При выходе за допустимые пределы индексов могут произойти сбой или повреждение данных.

  • Используйте макрос @fastmath, чтобы разрешить оптимизацию чисел с плавающей запятой, которая осуществляется корректно для вещественных чисел, но приводит к отклонениям для чисел IEEE. Используйте эту возможность с осторожностью, так как возможно изменение числовых результатов. Этот макрос соответствует параметру -ffast-math в clang.

  • Указывайте @simd перед циклами for, чтобы заявить, что итерации не зависят друг от друга и их очередность можно изменять. Обратите внимание, что зачастую Julia может автоматически векторизировать код без макроса @simd; он полезен только в случаях, когда такое преобразование иначе было бы недопустимым, в том числе для обеспечения реассоциативности чисел с плавающей запятой и игнорирования зависимых обращений к памяти (@simd ivdep). Опять-таки будьте очень осторожны при использовании утверждения @simd, так как ошибочное аннотирование цикла с зависимыми операциями может привести к непредвиденным результатам. В частности, имейте в виду, что операция setindex! применительно к некоторым подтипам AbstractArray по своей сути зависит от очередности итераций. Эта функция является экспериментальной и может измениться или исчезнуть в будущих версиях Julia.

Стандартная идиома 1:n для индексирования массива AbstractArray небезопасна, если в массиве применяется нестандартное индексирование, и может привести к ошибке сегментации, если проверка границ отключена. Используйте вместо этого LinearIndices(x) или eachindex(x) (см. также главу Массивы с пользовательскими индексами).

Хотя @simd необходимо указывать непосредственно перед наиболее глубоко вложенным циклом for, макросы @inbounds и @fastmath можно применять либо к отдельным выражениям, либо ко всем выражениям во вложенных блоках кода, например с помощью @inbounds begin или @inbounds for ....

Вот пример одновременного использования аннотаций @inbounds и @simd (в этом случае применяется @noinline, чтоб оптимизатор не нарушал ход тестирования):

@noinline function inner(x, y)
    s = zero(eltype(x))
    for i=eachindex(x)
        @inbounds s += x[i]*y[i]
    end
    return s
end

@noinline function innersimd(x, y)
    s = zero(eltype(x))
    @simd for i = eachindex(x)
        @inbounds s += x[i] * y[i]
    end
    return s
end

function timeit(n, reps)
    x = rand(Float32, n)
    y = rand(Float32, n)
    s = zero(Float64)
    time = @elapsed for j in 1:reps
        s += inner(x, y)
    end
    println("GFlop/sec        = ", 2n*reps / time*1E-9)
    time = @elapsed for j in 1:reps
        s += innersimd(x, y)
    end
    println("GFlop/sec (SIMD) = ", 2n*reps / time*1E-9)
end

timeit(1000, 1000)

На компьютере с процессором Intel Core i5 быстродействием 2,4 ГГц результат будет следующим:

GFlop/sec        = 1.9467069505224963
GFlop/sec (SIMD) = 17.578554163920018

(Производительность измеряется в GFlop/sec, и чем больше показатель, тем лучше.)

Вот пример использования всех трех аннотаций. Эта программа сначала вычисляет конечную разность одномерного массива, а затем L2-норму для результата:

function init!(u::Vector)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds @simd for i in 1:n #утверждая, что `u` — это `Vector`, можно предположить, что индексирование начинается с 1
        u[i] = sin(2pi*dx*i)
    end
end

function deriv!(u::Vector, du)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds du[1] = (u[2] - u[1]) / dx
    @fastmath @inbounds @simd for i in 2:n-1
        du[i] = (u[i+1] - u[i-1]) / (2*dx)
    end
    @fastmath @inbounds du[n] = (u[n] - u[n-1]) / dx
end

function mynorm(u::Vector)
    n = length(u)
    T = eltype(u)
    s = zero(T)
    @fastmath @inbounds @simd for i in 1:n
        s += u[i]^2
    end
    @fastmath @inbounds return sqrt(s)
end

function main()
    n = 2000
    u = Vector{Float64}(undef, n)
    init!(u)
    du = similar(u)

    deriv!(u, du)
    nu = mynorm(du)

    @time for i in 1:10^6
        deriv!(u, du)
        nu = mynorm(du)
    end

    println(nu)
end

main()

На компьютере с процессором Intel Core i7 быстродействием 2,7 ГГц результат будет следующим:

$ julia wave.jl;
  1.207814709 seconds
4.443986180758249

$ julia --math-mode=ieee wave.jl;
  4.487083643 seconds
4.443986180758249

Здесь параметр --math-mode=ieee отключает макрос @fastmath, чтобы можно было сравнить результаты.

В данном случае благодаря использованию @fastmath скорость увеличивается примерно в 3,7 раза. Такой прирост необычен — как правило, он бывает меньше. (В данном конкретном примере рабочий набор тестовых данных достаточно мал и полностью помещается в кэше L1 процессора, поэтому задержки при доступе к памяти не имеют значения и время, затрачиваемое на вычисление, в первую очередь определяется загрузкой ЦП. В большинстве реальных программ ситуация иная.) Кроме того, в данном случае эта оптимизация не влияет на результат, но обычно результат получается немного иным. В ряде случаев, особенно при использовании численно неустойчивых алгоритмов, результат может отличаться очень сильно.

Аннотация @fastmath переупорядочивает выражения с плавающей запятой, например, меняя порядок вычислений или предполагая, что некоторые особые случаи (inf, nan) невозможны. В таком случае (на данном конкретном компьютере) главное отличие в том, что выражение 1 / (2*dx) в функции deriv выводится из цикла (вычисляется за его пределами), как если бы мы записали idx = 1 / (2*dx). Внутри цикла выражение ... / (2*dx) преобразуется в ... * idx, которое вычисляется гораздо быстрее. Безусловно, как фактическая оптимизация, применяемая компилятором, так и итоговый прирост скорости во многом зависят от оборудования. Узнать, как изменился сгенерированный код, можно с помощью функции Julia code_native.

Обратите внимание, что макрос @fastmath также предполагает отсутствие значений NaN во время вычисления, которые могут привести к непредвиденному поведению:

julia> f(x) = isnan(x);

julia> f(NaN)
true

julia> f_fast(x) = @fastmath isnan(x);

julia> f_fast(NaN)
false

Обрабатывайте поднормальные числа как нули

Поднормальные числа, ранее называвшиеся денормализованными, полезны во многих ситуациях, но на некотором оборудовании влекут снижение производительности. Вызов set_zero_subnormals(true) разрешает операциям с плавающей запятой интерпретировать поднормальные входные или выходные значения как нулевые, что может улучшить производительность на некотором оборудовании. Вызов set_zero_subnormals(false) строго предписывает обрабатывать поднормальные значения согласно стандарту IEEE.

Ниже приведен пример существенного влияния поднормальных чисел на производительность на некотором оборудовании.

function timestep(b::Vector{T}, a::Vector{T}, Δt::T) where T
    @assert length(a)==length(b)
    n = length(b)
    b[1] = 1                            # Предельное условие
    for i=2:n-1
        b[i] = a[i] + (a[i-1] - T(2)*a[i] + a[i+1]) * Δt
    end
    b[n] = 0                            # Предельное условие
end

function heatflow(a::Vector{T}, nstep::Integer) where T
    b = similar(a)
    for t=1:div(nstep,2)                # Предполагаем, что nstep — четное число
        timestep(b,a,T(0.1))
        timestep(a,b,T(0.1))
    end
end

heatflow(zeros(Float32,10),2)           # Принудительная компиляция
for trial=1:6
    a = zeros(Float32,1000)
    set_zero_subnormals(iseven(trial))  # Нечетные попытки вычисляются в строгом соответствии с IEEE
    @time heatflow(a,1000)
end

Результат будет примерно таким:

  0.002202 seconds (1 allocation: 4.063 KiB)
  0.001502 seconds (1 allocation: 4.063 KiB)
  0.002139 seconds (1 allocation: 4.063 KiB)
  0.001454 seconds (1 allocation: 4.063 KiB)
  0.002115 seconds (1 allocation: 4.063 KiB)
  0.001455 seconds (1 allocation: 4.063 KiB)

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

В этом примере генерируется много поднормальных чисел, потому что значения в a образуют экспоненциально убывающую кривую, которая постепенно выравнивается.

Обработку поднормальных чисел как нулей следует использовать с осторожностью, потому что она может нарушить работу некоторых выражений, например, x-y == 0 подразумевает x == y:

julia> x = 3f-38; y = 2f-38;

julia> set_zero_subnormals(true); (x - y, x == y)
(0.0f0, false)

julia> set_zero_subnormals(false); (x - y, x == y)
(1.0000001f-38, false)

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

a = rand(Float32,1000) * 1.f-9

@code_warntype

Макрос @code_warntype (или его вариант в виде функции code_warntype) может иногда быть полезен для диагностики проблем с типами. Вот пример:

julia> @noinline pos(x) = x < 0 ? 0 : x;

julia> function f(x)
           y = pos(x)
           return sin(y*x + 1)
       end;

julia> @code_warntype f(3.2)
MethodInstance for f(::Float64)
  from f(x) @ Main REPL[9]:1
Arguments
  #self#::Core.Const(f)
  x::Float64
Locals
  y::Union{Float64, Int64}
Body::Float64
1 ─      (y = Main.pos(x))
│   %2 = (y * x)::Float64
│   %3 = (%2 + 1)::Float64
│   %4 = Main.sin(%3)::Float64
└──      return %4

Интерпретация выходных данных макроса @code_warntype, а также родственных ему макросов @code_lowered, @code_typed, @code_llvm и @code_native требует некоторого опыта. Ваш код представляется в форме, которая была существенно переработана в процессе генерирования скомпилированного машинного кода. Большинство выражений имеют аннотации типов в виде ::T (где T может быть, например, Float64). Самая важная особенность макроса @code_warntype заключается в том, что неконкретные типы отображаются красным цветом. Так как данный документ составлен в формате Markdown, не предполагающем выделения цветом, красный текст в нем обозначается прописными буквами.

В первой строке отображается выводимый возвращаемый тип функции — Body::Float64. В следующих строках представлено тело функции f в формате промежуточного представления SSA Julia. Нумерованные поля — это метки, представляющие цели переходов (посредством goto) в коде. Если посмотреть на тело функции, то можно увидеть, что сначала вызывается pos и выводится возвращаемый тип Union Union{Float64, Int64}. Он отображается прописными буквами, так как не является конкретным. Это означает, что точно определить возвращаемый тип pos исходя из входных типов невозможно. Однако результат операции y*x имеет тип Float64 независимо от того, имеет ли y тип Float64 или Int64. В итоге результат f(x::Float64) не будет нестабильным по типу, даже если некоторые промежуточные вычисления являются таковыми.

То, как вы распорядитесь этой информацией, решать вам. Очевидно, лучше всего будет исправить pos так, чтобы обеспечить стабильность типов. В таком случае все переменные в f будут конкретными и производительность будет оптимальной. Однако в определенных ситуациях такая временная нестабильность типов не имеет особого значения. Например, если pos никогда не используется изолированно, тот факт, что выходной тип f стабильный (для входных значений типа Float64), защитит дальнейший код от последствий нестабильности типов. Это особенно важно в случаях, когда устранить нестабильность типов трудно или невозможно. В таких ситуациях приведенные выше советы (например, касательно добавления аннотаций типов или разделения функций) — самое лучшее средство для сдерживания последствий нестабильности типов. Также имейте в виду, что даже в модуле Base Julia есть функции с нестабильностью типов. Например, функция findfirst возвращает индекс позиции в массиве, в которой найден ключ, или nothing, если ключ не найден — явный случай нестабильности типов. Чтобы было проще выявлять случаи нестабильности типов, которые могут иметь значение, объединения Union, содержащие missing или nothing, выделяются желтым цветом, а не красным.

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

  • Тело функции, начинающееся с Body::Union{T1,T2})

    • Интерпретация: функция с нестабильным возвращаемым типом

    • Рекомендация: сделайте возвращаемый тип стабильным, даже если для этого его нужно аннотировать.

  • invoke Main.g(%%x::Int64)::Union{Float64, Int64}

    • Интерпретация: вызов функции g, нестабильной по типу.

    • Рекомендация: исправьте функцию или при необходимости проаннотируйте возвращаемое значение.

  • invoke Base.getindex(%%x::Array{Any,1}, 1::Int64)::Any

    • Интерпретация: доступ к элементам плохо типизированного массива

    • Рекомендация: используйте массивы с более четко определенными типами или при необходимости проаннотируйте типы отдельных обращений к элементам.

  • Base.getfield(%%x, :(:data))::Array{Float64,N} where N

    • Интерпретация: получение поля неконечного типа. В данном случае тип x, допустим ArrayContainer, имеет поле data::Array{T}. Однако для того, чтобы тип Array был конкретным, также требуется измерение N.

    • Рекомендация: используйте конкретные типы, например Array{T,3} или Array{T,N}, где N — параметр ArrayContainer.

Производительность захваченных переменных

Рассмотрим следующий пример, в котором определяется внутренняя функция.

function abmult(r::Int)
    if r < 0
        r = -r
    end
    f = x -> x * r
    return f
end

Функция abmult возвращает функцию f, которая умножает свой аргумент на абсолютное значение r. Внутренняя функция, присваиваемая f, называется замыканием. Внутренние функции также применяются в языке для блоков do и генераторных выражений.

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

В предыдущем абзаце упоминался анализатор. Под этим имеется в виду этап компиляции, на котором впервые загружается модуль, содержащий abmult, в отличие от более позднего этапа, когда он впервые вызывается. Анализатору неизвестно, что Int — это фиксированный тип или что оператор r = -r преобразовывает значение Int в другое значение Int. Магия вывода типов происходит на более позднем этапе компиляции.

Таким образом, анализатору неизвестно, что r имеет фиксированный тип (Int) или что значение r не меняется после создания внутренней функции (и что, следовательно, контейнер не нужен). Поэтому анализатор создает код для контейнера, содержащего объект абстрактного типа, например Any, из-за чего для каждого экземпляра r во время выполнения требуется диспетчеризация типов. Это можно проверить, применив макрос @code_warntype к приведенной выше функции. Как упаковка, так и диспетчеризация типов во время выполнения могут снижать производительность.

Если захваченные переменные используются в части кода, для которой производительность крайне важна, то можно воспользоваться приведенными ниже советами. Во-первых, если известно, что тип захваченной переменной не меняется, его можно объявить явным образом с помощью аннотации типа (применительно к самой переменной, а не в правой части):

function abmult2(r0::Int)
    r::Int = r0
    if r < 0
        r = -r
    end
    f = x -> x * r
    return f
end

Аннотация типа частично компенсирует снижение производительности вследствие захвата, так как анализатор может связать конкретный тип с объектом в контейнере. Далее, если захваченную переменную вовсе не требуется упаковывать (так как она не будет переприсваиваться после создания замыкания), на это можно указать с помощью блока let следующим образом.

function abmult3(r::Int)
    if r < 0
        r = -r
    end
    f = let r = r
            x -> x * r
    end
    return f
end

Блок let создает новую переменную r, область которой ограничена внутренней функцией. Второй прием позволяет полностью восстановить производительность при наличии захваченных переменных. Имейте в виду: этот аспект компилятора быстро развивается, и вполне возможно, что в будущих версиях не потребуется так подробно аннотировать типы для обеспечения производительности. Пока же можно автоматизировать добавление операторов let, как в abmult3, с помощью пакетов от участников сообщества, таких как FastClosures.