🤴🏿 👵🏾 🎚️ Julia: types, multiméthodes et arithmétique sur polynômes 🌧️ 🔝 🚤

Cette publication se concentrera sur le principal, à mon avis, trait distinctif de la langue Julia - la présentation des fonctions sous la forme de méthodes avec répartition multiple. Cela vous permet d'augmenter les performances des calculs sans réduire la lisibilité du code et sans gâcher l'abstraction, d'une part, et vous permet de travailler avec des concepts mathématiques dans une notation plus familière, d'autre part. A titre d'exemple, la question de l'uniforme (du point de vue des opérations linéaires) travaille avec des polynômes dans la représentation de la liste des coefficients et avec des polynômes d'interpolation.

Syntaxe de base

Une brève introduction pour ceux qui ne sont pas au courant. Julia est un langage de script, il a REPL (boucle lecture-évaluation-impression, c'est-à-dire un shell interactif). À première vue, il ressemble assez, par exemple, à Python ou MATLAB.

Opérations arithmétiques

L'arithmétique est à peu près la même que partout: +, -, *, /, ^ pour l'exponentiation, etc.
Comparaison:>, <,> =, <=, == ,! = Etc.
Affectation: =.
Caractéristiques: la division par / renvoie toujours un nombre fractionnaire; si vous avez besoin de la partie entière de la division de deux entiers, vous devez utiliser l'opération div(m, n) ou l'équivalent infixe m ÷ n .

Les types

Types numériques:

Entiers ( Int ) - 2 , 3 , -42
UInt non signés ( UInt ) - 0x12345
Virgule flottante ( Float32 , Float64 ) - 1.0 , 3.1415 , -Inf , NaN
Rationnel ( Rational ) - 3//3 , 7//2
Réel ( Real ) - tout ce qui précède
Complexe ( Complex ) - 3+4*im , 2//3+2//3*im , 3.0+0.0*im ( im est une unité imaginaire, seul un nombre avec une partie imaginaire explicitement écrite est considéré comme complexe)
Number - toutes ces réponses

Cordes et caractères:

'a' - caractère ( Char )
"a" est une chaîne ( String )

NB: les chaînes, comme maintenant dans de nombreuses langues, sont immuables.
NB: les chaînes (ainsi que les noms de variables) supportent Unicode, y compris les emoji.

Tableaux:

x = [1, 2, 3] - spécification d'un tableau par énumération directe des éléments
constructeurs spéciaux: zeros(length) pour un tableau de zéros, ones(length) pour un tableau de uns, rand(length) pour un tableau de nombres aléatoires, etc.
prise en charge des baies multidimensionnelles
prise en charge des opérations d'algèbre linéaire (ajout de tableaux, multiplication scalaire, multiplication vectorielle matricielle, etc.) dans la bibliothèque standard

NB: toutes les collections sont indexées à partir d'une seule.
NB: car le langage est destiné aux tâches de calcul, les tableaux sont l'un des types les plus importants, vous devrez revenir plus d'une fois aux principes de leur travail.

Tuples (ensemble ordonné d'éléments, immuable):

(2, 5.3, "k") est un tuple régulier
(a = 3, b = 4) - tuple nommé

NB: les champs d'un tuple nommé sont accessibles à la fois par nom sur une période et par index via []

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

Dictionnaires:

 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

Constructions de base du langage de contrôle

1. Les variables sont créées automatiquement lors de l'affectation. Le type est facultatif.

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

2. Le bloc de saut conditionnel commence par l'expression if <condition> et se termine par le mot end . Vous pouvez également avoir une else lumière ou d'autres lumières:

 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. Il existe deux constructions de boucle: while et for . Le second fonctionne comme en Python, c'est-à-dire Itère sur la collection. Une utilisation courante est l'itération sur une plage de valeurs dont la syntaxe est start[:increment]:end . Contrairement à Python, une plage comprend à la fois des valeurs de début et de fin, c'est-à-dire la plage vide ne sera pas 1:1 (c'est une plage de 1), mais 1:0 . La fin du corps de la boucle est marquée par le mot 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. Les fonctions sont définies par le mot-clé function , la définition de la fonction se termine également par le mot end . Les arguments avec des valeurs par défaut et des arguments nommés sont pris en charge.

 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

En général, tout cela est assez similaire à Python, à l'exception de différences mineures dans la syntaxe et du fait que les blocs de code ne sont pas alloués avec des espaces, mais toujours avec des mots clés. Dans les cas simples, les programmes Python se traduisent même en Julia presque un à un.
Mais il y a une différence significative dans le fait que dans Julia, vous pouvez spécifier explicitement les types de variables, ce qui vous permet de compiler des programmes, d'obtenir du code rapide.
La deuxième différence significative est que Python implémente un modèle OOP «classique» avec des classes et des méthodes, tandis que Julia implémente un modèle multi-dispatch.

Annotations de type et répartition multiple

Voyons ce qu'est une fonction intégrée:

 julia> sqrt sqrt (generic function with 19 methods)

Comme REPL nous le montre, sqrt est une fonction générique avec 19 méthodes. Quel genre de fonction généralisée et quel genre de méthodes?

Et cela signifie qu'il existe plusieurs fonctions sqrt qui s'appliquent à différents types d'arguments et, en conséquence, calculent la racine carrée à l'aide de divers algorithmes. Vous pouvez voir les options disponibles en tapant

 julia> methods(sqrt)

On peut voir que la fonction est définie pour différents types de nombres, ainsi que pour les matrices.

Contrairement à la POO «classique», où l'implémentation concrète de la méthode n'est déterminée que par la classe appelante (répartition par le premier argument), dans Julia le choix d'une fonction est déterminé par les types (et le nombre) de tous ses arguments.
Lors de l'appel d'une fonction avec des arguments spécifiques de toutes ses méthodes, une est sélectionnée qui décrit le plus précisément l'ensemble spécifique de types avec lesquels la fonction est appelée, et c'est elle qui est utilisée.

Un trait distinctif est qu'une approche appelée compilation «juste à l'avance» par les auteurs de la langue est appliquée. C'est-à-dire les fonctions sont compilées pour les types de données donnés lors du premier appel, après quoi les appels suivants sont effectués beaucoup plus rapidement. La différence entre le premier appel et les appels suivants peut être très importante:

 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

Dans le mauvais cas, chaque appel de fonction est une vérification du type des arguments reçus et une recherche de la méthode souhaitée dans la liste. Cependant, si vous donnez des conseils au compilateur, vous pouvez éliminer les vérifications, ce qui entraînera un code plus rapide.

Par exemple, envisagez de calculer la somme

s u m_{k = 1}^{N} s q r t (- 1)^{k}

$\ 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

Le benchmark montre que la fonction S_typed() fonctionne non seulement plus rapidement, mais ne nécessite pas non plus d'allocation de mémoire pour chaque appel, contrairement à S() . Le problème ici est que le type de la mysqrt() retournée par mysqrt() n'est pas défini, tout comme le type du côté droit de l'expression

 sum = sum + mysqrt(sgn)

Par conséquent, le compilateur ne peut même pas déterminer quel type de sum sera à chaque itération. Ainsi, la boxe (accrochage d'étiquette de type) est une variable et la mémoire est allouée.
Pour la fonction S_typed() , le compilateur sait à l'avance que sum est une valeur complexe, donc le code est plus optimisé (en particulier, appeler mysqrt() peut être effectivement en ligne, renvoyant toujours la valeur de retour à Complex ).

Plus important encore, pour S_typed() compilateur sait que la valeur de retour est de type Complex , mais pour S() type de la valeur de sortie n'est pas défini à nouveau, ce qui ralentira toutes les fonctions où S() sera appelé.
Vous pouvez vérifier que le compilateur pense aux types renvoyés par l'expression à l'aide de la macro @code_warntype :

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

Si une fonction est appelée quelque part dans la boucle pour laquelle @code_warntype ne peut pas imprimer le type de retour, ou pour laquelle elle quelque part dans le corps montre la réception d'une valeur de type Any , alors l'optimisation de ces appels donnera très probablement une amélioration des performances très tangible.

Types de composés

Un programmeur peut définir des types de données composites pour ses besoins en utilisant la construction 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])

Les structures dans Julia sont immuables, c'est-à-dire qu'en créant une instance de la structure, il n'est plus possible de modifier les valeurs des champs (plus précisément, vous ne pouvez pas modifier l'adresse des champs en mémoire - des éléments de champs modifiables, tels que sv dans l'exemple ci-dessus, peuvent être modifiés). Les structures mutable struct sont créées par la construction de mutable struct , dont la syntaxe est la même que pour les structures régulières.

L'héritage des structures au sens «classique» n'est pas pris en charge, mais il existe la possibilité «d'hériter» du comportement en combinant des types composites en supertypes ou, comme on les appelle dans Julia, des types abstraits. Les relations de type sont exprimées comme A<:B (A est un sous-type de B) et A>:B (A est un sous-type de B). Cela ressemble à ceci:

 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

Étude de cas: polynômes

Un système de type couplé à une répartition multiple est pratique pour exprimer des concepts mathématiques. Voyons un exemple de bibliothèque simple pour travailler avec des polynômes.
Nous introduisons deux types de polynômes: «canoniques», définis par des coefficients aux puissances, et «interpolation», définis par un ensemble de paires (x, f (x)). Par souci de simplicité, nous ne considérerons que des arguments valides.

Pour stocker un polynôme dans une notation habituelle, une structure ayant un tableau ou un tuple de coefficients comme champ convient. Pour être complètement immuable, qu'il y ait un cortège. Ainsi, le code pour définir le type abstrait, la structure du polynôme et calculer la valeur du polynôme en un point donné est assez simple:

 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

Les polynômes d'interpolation nécessitent une structure de représentation et une méthode de calcul différentes. En particulier, si l'ensemble des points d'interpolation est connu à l'avance et qu'il est prévu de calculer le même polynôme à différents points, la formule d'interpolation de Newton est pratique:

P (x) = s u m_{k = 0}^{N} c_{k} n_{k} (x),

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

où n _k ( x ) sont des polynômes de base, n ₀ ( x ) et pour k > 0

n_{k} (x) = p r o d_{i = 0}^{k - 1} (x - x_{i}),

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

où x _i sont les nœuds d'interpolation.

À partir des formules ci-dessus, on peut voir que le stockage est commodément organisé sous la forme d'un ensemble de nœuds d'interpolation x _i et de coefficients c _i , et le calcul peut être effectué d'une manière similaire au schéma 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

La fonction de calcul de la valeur du polynôme dans les deux cas est appelée la même - evpoly() - mais accepte différents types d'arguments.

En plus de la fonction de calcul, il serait bien d'écrire une fonction qui crée un polynôme à partir de données connues.

Il existe deux techniques pour cela dans Julia: les constructeurs externes et les constructeurs internes. Un constructeur externe est simplement une fonction qui renvoie un objet du type approprié. Un constructeur interne est une fonction qui est introduite dans la description de la structure et remplace le constructeur standard. Il est conseillé d'utiliser le constructeur interne pour construire des polynômes d'interpolation, car

il est plus pratique d'obtenir un polynôme non pas à travers les nœuds d'interpolation et les coefficients, mais à travers les nœuds et les valeurs de la fonction interpolée
les nœuds d'interpolation doivent être distincts
le nombre de nœuds et de coefficients doit correspondre

L'écriture d'un constructeur interne dans lequel ces règles sont garanties d'être respectées garantit que toutes les variables créées du type InterpPolynomial , au moins, peuvent être correctement traitées par la fonction evpoly() .

Nous écrivons un constructeur de polynômes ordinaires qui prend comme entrée un tableau unidimensionnel ou un tuple de coefficients. Le constructeur du polynôme d'interpolation reçoit les nœuds d'interpolation et les valeurs souhaitées en eux et utilise la méthode des différences divisées pour calculer les coefficients.

 """ 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

Outre la génération réelle de polynômes, il serait bien de pouvoir effectuer des opérations arithmétiques avec eux.

Étant donné que les opérateurs arithmétiques de Julia sont des fonctions ordinaires, auxquelles une notation infixe est ajoutée en tant que sucre syntaxique (les expressions a + b et +(a, b) sont à la fois valides et absolument identiques), leur surcharge se fait de la même manière que l'écriture méthodes supplémentaires à leurs fonctions.

Le seul point subtil est que le code utilisateur est lancé à partir du module Main (espace de noms) et que les fonctions de la bibliothèque standard sont dans le module de Base , donc lors d'une surcharge, vous devez soit importer le module de Base soit écrire le nom complet de la fonction.

Donc, nous ajoutons l'ajout d'un polynôme avec un nombre:

 # -   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

Pour ajouter deux polynômes ordinaires, il suffit d'ajouter les coefficients, et lorsque vous ajoutez le polynôme d'interpolation à l'autre, vous pouvez trouver les valeurs de somme à plusieurs points et construire une nouvelle interpolation à partir d'eux.

 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

De la même manière, vous pouvez ajouter d'autres opérations arithmétiques sur des polynômes, résultant en leur représentation dans le code dans une notation mathématique naturelle.

C'est tout pour l'instant. J'essaierai d'écrire davantage sur la mise en œuvre d'autres méthodes numériques.

En préparation, les matériaux suivants ont été utilisés:

Documentation en langue Julia: docs.julialang.org
Plateforme de discussion linguistique Julia: discourse.julialang.org
J. Stoer, W. Bulirsch. Introduction à l'analyse numérique
Julia Hub: habr.com/en/hub/julia
Pensez Julia: benlauwens.imtqy.com/ThinkJulia.jl/latest/book.html

Julia: types, multiméthodes et arithmétique sur polynômes