Spec-Zone.ru › Julia 1.8

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

В следующих разделах мы кратко рассмотрим несколько техник, которые помогут сделать ваш код 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).

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

Если вместо этого мы передаем 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

Единственное выделение памяти связано с выполнением макроса @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. Однако, если мы вместо этого будем разрешать хранение только чисел одного типа, например, Float64, в a, эти числа можно хранить более эффективно:

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)  # or 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)
  • Инициализировать с первой итерацией цикла, в x = 1 / rand(), затем цикл for i = 2:10

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

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

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.

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

Этот шаблон используется в нескольких местах в Julia Base. Например, см. 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-компиляции выполнять в вашем собственном коде эквивалент инструкции switch или поиска в словаре.

Примеры бенчмарков, сравнивающих (1) диспетчерирование по типам, (2) поиск в словаре и (3) инструкцию «switch», можно найти на форуме рассылки.

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

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

Многомерные массивы в 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

Эта конвенция для упорядочения массивов распространена во многих языках, таких как Fortran, 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; # equivalent to 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) в десять раз быстрее и потребляет в 6 раз меньше памяти, чем 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)

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

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

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

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

julia> using Random

julia> x = randn(1_000_000);

julia> inds = shuffle(1:1_000_000)[1:800000];

julia> A = randn(50, 1_000_000);

julia> xtmp = zeros(800_000);

julia> Atmp = zeros(50, 800_000);

julia> @time sum(view(A, :, inds) * view(x, inds))
  0.412156 seconds (14 allocations: 960 bytes)
-4256.759568345458

julia> @time begin
           copyto!(xtmp, view(x, inds))
           copyto!(Atmp, view(A, :, inds))
           sum(Atmp * xtmp)
       end
  0.285923 seconds (14 allocations: 960 bytes)
-4256.759568345134

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

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

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

Например, если вы выполняете вычисления с 2d геометриями, у вас может быть много вычислений с векторами из 2 компонентов. Используя тип 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. В общем случае старайтесь переписать код, чтобы использовать 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; он полезен только в тех случаях, когда такая трансформация по другим причинам была бы незаконной, включая случаи, когда это позволяет перегруппировать операции с плавающей запятой и игнорировать зависимые обращения к памяти (@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 #by asserting that `u` is a `Vector` we can assume it has 1-based indexing
        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                            # 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

Это даёт вывод, похожий на

  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)
Variables
  #self#::Core.Const(f)
  x::Float64
  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 IR Julia. Номерные блоки являются метками и представляют цели для переходов (через goto) в вашем коде. Рассмотрев тело, вы увидите, что в первую очередь вызывается pos, а возвращаемое значение выведено как тип Union UNION{FLOAT64, INT64}, показанный заглавными буквами, так как это неконкретный тип. Это означает, что мы не можем знать точный возвращаемый тип pos на основе входных типов. Однако результат y*x является Float64, независимо от того, является ли y типом Float64 или Int64. Конечный результат заключается в том, что f(x::Float64) не будет иметь проблем с типом при выводе, даже если некоторые промежуточные вычисления нестабильны по типу.

Как использовать эту информацию — решать вам. Очевидно, что было бы намного лучше исправить pos на стабильный тип: в этом случае все переменные в f были бы конкретными, а его производительность была бы оптимальной. Однако существуют обстоятельства, когда такая временная нестабильность типа может не иметь большого значения: например, если pos никогда не используется изолированно, то тот факт, что вывод f стабилен по типу (для входных данных Float64) защитит последующий код от распространения эффектов нестабильности типа. Это особенно актуально в тех случаях, когда исправление нестабильности типа затруднено или невозможно. В таких случаях указанные выше советы (например, добавление аннотаций типов и/или разделение функций) являются лучшими инструментами для ограничения «ущерба» от нестабильности типов. Также обратите внимание, что даже в Julia Base есть функции, которые нестабильны по типу. Например, функция 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, область действия которой ограничена только внутренней функцией. Второй метод восстанавливает полную производительность языка в присутствии захваченных переменных. Обратите внимание, что этот аспект компилятора быстро развивается, и, вероятно, будущие версии не будут требовать такой степени аннотаций программиста для достижения производительности. Тем временем, некоторые пакеты, созданные пользователями, такие как FastClosures, автоматизируют вставку инструкций let как в abmult3.

© 2009–2022 Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors
Licensed under the MIT License.
https://docs.julialang.org/en/v1.8/manual/performance-tips/

Spec-Zone.ru

Настройки Оффлайн Что нового Помощь О нас
Spec-Zone .ru
спецификации, руководства, описания, API