--- engine: julia --- # A Case Study: The Parametric Data Type PComplex ```{julia} #| error: false #| echo: false #| output: false using InteractiveUtils import QuartoNotebookWorker Base.stdout = QuartoNotebookWorker.with_context(stdout) myactive_module() = Main.Notebook Base.active_module() = myactive_module() # https://github.com/JuliaLang/julia/blob/master/base/show.jl#L516-L520 # https://github.com/JuliaLang/julia/blob/master/base/show.jl#L3073-L3077 ``` We want to introduce a new numeric type **complex numbers in polar representation $z=r e^{i\phi}=(r,\phi)$**. - The type should integrate into the type hierarchy as a subtype of 'Number'. - $r$ and $\phi$ should be floating point numbers. (Unlike complex numbers in 'Cartesian' coordinates, restricting to integer values of r or $\phi$ makes little mathematical sense.) ## The Definition of `PComplex` A first attempt could look like this: ```{julia} struct PComplex1{T <: AbstractFloat} <: Number r :: T ϕ :: T end z1 = PComplex1(-32.0, 33.0) z2 = PComplex1{Float32}(12, 13) @show z1 z2; ``` :::{.callout-warning collapse="true" .titlenormal} ## It is not possible to redefine a `struct` once it has been defined in a Julia session. Therefore, I use different names. Another possibility is, for example, the use of [`ProtoStructs.jl`](https://juliahub.com/ui/Packages/General/ProtoStructs). ::: Julia automatically provides *default constructors*: - the constructor `PComplex1`, where the type `T` is inferred from the passed arguments, and - constructors `PComplex{Float64},...` with explicit type specification. Here, the arguments are attempted to be converted to the requested type. ------ We now want the constructor to do even more. In the polar representation, we want $0\le r$ and $0\le \phi<2\pi$ to hold. If the passed arguments do not satisfy this, they should be recalculated accordingly. To this end, we define an _inner constructor_ that replaces the _default constructor_. - An _inner constructor_ is a function within the `struct` definition. - In an _inner constructor_, one can use the special function `new`, which acts like the _default constructor_. ```{julia} struct PComplex{T <: AbstractFloat} <: Number r :: T ϕ :: T function PComplex{T}(r::T, ϕ::T) where T<:AbstractFloat if r<0 # flip the sign of r and correct phi r = -r ϕ += π end if r==0 ϕ=0 end # normalize r=0 case to phi=0 ϕ = mod(ϕ, 2π) # map phi into interval [0,2pi) new(r, ϕ) # new() is special function, end # available only inside inner constructors end ``` ```{julia} #| echo: false #| output: false #= in the whole quarto-runs we want to use the default show here =# zz = @which Base.show(stdout, PComplex{Float64}(2.,3.)) if zz.module != Base Base.delete_method(zz) end ``` ```{julia} z1 = PComplex{Float64}(-3.3, 7π+1) ``` For the explicit specification of an *inner constructor*, we have to pay a price: The *default constructors* provided by Julia are missing. The constructor without explicit type specification in curly braces, which takes over the type of the arguments, is also desired: ```{julia} PComplex(r::T, ϕ::T) where {T<:AbstractFloat} = PComplex{T}(r,ϕ) z2 = PComplex(2.0, 0.3) ``` ## A New Notation Julia uses `//` as an infix constructor for the type `Rational`. We want something equally nice. In electronics/electrical engineering, [AC quantities are described by complex numbers.](https://en.wikipedia.org/wiki/Phasor_analysis) A representation of complex numbers by "magnitude" and "phase" is common and is often represented in so-called [phasor form](https://en.wikipedia.org/wiki/Phasor): $$ z= r\enclose{phasorangle}{\phi} = 3.4\;\enclose{phasorangle}{45^\circ} $$ where the angle is usually noted in degrees. :::{.callout-note .titlenormal collapse="true"} ## Possible Infix Operators in Julia In Julia, a large number of Unicode characters are reserved for use as operators. The definitive list is in the [parser source code.](https://github.com/JuliaLang/julia/blob/eaa2c58aeb12f27c1d8c116ab111773a4fc4495f/src/julia-parser.scm#L13-L31) Details will be discussed in a later chapter. ::: The angle bracket symbol `∠` is not available as an operator symbol. We use `⋖` as an alternative. This can be entered in Julia as `\lessdot`. ```{julia} ⋖(r::Real, ϕ::Real) = PComplex(r, π*ϕ/180) z3 = 2. ⋖ 90. ``` (The type annotation -- `Real` instead of `AbstractFloat` -- is a preview of further constructors to come. For now, the operator `⋖` only works with `Float`s.) Of course, we also want the output to look nice. Details can be found in the [documentation](https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing). ```{julia} using Printf function Base.show(io::IO, z::PComplex) # we print the phase in degrees, rounded to tenths of a degree, p = z.ϕ * 180/π sp = @sprintf "%.1f" p print(io, z.r, "⋖", sp, '°') end @show z3; ``` ## Methods for `PComplex` For our type to be a proper member of the family of types derived from `Number`, we need a whole lot more. Arithmetic, comparison operators, conversions, etc. must be defined. We limit ourselves to multiplication and square roots. :::{.callout-note collapse="true"} ## Modules - To add to the `methods` of existing functions and operations, one must address them with their 'full name'. - All objects belong to a namespace or `module`. - Most basic functions belong to the module `Base`, which is always loaded without explicit `using ...` by default. - As long as one does not define own modules, own definitions are in the module `Main`. - The macro `@which`, applied to a name, shows in which module the name is defined. ```{julia} f(x) = 3x^3 @which f ``` ```{julia} wp = @which + ws = @which(sqrt) println("Module for addition: $wp, Module for sqrt: $ws") ``` ::: ```{julia} qwurzel(z::PComplex) = PComplex(sqrt(z.r), z.ϕ / 2) ``` ```{julia} #| echo: false #| output: false #= to make length(methods(sqrt)) work =# if hasmethod(sqrt, (PComplex,)) zz = @which Base.sqrt(PComplex{Float64}(1.,1.)) Base.delete_method(zz) end ``` The function `sqrt()` already has some methods: ```{julia} length(methods(sqrt)) ``` Now it will have one more method: ```{julia} Base.sqrt(z::PComplex) = qwurzel(z) length(methods(sqrt)) ``` ```{julia} sqrt(z2) ``` and now for multiplication: ```{julia} Base.:*(x::PComplex, y::PComplex) = PComplex(x.r * y.r, x.ϕ + y.ϕ) @show z1 * z2; ``` (Since the operator symbol is not a normal name, the colon must be with `Base.` in the composition.) We can, however, not yet multiply with other numeric types. One could now define a large number of corresponding methods. Julia provides one more mechanism for *numeric types* that simplifies this somewhat. ## Type Promotion and Conversion In Julia, one can naturally use the most diverse numeric types side by side. ```{julia} 1//3 + 5 + 5.2 + 0xff ``` If one looks at the numerous methods defined, for example, for `+` and `*`, one finds among them a kind of 'catch-all definition' ```julia +(x::Number, y::Number) = +(promote(x,y)...) *(x::Number, y::Number) = *(promote(x,y)...) ``` (The 3 dots are the splat operator, which decomposes the tuple returned by promote() back into its components.) Since the method with the types `(Number, Number)` is very general, it is only used when more specific methods do not apply. What happens here? ### The Function `promote(x,y,...)` This function attempts to convert all arguments to a common type that can represent all values (as precisely as possible). ```{julia} promote(12, 34.555, 77/99, 0xff) ``` ```{julia} z = promote(BigInt(33), 27) @show z typeof(z); ``` The function `promote()` uses two helpers, the functions `promote_type(T1, T2)` and `convert(T, x)` As usual in Julia, one can extend this mechanism with [custom *promotion rules* and `convert(T,x)` methods.](https://docs.julialang.org/en/v1/manual/conversion-and-promotion/) ### The Function `promote_type(T1, T2,...)` It determines to which type conversion should take place. Arguments are types, not values. ```{julia} @show promote_type(Rational{Int64}, ComplexF64, Float32); ``` ### The Function `convert(T,x)` The methods of `convert(T, x)` convert `x` into an object of type `T`. Such a conversion should be lossless. ```{julia} z = convert(Float64, 3) ``` ```{julia} z = convert(Int64, 23.00) ``` ```{julia} z = convert(Int64, 2.3) ``` The special role of `convert()` lies in the fact that it is used *implicitly* and automatically at various points: > [The following language constructs call convert](https://docs.julialang.org/en/v1/manual/conversion-and-promotion/#When-is-convert-called?): > - Assigning to an array converts to the array's element type. - Assigning to a field of an object converts to the declared type of the field. - Constructing an object with new converts to the object's declared field types. - Assigning to a variable with a declared type (e.g. local x::T) converts to that type. - A function with a declared return type converts its return value to that type. -- and of course in `promote()` For self-defined data types, one can extend convert() with further methods. For data types within the Number hierarchy, there is again a 'catch-all definition' ```julia convert(::Type{T}, x::Number) where {T<:Number} = T(x) ``` So: If for a type `T` from the hierarchy `T<:Number` there exists a constructor `T(x)` with a numeric argument `x`, then this constructor `T(x)` is automatically used for conversions. (Of course, more specific methods for `convert()` can also be defined, which then have priority.) ### Further Constructors for `PComplex` ```{julia} ## (a) r, ϕ arbitrary reals, e.g. Integers, Rationals PComplex{T}(r::T1, ϕ::T2) where {T<:AbstractFloat, T1<:Real, T2<: Real} = PComplex{T}(convert(T, r), convert(T, ϕ)) PComplex(r::T1, ϕ::T2) where {T1<:Real, T2<: Real} = PComplex{promote_type(Float64, T1, T2)}(r, ϕ) ## (b) For conversion from reals: constructor with ## only one argument r PComplex{T}(r::S) where {T<:AbstractFloat, S<:Real} = PComplex{T}(convert(T, r), convert(T, 0)) PComplex(r::S) where {S<:Real} = PComplex{promote_type(Float64, S)}(r, 0.0) ## (c) Conversion Complex -> PComplex PComplex{T}(z::Complex{S}) where {T<:AbstractFloat, S<:Real} = PComplex{T}(abs(z), angle(z)) PComplex(z::Complex{S}) where {S<:Real} = PComplex{promote_type(Float64, S)}(abs(z), angle(z)) ``` A test of the new constructors: ```{julia} 3//5 ⋖ 45, PComplex(Complex(1,1)), PComplex(-13) ``` We now still need *promotion rules* that determine which type should result from `promote(x::T1, y::T2)`. This internally extends `promote_type()` with the necessary further methods. ### *Promotion rules* for `PComplex` ```{julia} Base.promote_rule(::Type{PComplex{T}}, ::Type{S}) where {T<:AbstractFloat,S<:Real} = PComplex{promote_type(T,S)} Base.promote_rule(::Type{PComplex{T}}, ::Type{Complex{S}}) where {T<:AbstractFloat,S<:Real} = PComplex{promote_type(T,S)} ``` 1. Rule: : If a `PComplex{T}` and an `S<:Real` meet, then both should be converted to `PComplex{U}`, where `U` is the type to which `S` and `T` can both be converted (_promoted_). 2. Rule : If a `PComplex{T}` and a `Complex{S}` meet, then both should be converted to `PComplex{U}`, where `U` is the type to which `S` and `T` can be converted. Now multiplication with arbitrary numeric types works. ```{julia} z3, 3z3 ``` ```{julia} (3.0+2im) * (12⋖30.3), 12sqrt(z2) ``` :::{.callout-caution icon="false" collapse="true" .titlenormal} ## Summary: our type `PComplex` ```julia struct PComplex{T <: AbstractFloat} <: Number r :: T ϕ :: T function PComplex{T}(r::T, ϕ::T) where T<:AbstractFloat if r<0 # flip the sign of r and correct phi r = -r ϕ += π end if r==0 ϕ=0 end # normalize r=0 case to phi=0 ϕ = mod(ϕ, 2π) # map phi into interval [0,2pi) new(r, ϕ) # new() is special function, end # available only inside inner constructors end # additional constructors PComplex(r::T, ϕ::T) where {T<:AbstractFloat} = PComplex{T}(r,ϕ) PComplex{T}(r::T1, ϕ::T2) where {T<:AbstractFloat, T1<:Real, T2<: Real} = PComplex{T}(convert(T, r), convert(T, ϕ)) PComplex(r::T1, ϕ::T2) where {T1<:Real, T2<: Real} = PComplex{promote_type(Float64, T1, T2)}(r, ϕ) PComplex{T}(r::S) where {T<:AbstractFloat, S<:Real} = PComplex{T}(convert(T, r), convert(T, 0)) PComplex(r::S) where {S<:Real} = PComplex{promote_type(Float64, S)}(r, 0.0) PComplex{T}(z::Complex{S}) where {T<:AbstractFloat, S<:Real} = PComplex{T}(abs(z), angle(z)) PComplex(z::Complex{S}) where {S<:Real} = PComplex{promote_type(Float64, S)}(abs(z), angle(z)) # nice input ⋖(r::Real, ϕ::Real) = PComplex(r, π*ϕ/180) # nice output using Printf function Base.show(io::IO, z::PComplex) # we print the phase in degrees, rounded to tenths of a degree, p = z.ϕ * 180/π sp = @sprintf "%.1f" p print(io, z.r, "⋖", sp, '°') end # arithmetic Base.sqrt(z::PComplex) = PComplex(sqrt(z.r), z.ϕ / 2) Base.:*(x::PComplex, y::PComplex) = PComplex(x.r * y.r, x.ϕ + y.ϕ) # promotion rules Base.promote_rule(::Type{PComplex{T}}, ::Type{S}) where {T<:AbstractFloat,S<:Real} = PComplex{promote_type(T,S)} Base.promote_rule(::Type{PComplex{T}}, ::Type{Complex{S}}) where {T<:AbstractFloat,S<:Real} = PComplex{promote_type(T,S)} ``` ::: :::{.content-hidden unless-format="xxx"} Now something like `PComplex(1, 0)` does not work yet. We also want to allow other real types for `r` and `ϕ`. For simplicity, we convert everything to `Float64` here. We proceed analogously if only one real or complex argument is used. ```julia PComplex(r::Real, ϕ::Real) = PComplex(Float64(r), Float64(ϕ)) PComplex(r::Real) = PComplex(Float64(r), 0.0) PComplex(z::Complex) = PComplex(abs(z), angle(z)) z3 = PComplex(-2); z4 = PComplex(3im) @show z3 z4; ``` :::