Рекомендации по производительности

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

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

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

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

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

const DEFAULT_VAL = 0

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

global x
y = f(x::Int + 1)

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

В следующем сеансе REPL:

julia> x = 1.0

эквивалентно:

julia> global x = 1.0

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

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

Полезным инструментом для измерения производительности является макрос @time. Следующий пример иллюстрирует хороший стиль работы:

julia> function f(n)
           s = 0
           for i = 1:n
               s += i/2
           end
           s
       end
f (generic function with 1 method)

julia> @time f(1)
  0.012686 seconds (2.09 k allocations: 103.421 KiB)
0.5

julia> @time f(10^6)
  0.021061 seconds (3.00 M allocations: 45.777 MiB, 11.69% gc time)
2.5000025e11

При первом вызове (@time f(1)), f компилируется. (Если вы ещё не использовали @time в этом сеансе, он также скомпилирует функции, необходимые для измерения времени.) Результаты этого запуска не следует принимать всерьёз. При втором запуске, помимо отчёта о времени, указывается, что выделено большое количество памяти. Это единственное преимущество @time по сравнению с функциями, такими как tic() и toc(), которые только сообщают о времени.

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

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

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

julia> @time f_improved(1)
  0.007008 seconds (1.32 k allocations: 63.640 KiB)
0.5

julia> @time f_improved(10^6)
  0.002997 seconds (6 allocations: 192 bytes)
2.5000025e11

Ниже вы узнаете, как выявить проблему с f и как её исправить.

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

Инструменты

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

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

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

• @code_warntype генерирует представление вашего кода, которое может быть полезно для поиска выражений, которые приводят к неопределённости типов. См. @code_warntype ниже.

• Пакет Lint также может предупредить вас о некоторых типах программистских ошибок.

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

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

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

a = Real[]    # typeof(a) = Array{Real,1}
if (f = rand()) < .8
    push!(a, f)
end

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

a = Float64[] # typeof(a) = Array{Float64,1}

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

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

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

Во многих языках с необязательными объявлениями типов добавление объявлений является основным способом ускорения работы кода. Это не относится к 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}})
code_llvm(func,Tuple{MyType})

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

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

Те же принципы хорошей практики также работают для типов контейнеров:

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

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

Например:

julia> c = MySimpleContainer(1:3);

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

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

julia> typeof(c)
MySimpleContainer{Array{Int64,1}}

julia> b = MyAmbiguousContainer(1:3);

julia> typeof(b)
MyAmbiguousContainer{Int64}

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

julia> typeof(b)
MyAmbiguousContainer{Int64}

Для 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)

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

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

function myfun(c::MySimpleContainer{Vector{T}}) where T<:AbstractFloat
    ...
end
function myfun(c::MySimpleContainer{Vector{T}}) where T<:Integer
    ...
end

Это работает нормально для Vector{T}, но нам также пришлось бы написать явные версии для UnitRange{T} или других абстрактных типов. Чтобы предотвратить подобную рутину, вы можете использовать два параметра в объявлении MyContainer:

julia> mutable struct MyContainer{T, A<:AbstractVector}
           a::A
       end

julia> MyContainer(v::AbstractVector) = MyContainer{eltype(v), typeof(v)}(v)
MyContainer

julia> b = MyContainer(1:5);

julia> typeof(b)
MyContainer{Int64,UnitRange{Int64}}

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

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

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

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

julia> myfunc(MyContainer(1:3))
2

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

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

Как вы можете видеть, с этим подходом можно специализироваться как на типе элемента T, так и на типе массива A.

Однако, есть одна оставшаяся проблема: мы не потребовали, чтобы A имел тип элемента T, поэтому вполне возможно создать такой объект:

julia> b = MyContainer{Int64, UnitRange{Float64}}(UnitRange(1.3, 5.0));

julia> typeof(b)
MyContainer{Int64,UnitRange{Float64}}

Чтобы предотвратить это, можно добавить внутренний конструктор:

julia> mutable struct MyBetterContainer{T<:Real, A<:AbstractVector}
           a::A
           MyBetterContainer{T,A}(v::AbstractVector{T}) where {T,A} = new(v)
       end

julia> MyBetterContainer(v::AbstractVector) = MyBetterContainer{eltype(v),typeof(v)}(v)
MyBetterContainer

julia> b = MyBetterContainer(UnitRange(1.3, 5.0));

julia> typeof(b)
MyBetterContainer{Float64,UnitRange{Float64}}

julia> b = MyBetterContainer{Int64, UnitRange{Float64}}(UnitRange(1.3, 5.0));
ERROR: MethodError: Cannot `convert` an object of type UnitRange{Float64} to an object of type MyBetterContainer{Int64,UnitRange{Float64}}
[...]

Внутренний конструктор требует, чтобы тип элемента A был T.

Аннотирование значений, взятых из нетипизированных расположений

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

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

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

Объявление типов ключевых аргументов

Ключевые аргументы могут иметь объявленные типы:

function with_keyword(x; name::Int = 1)
    ...
end

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

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

Передача динамических списков ключевых аргументов, как в f(x; keywords...), может быть медленной и следует избегать в производительно-чувствительном коде.

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

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

Вот пример "сложной функции", которую следует написать как множество определений:

function norm(A)
    if isa(A, Vector)
        return sqrt(real(dot(A,A)))
    elseif isa(A, Matrix)
        return maximum(svd(A)[2])
    else
        error("norm: invalid argument")
    end
end

Это можно записать более лаконично и эффективно как:

norm(x::Vector) = sqrt(real(dot(x,x)))
norm(A::Matrix) = maximum(svd(A)[2])

Написание "стабильных по типу" функций

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

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

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

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

Также существует функция one() и более общая функция oftype(x, y), которая возвращает y, преобразованное к типу x.

Избегайте изменения типа переменной

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

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

Локальная переменная x начинается как целое число, а после одной итерации цикла становится числом с плавающей точкой (результат оператора /). Это затрудняет оптимизацию тела цикла компилятором. Существует несколько возможных исправлений:

• Инициализировать x значением x = 1.0

• Объявить тип x: x::Float64 = 1

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

• Инициализировать переменную первым значением итерации, т.е. x = 1/bar(), а затем проитерировать for i = 2:10

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

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

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

julia> strange_twos(3)
3-element Array{Float64,1}:
 2.0
 2.0
 2.0

Это следует записать как:

julia> function fill_twos!(a)
           for i=1:length(a)
               a[i] = 2
           end
       end
fill_twos! (generic function with 1 method)

julia> function strange_twos(n)
           a = Array{rand(Bool) ? Int64 : Float64}(n)
           fill_twos!(a)
           return a
       end
strange_twos (generic function with 1 method)

julia> strange_twos(3)
3-element Array{Float64,1}:
 2.0
 2.0
 2.0

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

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

Этот шаблон используется в нескольких местах стандартной библиотеки. Например, см. hvcat_fill в abstractarray.jl или функцию fill!, которую мы могли бы использовать вместо написания собственной fill_twos!.

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

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

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

julia> A = fill(5.0, (3, 3))
3×3 Array{Float64,2}:
 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 Array{Float64,2}:
 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, ::Type{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 Array{Float64,2}:
 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 во время компиляции.

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

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

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

Некоторые эталонные тесты времени выполнения, сравнивающие (1) диспетчерирование типов, (2) поиск в словаре и (3) оператор "switch", можно найти на списке рассылки.

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

Доступ к массивам в порядке памяти, по столбцам

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

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

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

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

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

function copy_cols(x::Vector{T}) where T
    n = size(x, 1)
    out = Array{T}(n, n)
    for i = 1:n
        out[:, i] = x
    end
    out
end

function copy_rows(x::Vector{T}) where T
    n = size(x, 1)
    out = Array{T}(n, n)
    for i = 1:n
        out[i, :] = x
    end
    out
end

function copy_col_row(x::Vector{T}) where T
    n = size(x, 1)
    out = Array{T}(n, n)
    for col = 1:n, row = 1:n
        out[row, col] = x[row]
    end
    out
end

function copy_row_col(x::Vector{T}) where T
    n = size(x, 1)
    out = Array{T}(n, n)
    for row = 1:n, col = 1:n
        out[row, col] = x[col]
    end
    out
end

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

julia> x = randn(10000);

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

julia> map(fmt, Any[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 или какой-либо другой сложный тип, ей может потребоваться выделить память. К сожалению, выделение памяти и обратный процесс — сборка мусора — часто являются существенными узкими местами.

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

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

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

с

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

function loopinc_prealloc()
    ret = Array{Int}(3)
    y = 0
    for i = 1:10^7
        xinc!(ret, i)
        y += ret[2]
    end
    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, так как полученные циклы могут быть объединены с окружающими вычислениями. Например, рассмотрим две функции:

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

fdot(x) = @. 3x^2 + 4x + 7x^3 # equivalent to 3 .* x.^2 .+ 4 .* x .+ 7 .* x.^3

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

julia> x = rand(10^6);

julia> @time f(x);
  0.010986 seconds (18 allocations: 53.406 MiB, 11.45% gc time)

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

julia> @time f.(x);
  0.003297 seconds (30 allocations: 7.631 MiB)

То есть, fdot(x) в три раза быстрее и выделяет в 7 раз меньше памяти, чем 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 (7 allocations: 7.630 MB)

julia> @time fview(x);
  0.001020 seconds (6 allocations: 224 bytes)

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

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

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

println(file, "$a $b")

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

println(file, a, " ", b)

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

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

против:

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

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

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

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

быстрее, чем:

refs = Vector{Any}(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. В общем, старайтесь переписывать код, чтобы использовать abs2() вместо abs() для комплексных аргументов.

• Используйте 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 необходимо размещать непосредственно перед циклом, @inbounds и @fastmath могут быть применены к нескольким операторам сразу, например, с помощью begin ... end, или даже к всей функции.

Вот пример с аннотациями @inbounds и @simd.

function inner(x, y)
    s = zero(eltype(x))
    for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end

function innersimd(x, y)
    s = zero(eltype(x))
    @simd for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    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        = ",2.0*n*reps/time*1E-9)
    time = @elapsed for j in 1:reps
        s+=innersimd(x,y)
    end
    println("GFlop/sec (SIMD) = ",2.0*n*reps/time*1E-9)
end

timeit(1000,1000)

На компьютере с процессором Intel Core i5 2,4 ГГц это даёт:

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

(GFlop/sec измеряет производительность, и более высокие числа лучше.) Диапазон для цикла @simd for должен быть одномерным. Переменная, используемая для накопления, такая как s в примере, называется переменной редукции. Используя @simd, вы утверждаете несколько свойств цикла:

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

• Действия с плавающей точкой над переменными редукции могут быть переупорядочены, что может привести к результатам, отличным от результатов без @simd.

• Ни одна итерация не ожидает другой итерации для продвижения вперед.

Цикл, содержащий break, continue или @goto, вызовет ошибку на этапе компиляции.

Использование @simd просто даёт компилятору разрешение на векторизацию. Будет ли она действительно применена, зависит от компилятора. Чтобы действительно извлечь выгоду из текущей реализации, ваш цикл должен обладать следующими дополнительными свойствами:

• Цикл должен быть внутренним циклом.

• Тело цикла должно быть кодом без ветвлений. Вот почему @inbounds в настоящее время необходим для всех обращений к массивам. Компилятор иногда может преобразовать короткие &&, || и ?: выражения в код без ветвлений, если безопасно вычислять все операнды безусловно. Рассмотрите использование ifelse() вместо ?: в цикле, если это безопасно.

• Доступы должны иметь шаблон шага и не могут быть «gather» (чтение с произвольными индексами) или «scatter» (запись с произвольными индексами).

• Шаг должен быть единичным.

• В некоторых простых случаях, например, при обращении к 2-3 массивам в цикле, автоматическая векторизация LLVM может активироваться, что не приведёт к дальнейшему ускорению с помощью @simd.

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

function init!(u)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds @simd for i in 1:n
        u[i] = sin(2pi*dx*i)
    end
end

function deriv!(u, 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 norm(u)
    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/n)
end

function main()
    n = 2000
    u = Array{Float64}(n)
    init!(u)
    du = similar(u)

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

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

    println(nu)
end

main()

На компьютере с процессором Intel Core i7 с тактовой частотой 2,7 ГГц это даёт:

$ julia wave.jl;
elapsed time: 1.207814709 seconds (0 bytes allocated)

$ julia --math-mode=ieee wave.jl;
elapsed time: 4.487083643 seconds (0 bytes allocated)

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

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

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

Обработка чисел с малой нормой как нулей

Числа с малой нормой, ранее называемые числами с малой нормой, полезны во многих контекстах, но на некоторых аппаратных платформах приводят к снижению производительности. Вызов 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                            # Boundary condition
    for i=2:n-1
        b[i] = a[i] + (a[i-1] - T(2)*a[i] + a[i+1]) * Δt
    end
    b[n] = 0                            # Boundary condition
end

function heatflow(a::Vector{T}, nstep::Integer) where T
    b = similar(a)
    for t=1:div(nstep,2)                # Assume nstep is even
        timestep(b,a,T(0.1))
        timestep(a,b,T(0.1))
    end
end

heatflow(zeros(Float32,10),2)           # Force compilation
for trial=1:6
    a = zeros(Float32,1000)
    set_zero_subnormals(iseven(trial))  # Odd trials use strict IEEE arithmetic
    @time heatflow(a,1000)
end

В этом примере генерируется много чисел с малой нормой, поскольку значения в 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()) иногда может быть полезен для диагностики проблем, связанных с типами. Вот пример:

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

function f(x)
    y = pos(x)
    sin(y*x+1)
end

julia> @code_warntype f(3.2)
Variables:
  #self#::#f
  x::Float64
  y::UNION{FLOAT64,INT64}
  fy::Float64
  #temp#@_5::UNION{FLOAT64,INT64}
  #temp#@_6::Core.MethodInstance
  #temp#@_7::Float64

Body:
  begin
      $(Expr(:inbounds, false))
      # meta: location REPL[1] pos 1
      # meta: location float.jl < 487
      fy::Float64 = (Core.typeassert)((Base.sitofp)(Float64,0)::Float64,Float64)::Float64
      # meta: pop location
      unless (Base.or_int)((Base.lt_float)(x::Float64,fy::Float64)::Bool,(Base.and_int)((Base.and_int)((Base.eq_float)(x::Float64,fy::Float64)::Bool,(Base.lt_float)(fy::Float64,9.223372036854776e18)::Bool)::Bool,(Base.slt_int)((Base.fptosi)(Int64,fy::Float64)::Int64,0)::Bool)::Bool)::Bool goto 9
      #temp#@_5::UNION{FLOAT64,INT64} = 0
      goto 11
      9:
      #temp#@_5::UNION{FLOAT64,INT64} = x::Float64
      11:
      # meta: pop location
      $(Expr(:inbounds, :pop))
      y::UNION{FLOAT64,INT64} = #temp#@_5::UNION{FLOAT64,INT64} # line 3:
      unless (y::UNION{FLOAT64,INT64} isa Int64)::ANY goto 19
      #temp#@_6::Core.MethodInstance = MethodInstance for *(::Int64, ::Float64)
      goto 28
      19:
      unless (y::UNION{FLOAT64,INT64} isa Float64)::ANY goto 23
      #temp#@_6::Core.MethodInstance = MethodInstance for *(::Float64, ::Float64)
      goto 28
      23:
      goto 25
      25:
      #temp#@_7::Float64 = (y::UNION{FLOAT64,INT64} * x::Float64)::Float64
      goto 30
      28:
      #temp#@_7::Float64 = $(Expr(:invoke, :(#temp#@_6), :(Main.*), :(y), :(x)))
      30:
      return $(Expr(:invoke, MethodInstance for sin(::Float64), :(Main.sin), :((Base.add_float)(#temp#@_7,(Base.sitofp)(Float64,1)::Float64)::Float64)))
  end::Float64

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

Верхняя часть вывода обобщает информацию о типе для различных внутренних переменных функции. Вы можете видеть, что y, одна из переменных, созданных вами, является Union{Int64,Float64}, из-за нестабильности типа pos. Существует ещё одна переменная, _var4, тип которой также такой же.

Следующие строки представляют тело f. Строки, начинающиеся с числа, за которым следует двоеточие (1:, 2: ), являются метками и представляют цели переходов (через goto) в вашем коде. Смотря на тело, вы можете видеть, что pos был встроен в f — всё, что перед 2:, происходит из кода, определённого в pos.

Начиная с 2:, определена переменная y и снова помечена как тип Union . Далее мы видим, что компилятор создал временную переменную _var1 для хранения результата y*x . Поскольку Float64 умножается либо на Int64, либо на Float64, результат будет типа Float64, поэтому здесь заканчивается вся нестабильность типа. В итоге f(x::Float64) не будет иметь нестабильности типа в своём выводе, даже если некоторые промежуточные вычисления будут иметь нестабильность типа.

Как вы будете использовать эту информацию, решать вам. Очевидно, было бы наилучшим способом исправить pos на стабильный тип: если вы это сделаете, все переменные в f будут конкретными, и его производительность будет оптимальной. Однако существуют обстоятельства, при которых этот вид эпизодической нестабильности типа может не иметь большого значения: например, если pos никогда не используется изолированно, то тот факт, что вывод f стабилен (для входных данных Float64), защитит последующий код от распространения эффектов нестабильности типа. Это особенно актуально в тех случаях, когда исправление нестабильности типа сложно или невозможно: например, в настоящее время нельзя вывести тип возвращаемого значения анонимной функции. В таких случаях вышеуказанные советы (например, добавление аннотаций типа и/или разделение функций) являются вашими лучшими инструментами для ограничения «вреда» от нестабильности типа.

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

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

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

  • Предложение: сделайте тип возвращаемого значения стабильным, даже если вам придётся это аннотировать

• f(x::T)::Union{T1,T2}

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

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

• (top(arrayref))(A::Array{Any,1},1)::Any

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

  • Предложение: используйте массивы с лучше определёнными типами или, при необходимости, аннотируйте тип отдельных обращений к элементам

• (top(getfield))(A::ArrayContainer{Float64},:data)::Array{Float64,N}

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

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