Несмотря на то, что в языке Julia по замыслу отсутствует «классическое» объектно-ориентированное программирование с классами и методами, язык предоставляет средства абстрагирования, ключевую роль в которых играет система типов и элементы функционального программирования. Рассмотрим подробнее второй пункт.

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

Уровень 1. Функции как подпрограммы

Выделение подпрограмм и присваивание им собственных имён идёт ещё с доисторических времён, когда Фортран считался языком высокого уровня, а Си ещё не было.

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

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

"""
    sum_all(collection)

Sum all elements of a collection and return the result
"""
function sum_all(collection)
    sum = 0
    for item in collection
        sum += collection
    end
    sum
end

Отличительностью синтаксиса является унаследованное от Лиспа поведение: для «нормального» возврата значения из функции слово return не обязательно: возвращается значение последнего вычисленного перед end выражения. В примере выше будет возвращено значение переменной sum. Таким образом, return может использоваться в качестве маркера особого поведения функции:

function safe_division(number, divisor)
    if divisor == 0
        return 0
    end
    number / divisor
end

# то же самое
function safe_division1(number, divisor)
    if divisor == 0
        0 # при попадании в эту ветку до выхода из функции нет других инструкций
    else
        number / divisor
    end
end

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

hypotenuse(a, b) = sqrt(a^2 + b^2)

«Безопасное» деление с помощью тернарного оператора можно записать так:

safe_division(number, divisor) = divisor == 0 ? 0 : number / divisor

Как можно заметить, для аргументов функции не обязательно указывать типы. С учётом того, как работает JIT-компилятор Julia, «утиная типизация» даже не всегда будет приводить к падению производительности.

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

function safe_division(number, divisor)
    if divisor == 0
        return zero(number / divisor)
    end
    number / divisor
end

Теперь, если типы обоих аргументов известны на этапе компиляции, тип возвращаемого результата также однозначно выводится, т.к. функция zero(x) возвращает нулевое значение того же типа, что её аргумент (а деление на ноль, согласно IEEE 754, имеет вполне представимое значение в формате чисел с плавающей точкой).

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

# фиксированный список аргументов
function hello(name)
    println("Hello, ", name)
end

# аргумент со значением по умолчанию
# аргументы со значением по умолчанию должны идти после обязательных аргументов
function greeting_d(name, greeting = "Hello")
    println(greeting, ", ", name)
end

# именованный аргумент со значением по умолчанию
# именованные аргументы отделяются точкой с запятой и идут после
# обязательных аргументов и опциональных позиционных аргументов
function greeting_kw(name; greeting = "Hello")
    println(greeting, ", ", name)
end

# аргумент greeting хоть и именованный, но должен обязательно быть предоставлен при вызове
function greeting_oblkw(name; greeting)
    println(greeting, ", ", name)
end

# функция с неопределённым числом аргументов
# все аргументы, начиная со второго, в теле интерпретируются как кортеж names
function greeting_nary(greeting, names...)
    print(greeting)
    for name in names
        print(", ", name)
    end
    print('\n')
end

julia> hello("world")
Hello, world

julia> greeting_d("world")
Hello, world

julia> greeting_d("Mr. Smith", "How do you do")
How do you do, Mr. Smith

julia> greeting_kw("Mr. Smith")
Hello, Mr. Smith

julia> greeting_kw("mom", greeting = "Hi")
Hi, mom

julia> greeting_oblkw("world")
ERROR: UndefKeywordError: keyword argument greeting not assigned
Stacktrace:
 [1] greeting_oblkw(::String) at ./REPL[23]:3
 [2] top-level scope at none:0

julia> greeting_oblkw("mom", greeting = "Hi")
Hi, mom

julia> greeting_nary("Hi", "mom", "dad", "everyone")
Hi, mom, dad, everyone

Уровень 2. Функции как данные

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

function f_x_x(fn, x)
    fn(x, x)
end

julia> f_x_x(+, 3)
6 # +(3, 3) = 3+3 = 6
julia> f_x_x(*, 3)
9 # *(3, 3) = 9
julia> f_x_x(^, 3)
27 # ^(3, 3) = 3^3 = 27
julia> f_x_x(log, 3)
1.0 # log(3, 3) = 1

«Классические» функции, которые принимают функциональный аргумент, — это map, reduce и filter.

map(f, x...) применяет функцию f к значениям всех элементов из x(или кортежам из i-х элементов) и возвращает результаты как новую коллекцию:

julia> map(cos, [0, π/3, π/2, 2*π/3, π])
5-element Array{Float64,1}:
  1.0                  
  0.5000000000000001   
  6.123233995736766e-17
 -0.4999999999999998   
 -1.0

julia> map(+, (2, 3), (1, 1))
(3, 4)

reduce(f, x; init_val) «сводит» коллекцию к одному значению, «разворачивая» цепочку f(f(...f(f(init_val, x[1]), x[2])...), x[end]):

function myreduce(fn, values, init_val)
    accum = init_val
    for x in values
        accum = fn(accum, x)
    end
    accum
end

Поскольку на самом деле не определено, в каком порядке будет проход массива при редукции, а также то, будет вызываться fn(accum, x) или fn(x, accum), редукция будет давать предсказуемый результат только с коммутативными или ассоциативными операторами, типа сложения или умножения.

filter(predicate, x) возвращает массив из элементов x, которые удовлетворяют предикату predicate:

julia> filter(isodd, 1:10)
5-element Array{Int64,1}:
 1
 3
 5
 7
 9
julia> filter(iszero, [[0], 1, 0.0, 1:-1, 0im])
4-element Array{Any,1}:
    [0] 
   0.0  
    1:0 
 0 + 0im

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

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

Уровень 3. Функции как абстракции

Нередко в map() или filter() нужно использовать функцию, которой не было присвоено собственное имя. Julia в этом случае позволяет выразить абстракцию выполнения операций над аргументом, не вводя для этой последовательности собственное имя. Такая абстракция называется анонимной функцией, или лямбда-функцией (т.к. в математической традиции такие функции обозначаются буквой лямбда). Синтаксис этого представления такой:

# именованная функция
square(x) = x^2
# анонимный аналог
x -> x^2

# именованная функция
hypot(a, b) = sqrt(x^2 + y^2)
# анонимный аналог - если аргументов больше одного, вокруг них ставятся скобки,
# чтобы не спутать с кортежем из переменной и анонимной функции от одной переменной
(x, y) -> sqrt(x^2 + y^2)

# именованная функция
fortytwo() = 42
# анонимная функция
() -> 42

julia> map(i -> map(x -> x^i, 1:5), 1:5)
5-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4, 5]         
 [1, 4, 9, 16, 25]       
 [1, 8, 27, 64, 125]     
 [1, 16, 81, 256, 625]   
 [1, 32, 243, 1024, 3125]

Как именованные, так и анонимные функции можно присваивать переменным и возвращать в качестве значений:

julia> double_squared = x -> (2 * x)^2
#17 (generic function with 1 method)

julia> double_squared(5)
100

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

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

function normal(x, y)
    z = x + y
    x + y * z
end

function strange(x, y)
    x + y * z
end

Про функцию normal() можно сказать, что в её теле все имена переменных связанные, т.е. если мы везде (включая список аргументов) заменим «x» на «m» (или любой другой идентификатор), «y» на «n», а «z» на «sum_of_m_and_n» — смысл выражения не поменяется. В функции strange() имя z является несвязанным, т.е. а) смысл может поменяться, если это имя будет заменено на другое и б) корректность функции зависит от того, была ли на момент вызова функции определена переменная с именем «z».

Вообще говоря, с функцией normal() тоже не всё так чисто:

Что будет, если вне функции определена переменная с именем z?
Символы + и *, вообще-то, тоже являются несвязанными идентификаторами.

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

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

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

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

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

function f_with_counter(fn)
    call_count = 0
    ncalls() = call_count
    # invoke() будет подсчитывать, сколько раз она вызвана
    # получить это значение можно, вызвав ncalls()
    function invoke(args...)
        call_count += 1
        fn(args...)
    end
    # возвращаем обе функции как кортеж с именованными полями
    # call_count теперь будет извне недоступна напрямую,
    # но экспортированные invoke() и call_count() имеют к нему доступ и могут менять
    (call = invoke, call_count = ncalls)
end

julia> abscount = f_with_counter(abs)
(call = getfield(Main, Symbol("#invoke#22")){typeof(abs)}(abs, Core.Box(0)), call_count = getfield(Main, Symbol("#ncalls#21"))(Core.Box(0)))
julia> abscount.call_count()
0
julia> abscount.call(-20)
20
julia> abscount.call_count()
1
julia> abscount.call(im)
1.0
julia> abscount.call_count()
2

Case study: всё те же полиномы

В прошлой статье рассмотрено представление полиномов как структур. В частности, одна из структур хранения — это список коэффициентов, начиная с младшего. Для вычисления полинома p в точке x предлагалось вызвать функцию evpoly(p, x), которая вычисляет полином по схеме Горнера.

Полный код определений

abstract type AbstractPolynomial end

"""
    Polynomial <: AbstractPolynomial

Polynomials written in the canonical form

---

    Polynomial(v::T) where T<:Union{Vector{<:Real}, NTuple{<:Any, <:Real}})

Construct a `Polynomial` from the list of the coefficients. The coefficients are assumed to go from power 0 in the ascending order. If an empty collection is provided, the constructor returns a zero polynomial.
"""
struct Polynomial<:AbstractPolynomial
    degree::Int
    coeff::NTuple{N, Float64} where N
    function Polynomial(v::T where T<:Union{Vector{<:Real},
                                            NTuple{<:Any, <:Real}})
        coeff = isempty(v) ? (0.0,) : tuple([Float64(x) for x in v]...)
        return new(length(coeff)-1, coeff)
    end
end

"""
    InterpPolynomial <: AbstractPolynomial

Interpolation polynomials in Newton's form

---

    InterpPolynomial(xsample::Vector{<:Real}, fsample::Vector{<:Real})

Construct an `InterpPolynomial` from a vector of points `xsample` and corresponding function values `fsample`. All values in `xsample` must be distinct.
"""
struct InterpPolynomial<:AbstractPolynomial
    degree::Int
    xval::NTuple{N, Float64} where N
    coeff::NTuple{N, Float64} where N
    function InterpPolynomial(xsample::X,
                              fsample::F) where {X<:Union{Vector{<:Real},
                                                          NTuple{<:Any, <:Real}},
                                                 F<:Union{Vector{<:Real},
                                                          NTuple{<:Any, <:Real}}}
        if !allunique(xsample)
            throw(DomainError("Cannot interpolate with duplicate X points"))
        end
        N = length(xsample)
        if length(fsample) != N
            throw(DomainError("Lengths of X and F are not the same"))
        end

        coeff = [Float64(f) for f in fsample]

        for i = 2:N
            for j = 1:(i-1)
                coeff[i] = (coeff[j] - coeff[i]) / (xsample[j] - xsample[i])
            end
        end
        new(N-1, ntuple(i -> Float64(xsample[i]), N), tuple(coeff...))
    end
end

function InterpPolynomial(fn, xsample::T) where {T<:Union{Vector{<:Real},
                                                          NTuple{<:Any, <:Real}}}
    InterpPolynomial(xsample, map(fn, xsample))
end

function evpoly(p::Polynomial, z::Real)
    ans = p.coeff[end]
    for idx = p.degree:-1:1
        ans = p.coeff[idx] + z * ans
    end
    return ans
end

function evpoly(p::InterpPolynomial, z::Real)
    ans = p.coeff[p.degree+1]
    for idx = p.degree:-1:1
        ans = ans * (z - p.xval[idx]) + p.coeff[idx]
    end
    return ans
end

function Base.:+(p1::Polynomial, p2::Polynomial)
    # при сложении нужно учесть, какой должна быть наивысшая степень
    deg = max(p1.degree, p2.degree)
    coeff = zeros(deg+1)

    coeff[1:p1.degree+1] .+= p1.coeff
    coeff[1:p2.degree+1] .+= p2.coeff

    Polynomial(coeff)
end

function Base.:+(p1::InterpPolynomial, p2::InterpPolynomial)
    xmax = max(p1.xval..., p2.xval...)
    xmin = min(p1.xval..., p2.xval...)
    deg = max(p1.degree, p2.degree)
    # для построения суммы строим чебышёвскую сетку между минимальным
    # и максимальным из узлов обоих полиномов
    xmid = 0.5 * xmax + 0.5 * xmin
    dx = 0.5 * (xmax - xmin) / cos(0.5 * π / (deg + 1))
    chebgrid = [xmid + dx * cos((k - 0.5) * π / (deg + 1)) for k = 1:deg+1]
    fsample = [evpoly(p1, x) + evpoly(p2, x) for x in chebgrid]
    InterpPolynomial(chebgrid, fsample)
end

function Base.:+(p1::InterpPolynomial, p2::Polynomial)
    xmax = max(p1.xval...)
    xmin = min(p1.xval...)
    deg = max(p1.degree, p2.degree)
    xmid = 0.5 * xmax + 0.5 * xmin
    dx = 0.5 * (xmax - xmin) / cos(0.5 * π / (deg + 1))
    chebgrid = [xmid + dx * cos((k - 0.5) * π / (deg + 1)) for k = 1:deg+1]
    fsample = [evpoly(p1, x) + evpoly(p2, x) for x in chebgrid]
    InterpPolynomial(chebgrid, fsample)
end

function Base.:+(p1::Polynomial, p2::InterpPolynomial)
    p2 + p1
end

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

struct Polynomial
    degree::Int
    coeff::NTuple{N, Float64} where N
    function Polynomial(v::T where T<:Union{Vector{<:Real},
                                            NTuple{<:Any, <:Real}})
        # в случае пустого массива / кортежа в аргументе возвращаем P(x) ≡ 0
        coeff = isempty(v) ? (0.0,) : tuple([Float64(x) for x in v]...)
        # возврат значения - специальным оператором new
        # аргументы - перечисление значений полей
        return new(length(coeff)-1, coeff)
    end
end

"""
    evpoly(p::Polynomial, z::Real)

Evaluate polynomial `p` at `z` using the Horner's rule
"""
function evpoly(p::Polynomial, z::Real)
    ans = p.coeff[end]
    for idx = p.degree:-1:1
        ans = p.coeff[idx] + z * ans
    end
    return ans
end

Переделаем это определение в функцию, которая принимает массив / кортеж коэффициентов и возвращает собственно функцию, вычисляющую полином:

function Polynomial_as_closure(v::T where T<:Union{Vector{<:Real},
                                        NTuple{<:Any, <:Real}})
    # в случае пустого массива / кортежа в аргументе возвращаем P(x) ≡ 0
    if isempty(v)
        return x::Real -> 0.0
    end
    
    coeff = tuple(map(float, v)...)
    degree = length(coeff) - 1
    
    function evpoly(z::Real)
        ans = coeff[end]
        for idx = degree:-1:1
            ans = coeff[idx] + z * ans
        end
        return ans
    end

    evpoly
end

julia> p = Polynomial_as_closure((0, 1, 1)) # x² + x
(::getfield(Main, Symbol("#evpoly#28")){Tuple{Float64,Float64,Float64},Int64}) (generic function with 1 method)

julia> p(1) # ура, можно обойтись без evpoly()!
2.0
julia> p(11)
132.0

Аналогично можно написать функцию и для интерполяционного полинома.

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

К счастью, в этом случае мы можем взять лучшее из двух миров, поскольку Julia позволяет создавать так называемые callable structs. Т.е. можно задать полином как структуру, но иметь возможность его вызова как функции! К определениям структур из предыдущей статьи нужно лишь добавить:

function (p::Polynomial)(z::Real)
    evpoly(p, z)
end

function (p::InterpPolynomial)(z::Real)
    evpoly(p, z)
end

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

function InterpPolynomial(fn, xsample::T) where {T<:Union{Vector{<:Real},
                                                          NTuple{<:Any, <:Real}}}
    InterpPolynomial(xsample, map(fn, xsample))
end

Проверим определение

julia> psin = InterpPolynomial(sin, [0, π/6, π/2, 5*π/6, π]) # интерполяция синуса
InterpPolynomial(4, (0.0, 0.5235987755982988, 1.5707963267948966, 2.6179938779914944, 3.141592653589793), (0.0, 0.954929658551372, -0.30396355092701327, -0.05805276197975913, 0.036957536116863636))
julia> pcos = InterpPolynomial(cos, [0, π/6, π/2, 5*π/6, π]) # интерполяция косинуса
InterpPolynomial(4, (0.0, 0.5235987755982988, 1.5707963267948966, 2.6179938779914944, 3.141592653589793), (1.0, -0.2558726308373678, -0.36358673785585766, 0.1388799037738005, 5.300924469105863e-17))
julia> psum = pcos + psin 
InterpPolynomial(4, (3.141592653589793, 2.5416018461576297, 1.5707963267948966, 0.5999908074321635, 0.0), (-1.0, -1.2354929267138448, 0.03888175053443867, 0.1969326657535598, 0.03695753611686364))
julia> for x = range(0, π, length = 20)
           println("Error at x = ", x, ": ", abs(psum(x) - (sin(x) + cos(x))))
       end
Error at x = 0.0: 0.0
Error at x = 0.3490658503988659: 0.002748366490382681
Error at x = 0.6981317007977318: 0.0031870524474437723
Error at x = 1.0471975511965976: 0.006538414090220712
Error at x = 1.3962634015954636: 0.0033647273630357244
Error at x = 1.7453292519943295: 0.003570894863996865
Error at x = 2.0943951023931953: 0.007820939854677023
Error at x = 2.443460952792061: 0.004305934583281101
Error at x = 2.792526803190927: 0.00420977797025246
Error at x = 3.141592653589793: 1.1102230246251565e-16

Заключение

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

Julia: функции и структуры-как-функции