⬆️ ⬇️

Julia: functions and structure-as-functions

Despite the fact that the Julia language by design lacks “classical” object-oriented programming with classes and methods, the language provides abstraction tools, in which the type system and elements of functional programming play a key role. Consider the second paragraph in more detail.



The concept of functions in Julia is probably the most similar to the languages ​​of the Lisp family (to be more precise, Lisp-1 branches), and functions can be viewed at three levels: as subroutines, as abstractions to a certain sequence of actions, and as data representing this abstraction .



Level 1. Functions as subroutines



Allocation of subroutines and assigning them their own names has been going on since prehistoric times, when Fortran was considered a high-level language, and C was not yet.



In this understanding, the means of Julia are standard. "Feature" can be called the fact that syntactically there is no division into procedures and functions. Regardless of whether a subroutine is called to get some value or only to perform some actions on data, it is called a function.

')

The function definition begins with the function keyword, followed by the list of arguments in brackets, a sequence of commands, and the word end ends with the definition:



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


The syntax is distinguished by the behavior inherited from Lisp: for a “normal” return of a value from a function, the word return not necessary: ​​the value of the last expression calculated before end returned. In the example above, the value of the variable variable will be returned. Thus, return can be used as a marker for the special behavior of a function:



 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 


For functions with a short definition, there is a shortened syntax, similar to a mathematical notation. Thus, the calculation of the length of the hypotenuse by the length of the legs can be determined as follows:



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


"Safe" division using the ternary operator can be written as:



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


As you can see, it is not necessary to specify types for function arguments. Taking into account how Julia's JIT compiler works, “duck typing” will not even always lead to a drop in performance.



As I tried to demonstrate in the previous article , the Julia compiler can display the type of the returned result by the types of input arguments. Therefore, for example, the safe_division function for quick work requires minimal modification:



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


Now, if the types of both arguments are known at the compilation stage, the type of the returned result is also unambiguously derived, since the zero(x) function returns a zero value of the same type as its argument (and the division by zero, according to IEEE 754 , has a quite representable value in the format of floating-point numbers).



Functions can have a fixed number of positional arguments, positional arguments with default values, named arguments, and a variable number of arguments. Syntax:



 #    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 


Level 2. Functions as data



The name of the function can be used not only in direct calls, but also as an identifier with which the procedure for obtaining a value is associated. For example:



 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 


The “classic” functions that take a functional argument are map , reduce and filter .



map(f, x...) applies the function f to the values ​​of all elements from x (or tuples from the i-th elements) and returns the results as a new collection:



 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) "reduces" the collection to a single value, "expanding" the chain 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 


Since it is not actually determined in what order the array will pass during the reduction, and also whether fn(accum, x) or fn(x, accum) will be fn(x, accum) , the reduction will give a predictable result only with commutative or associative operators, such as addition or multiplication.



filter(predicate, x) returns an array of x elements that satisfy the predicate 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 


Using higher-order functions for array operations instead of writing a loop has several advantages:



  1. the code gets shorter
  2. map() or reduce() show the semantics of the operation being performed, then the semantics of what is happening in the cycle must also be understood
  3. map() allows the compiler to understand that operations on array elements are data-independent, which allows for additional optimizations


Level 3. Functions as abstractions



Often in map() or filter() you need to use a function that has not been assigned its own name. Julia in this case allows you to express an abstraction of performing operations on an argument without entering your own name for this sequence. Such an abstraction is called an anonymous function , or a lambda function (since in the mathematical tradition such functions are denoted by the letter lambda). The syntax for this view is:



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


Both named and anonymous functions can be assigned to variables and returned as values:



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


Variable scope and lexical closures



Normally, functions are attempted to be written in such a way that all the data necessary for the calculation are obtained through formal arguments, i.e. Any variable names that are found in the body are either the names of formal arguments or the names of variables entered inside the function body.



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


About the normal() function, it can be said that in its body all variable names are related , i.e if we everywhere (including the argument list) replace “x” with “m” (or any other identifier), “y” with “n”, and “z” with “sum_of_m_and_n” - the meaning of the expression will not change. In the function strange() name z is unbound , i.e. a) the meaning may change if this name is replaced by another and b) the correctness of the function depends on whether a variable with the name “z” was defined at the time of the function call.



Generally speaking, the normal() function is also not so clean:



  1. What happens if a variable with the name z is defined outside the function?
  2. The + and * characters are also unbound identifiers.


With paragraph 2, nothing can be done except to agree - it is logical that the definitions of all functions used in the system should exist, and we hope that their real meaning meets our expectations.



Point 1 is less obvious than it seems. The point is that the answer depends on where the function is defined. If it is defined globally, then z inside normal() will be a local variable, i.e. even if there is a global variable z its value will not be overwritten. If the definition of a function is inside a block of code, then if there is an earlier definition of z in this block, the value of the external variable will be changed.



If the function body contains the name of an external variable, then this name is associated with the value that existed in the environment where the function was created. If the function itself is exported from this environment (for example, if it is returned from another function as a value), then it “captures” the variable from the internal environment, to which the new environment no longer has access. This is called lexical closure.



Closures are mainly useful in two situations: when you need to create a function according to specified parameters and when you need a function that has some internal state.



Consider the situation with the function that encapsulates the internal state:



 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: all the same polynomials



In the last article the presentation of polynomials as structures was considered. In particular, one of the storage structures is a list of coefficients, starting with the youngest. To calculate the polynomial p at the point x proposed to call the function evpoly(p, x) , which calculates the polynomial according to the Horner scheme.



Full definition code
 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 




The representation of a polynomial in the form of a structure does not quite correspond to its intuitive understanding as a mathematical function. But with the help of returning a functional value, you can set polynomials and directly as functions. So, it was:



 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 


Let us redo this definition into a function that accepts an array / tuple of coefficients and returns the function itself, which computes a polynomial:

 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 


Similarly, you can write a function for the interpolation polynomial.



An important question: was there anything in the previous definition that was lost in the new one? Unfortunately, yes - the task of a polynomial as a structure gave hints for the compiler, and for us the possibility of overloading arithmetic operators for this structure. For the functions of such a powerful type system, Julia, alas, does not provide.



Fortunately, in this case, we can take the best of both worlds, since Julia allows you to create so-called callable structs. Those. You can set a polynomial as a structure, but be able to call it as a function! To the definitions of structures from the previous article you need only add:



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


With the help of functional arguments, you can also add an external interpolation polynomial constructor for some function, built on a set of points:



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


Check the definition
 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 




Conclusion



The features borrowed from functional programming in Julia give a greater expressiveness of the language compared to the purely imperative style. Representation of structures in the form of functions is a way of more convenient and natural recording of mathematical concepts.

Source: https://habr.com/ru/post/451172/



All Articles