Julia: tipos, multimétodos e aritmética sobre polinômios

Nesta publicação, focaremos o principal, na minha opinião, característica distintiva da linguagem Julia - a representação de funções na forma de métodos com despacho múltiplo. Isso permite aumentar o desempenho dos cálculos sem reduzir a legibilidade do código e sem prejudicar a abstração, por um lado, e permite trabalhar com conceitos matemáticos em uma notação mais familiar, por outro. Como exemplo, a questão do uniforme (do ponto de vista das operações lineares) trabalha com polinômios na representação da lista de coeficientes e com polinômios de interpolação.

Sintaxe básica

Uma breve introdução para quem não conhece. Julia é uma linguagem semelhante a script, possui REPL (loop de leitura-avaliação-impressão, ou seja, um shell interativo). À primeira vista, parece bastante semelhante, por exemplo, ao Python ou ao MATLAB.

Operações aritméticas

A aritmética é quase a mesma de todos os lugares: +, -, *, /, ^ para exponenciação, etc.
Comparação:>, <,> =, <=, == ,! = Etc.
Tarefa: =.
Características: a divisão por / sempre retorna um número fracionário; se você precisar da parte inteira da divisão de dois inteiros, precisará usar a operação div(m, n) ou o infixo equivalente m ÷ n .

Tipos

Tipos numéricos:

• Inteiros ( Int ) - 2 , 3 , -42
• Inteiros não assinados ( UInt ) - 0x12345
• Ponto flutuante ( Float32 , Float64 ) - 1.0 , 3.1415 , -Inf , NaN
• Rational ( Rational ) - 3//3 , 7//2
• Real ( Real ) - todos os itens acima
• Complexo ( Complex ) - 3+4*im , 2//3+2//3*im , 3.0+0.0*im ( im é uma unidade imaginária, apenas um número com uma parte imaginária explicitamente escrita é considerado complexo)
• Number - todos os itens acima

Cordas e caracteres:

• 'a' ( Char )
• "a" é uma string ( String )

Nota: as strings, como agora em muitos idiomas, são imutáveis.
NB: strings (assim como nomes de variáveis) suportam Unicode, incluindo emoji.

Matrizes:

• x = [1, 2, 3] - especificando uma matriz por enumeração direta de elementos
• construtores especiais: zeros(length) para uma matriz de zeros, ones(length) para uma matriz de unidades, rand(length) para uma matriz de números aleatórios, etc.
• suporte a array multidimensional
• suporte para operações de álgebra linear (adição de matrizes, multiplicação escalar, multiplicação de vetores matriciais e muito mais) na biblioteca padrão

Nota: todas as coleções são indexadas a partir de uma.
NB: porque como a linguagem é destinada a tarefas computacionais, as matrizes são um dos tipos mais importantes; você precisará retornar aos princípios de seu trabalho mais de uma vez.

Tuplas (conjunto de elementos ordenado, imutável):

• (2, 5.3, "k") é uma tupla regular
• (a = 3, b = 4) - tupla nomeada

NB: os campos de uma tupla nomeada podem ser acessados ​​pelo nome por um período e pelo índice via []

 julia> x = (a = 5, b = 12) (a = 5, b = 12) julia> x[1] 5 julia> sqrt(xa^2 + x[2]^2) 13.0 

Dicionários:

 julia> x = Dict('a' => 5, 'b' => 12) Dict{Char,Int64} with 2 entries: 'a' => 5 'b' => 12 julia> x['c'] = 13 13 julia> x Dict{Char,Int64} with 3 entries: 'a' => 5 'c' => 13 'b' => 12 

Construções básicas da linguagem de controle

1. As variáveis ​​são criadas automaticamente na atribuição. O tipo é opcional.

 julia> x = 7; x + 2 9 julia> x = 42.0; x * 4 168.0 

2. O bloco de salto condicional começa com a expressão if <condition> e termina com a palavra end . Você também pode ter uma else luz ou outras luzes:

 if x > y println("X is more than Y") elseif x == y println("X and Y are equal") else println("X is less than Y") end 

3. Existem duas construções de loop: while e for . O segundo funciona como em Python, ou seja, Repete a coleção. Um uso comum é iterativo em um intervalo de valores cuja sintaxe é start[:increment]:end . Ao contrário do Python, um intervalo inclui valores iniciais e finais, ou seja, o intervalo vazio não será 1:1 (este é um intervalo de 1), mas 1:0 . O final do corpo do loop é marcado com a palavra end .

 julia> for i in 1:3; print(i, " "); end #   1  3   1 ( ) 1 2 3 julia> for i in 1:2:3; print(i, " "); end #   1  3   2 1 3 

4. As funções são definidas pela palavra-chave function , a definição da função também termina com a palavra end . Argumentos com valores padrão e argumentos nomeados são suportados.

 function square(x) return x * x end function cube(x) x * square(x) #       ; return   end function root(x, degree = 2) #  degree     return x^(1.0/degree) end function greeting(name; times = 42, greet = "hello") #       println(times, " times ", greet, " to ", name) end julia> greeting("John") 42 times hello to John julia> greeting("Mike", greet = "wassup", times = 100500) #           100500 times wassup to Mike 

Em geral, tudo isso é bastante semelhante ao Python, exceto por pequenas diferenças na sintaxe e pelo fato de os blocos de código não serem alocados com espaços, mas ainda com palavras-chave. Em casos simples, os programas Python se traduzem em Julia quase um a um.
Mas há uma diferença significativa no fato de que em Julia você pode especificar explicitamente tipos de variáveis, o que permite compilar programas, obtendo código rápido.
A segunda diferença significativa é que o Python implementa um modelo de POO “clássico” com classes e métodos, enquanto Julia implementa um modelo de despacho múltiplo.

Anotações de tipo e expedição múltipla

Vamos ver o que é uma função interna:

 julia> sqrt sqrt (generic function with 19 methods) 

Como o REPL nos mostra, o sqrt é uma função genérica com 19 métodos. Que tipo de função generalizada e que tipo de métodos?

E isso significa que existem várias funções sqrt que se aplicam a diferentes tipos de argumentos e, portanto, calculam a raiz quadrada usando vários algoritmos. Você pode ver quais opções estão disponíveis, digitando

 julia> methods(sqrt) 

Pode-se observar que a função é definida para diferentes tipos de números, bem como para matrizes.

Diferente do OOP “clássico”, em que a implementação concreta do método é determinada apenas pela classe de chamada (despachada pelo primeiro argumento), em Julia a escolha de uma função é determinada pelos tipos (e número) de todos os seus argumentos.
Ao chamar uma função com argumentos específicos de todos os seus métodos, é selecionado um que descreva com mais precisão o conjunto específico de tipos com os quais a função é chamada e é ela que é usada.

Uma característica distintiva é que é aplicada uma abordagem chamada compilação "antecipadamente" pelos autores da linguagem. I.e. As funções são compiladas para os tipos de dados fornecidos na primeira chamada, após o que as seguintes chamadas são realizadas muito mais rapidamente. A diferença entre a primeira e a subseqüente chamada pode ser muito significativa:

 julia> @time sqrt(8) #  @time -      0.006811 seconds (3.15 k allocations: 168.516 KiB) #   ,        2.8284271247461903 julia> @time sqrt(15) 0.000002 seconds (5 allocations: 176 bytes) # 5   -     @time 3.872983346207417 

No caso ruim, cada chamada de função é uma verificação de tipo dos argumentos recebidos e a procura do método desejado na lista. No entanto, se você der dicas ao compilador, poderá eliminar as verificações, o que levará a um código mais rápido.

Por exemplo, considere calcular a soma

$$

$$\ sum_ {k = 1} ^ N \ sqrt {(- 1) ^ k}$$

 function mysqrt(num) #    -     #   -           if num >= 0 return sqrt(num) else return sqrt(complex(num)) end end function S(n) #    sum = 0 sgn = -1 for k = 1:n sum += mysqrt(sgn) sgn = -sgn end return sum end function S_typed(n::Integer) # ..     ,      #     sum::Complex = 0.0 sgn::Int = -1 for k = 1:n sum += mysqrt(sgn) sgn = -sgn end return sum end 

O benchmark mostra que a função S_typed() não é executada apenas mais rapidamente, mas também não requer alocação de memória para cada chamada, ao contrário de S() . O problema aqui é que o tipo do mysqrt() retornado do mysqrt() não está definido, assim como o tipo do lado direito da expressão

 sum = sum + mysqrt(sgn) 

Como resultado, o compilador não consegue nem descobrir qual será a sum tipo a cada iteração. Portanto, o boxe (vinculação de etiqueta de tipo) é uma variável e a memória é alocada.
Para a função S_typed() , o compilador sabe de antemão que sum é um valor complexo, portanto o código é mais otimizado (em particular, chamar mysqrt() pode ser efetivamente incorporado, retornando sempre o valor de retorno para Complex ).

Mais importante, para S_typed() compilador sabe que o valor de retorno é do tipo Complex , mas para S() tipo de valor de saída não é definido novamente, o que desacelerará todas as funções em que S() será chamado.
Você pode verificar se o compilador pensa nos tipos retornados da expressão usando a macro @code_warntype :

 julia> @code_warntype S(3) Body::Any #     ,      ... julia> @code_warntype S_typed(3) Body::Complex{Float64} #      ... 

Se uma função é chamada em algum lugar do loop para o qual @code_warntype não pode imprimir o tipo de retorno, ou para o qual em algum lugar do corpo mostra o recebimento de um valor do tipo Any , a otimização dessas chamadas provavelmente dará um aumento de desempenho muito tangível.

Tipos de compostos

Um programador pode definir tipos de dados compostos para suas necessidades usando a construção struct :

 julia> struct GenericStruct #   struct    name b::Int c::Char v::Vector end #       #       ,        julia> s = GenericStruct("Name", 1, 'z', [3., 0]) GenericStruct("Name", 1, 'z', [3.0, 0.0]) julia> s.name, sb, sc, sv ("Name", 1, 'z', [3.0, 0.0]) 

As estruturas em Julia são imutáveis, ou seja, ao criar uma instância da estrutura, não é mais possível alterar os valores dos campos (mais precisamente, você não pode alterar o endereço dos campos na memória - elementos de campos mutáveis, como sv no exemplo acima, podem ser alterados). Estruturas mutable struct são criadas pela construção mutable struct , cuja sintaxe é igual à das estruturas regulares.

A herança de estruturas no sentido “clássico” não é suportada, no entanto, existe a possibilidade de “herdar” o comportamento combinando tipos compostos em supertipos ou, como são chamados em Julia, tipos abstratos. Os relacionamentos de tipo são expressos como A<:B (A é um subtipo de B) e A>:B (A é um subtipo de B). Parece algo como isto:

 abstract type NDimPoint end #   -     # ,    -     N  struct PointScalar<:NDimPoint x1::Real end struct Point2D<:NDimPoint x1::Real x2::Real end struct Point3D<:NDimPoint x1::Real x2::Real x3::Real end #     ;   Markdown """ mag(p::NDimPoint) Calculate the magnitude of the radius vector of an N-dimensional point `p` """ function mag(p::NDimPoint) sqrmag = 0.0 # ..   ,       #     T   fieldnames(T) for name in fieldnames(typeof(p)) sqrmag += getfield(p, name)^2 end return sqrt(sqrmag) end """ add(p1::T, p2::T) where T<:NDimPoint Calculate the sum of the radius vectors of two N-dimensional points `p1` and `p2` """ function add(p1::T, p2::T) where T<:NDimPoint #  -  , ..       #     list comprehension sumvector = [Float64(getfield(p1, name) + getfield(p1, name)) for name in fieldnames(T)] #     ,    #  ...      , .. # f([1, 2, 3]...) -   ,  f(1, 2, 3) return T(sumvector...) end 

Estudo de caso: Polinômios

Um sistema de tipos acoplado ao envio múltiplo é conveniente para expressar conceitos matemáticos. Vejamos um exemplo de uma biblioteca simples para trabalhar com polinômios.
Introduzimos dois tipos de polinômios: "canônico", definido através de coeficientes de potência, e "interpolação", definida por um conjunto de pares (x, f (x)). Para simplificar, consideraremos apenas argumentos válidos.

Para armazenar um polinômio em uma notação usual, é adequada uma estrutura com uma matriz ou uma tupla de coeficientes como um campo. Para ser completamente imutável, faça uma carreata. Assim, o código para definir o tipo abstrato, a estrutura do polinômio e calcular o valor do polinômio em um determinado ponto é bastante simples:

 abstract type AbstractPolynomial end """ Polynomial <: AbstractPolynomial Polynomials written in the canonical form """ struct Polynomial<:AbstractPolynomial degree::Int coeff::NTuple{N, Float64} where N # NTuple{N, Type} -    N    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 

Polinômios de interpolação precisam de uma estrutura de representação e método de cálculo diferentes. Em particular, se o conjunto de pontos de interpolação é conhecido antecipadamente, e está planejado calcular o mesmo polinômio em pontos diferentes, a fórmula de interpolação de Newton é conveniente:

$$

$$P (x) = \ sum_ {k = 0} ^ N {c_k n_k (x)},$$

onde n _k ( x ) são polinômios básicos, n ₀ ( x ) e para k > 0

$$

$$n_k (x) = \ prod_ {i = 0} ^ {k-1} {(x-x_i)},$$

onde x _i são os nós de interpolação.

A partir das fórmulas acima, pode-se observar que o armazenamento é convenientemente organizado na forma de um conjunto de nós de interpolação x _i e coeficientes c _i , e o cálculo pode ser feito de maneira semelhante ao esquema de Horner.

 """ InterpPolynomial <: AbstractPolynomial Interpolation polynomials in Newton's form """ struct InterpPolynomial<:AbstractPolynomial degree::Int xval::NTuple{N, Float64} where N coeff::NTuple{N, Float64} where N end """ evpoly(p::Polynomial, z::Real) Evaluate polynomial `p` at `z` using the Horner's rule """ 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 

A função para calcular o valor do polinômio em ambos os casos é chamada de mesma - evpoly() - mas aceita diferentes tipos de argumentos.

Além da função de cálculo, seria bom escrever uma função que cria um polinômio a partir de dados conhecidos.

Existem duas técnicas para isso em Julia: construtores externos e construtores internos. Um construtor externo é simplesmente uma função que retorna um objeto do tipo apropriado. Um construtor interno é uma função que é introduzida na descrição da estrutura e substitui o construtor padrão. É recomendável usar o construtor interno para construir polinômios de interpolação, pois

• é mais conveniente obter um polinômio não através dos nós e coeficientes de interpolação, mas através dos nós e valores da função interpolada
• nós de interpolação devem ser distintos
• o número de nós e coeficientes deve corresponder

Escrever um construtor interno no qual essas regras são garantidas para serem seguidas garante que todas as variáveis ​​criadas do tipo InterpPolynomial , pelo menos, possam ser processadas corretamente pela função evpoly() .

Escrevemos um construtor de polinômios comuns que usa uma matriz unidimensional ou uma tupla de coeficientes como entrada. O construtor do polinômio de interpolação recebe os nós de interpolação e os valores desejados neles e usa o método das diferenças divididas para calcular os coeficientes.

 """ 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}}) #     /     P(x) ≡ 0 coeff = isempty(v) ? (0.0,) : tuple([Float64(x) for x in v]...) #   -   new #  -    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}}} #   ,    ,   f  ,   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] #     (Stoer, Bulirsch, Introduction to Numerical Analysis, . 2.1.3) for i = 2:N for j = 1:(i-1) coeff[i] = (coeff[j] - coeff[i]) / (xsample[j] - xsample[i]) end end new(N-1, tuple([Float64(x) for x in xsample]...), tuple(coeff...)) end end 

Além da geração real de polinômios, seria bom poder executar operações aritméticas com eles.

Como os operadores aritméticos em Julia são funções comuns, às quais uma notação de infixo é adicionada como açúcar sintático (as expressões a + b e +(a, b) são válidas e absolutamente idênticas), sua sobrecarga é feita da mesma maneira que a escrita métodos adicionais para suas funções.

O único ponto sutil é que o código do usuário é iniciado a partir do módulo Main (namespace) e as funções da biblioteca padrão estão no módulo Base , portanto, ao sobrecarregar, você deve importar o módulo Base ou escrever o nome completo da função.

Então, adicionamos a adição de um polinômio com um número:

 # -   Base.+  , #    Base.:+,   " :+   Base" function Base.:+(p::Polynomial, x::Real) Polynomial(tuple(p.coeff[1] + x, p.coeff[2:end]...)) end function Base.:+(p::InterpPolynomial, x::Real) # ..           - #          . #       - #        fval::Vector{Float64} = [evpoly(p, xval) + x for xval in p.xval] InterpPolynomial(p.xval, fval) end #       function Base.:+(x::Real, p::AbstractPolynomial) return p + x end 

Para adicionar dois polinômios comuns, basta adicionar os coeficientes e, ao adicionar o polinômio de interpolação ao outro, você pode encontrar os valores da soma em vários pontos e criar uma nova interpolação a partir deles.

 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 

Da mesma forma, você pode adicionar outras operações aritméticas em polinômios, resultando em sua representação no código em uma notação matemática natural.

Por enquanto é tudo. Vou tentar escrever mais sobre a implementação de outros métodos numéricos.

Na preparação, foram utilizados os seguintes materiais:

1. Documentação no idioma Julia: docs.julialang.org
2. Plataforma de discussão de idiomas Julia: discourse.julialang.org
3. J. Stoer, W. Bulirsch. Introdução à Análise Numérica
4. Julia Hub: habr.com/en/hub/julia
5. Pense em Julia: benlauwens.imtqy.com/ThinkJulia.jl/latest/book.html