Overhaul standardization, complete the set of index conversions

This makes `standardize` more consistent; it's no longer mutating and
no longer requires a fully specified vector for Noll.

The new index types allow coefficient vector conversions from Fringe to
Noll and vice versa as well as from standard indexed order.

The documentation is getting pretty unruly, I should bring Documenter
into the mix and throw up a standard GitHub Page; it should be easy
enough to generate from the docstrings, but I kinda like everything
there in the README.

Also, I might need to write more tests under "standard".

README.md

@@ -175,12 +175,14 @@ This package uses the ANSI Z80.28-2004 standard sequential ordering scheme where
 
 * `noll_to_j(noll::Int)`: converts Noll indices to ANSI standard indices;
 * `j_to_noll(j::Int)`: converts ANSI standard indices to Noll indices;
-* `standardize!(noll::Vector)`: re-orders a Noll specified Zernike expansion coefficient vector according to the ANSI standard; this requires a full ordered vector up to `n_max`;
+* `standardize(noll::Noll)`: re-orders a Noll specified Zernike expansion coefficient vector according to the ANSI standard;
 * `fringe_to_j(fringe::Int)`: converts Fringe indices to ANSI standard indices; only indices 1:37 are valid;
 * `j_to_fringe(j::Int)`: converts ANSI standard indices to Fringe indices;
-* `standardize(fringe::Vector)`: formats a Fringe specified Zernike expansion coefficient vector according to the ANSI standard;
+* `standardize(fringe::Fringe)`: formats a Fringe specified Zernike expansion coefficient vector according to the ANSI standard;
 * `standardize(v_sub::Vector, orders::Vector{Tuple{Int, Int}})`: pads a subset Zernike expansion coefficient vector to the full standard length up to `n_max` (`1:j_max+1`).
 
+The `Noll` and `Fringe` types are used to wrap the input coefficient vectors for the `standardize` method arguments (e.g. `standardize(Fringe(v::Vector{Float64}))`). These can also be used to convert between the ordering schemes (e.g. `Noll(fringe::Fringe)`, `Fringe(s::Standard)`).
+
 The `standardize` fringe method expects unnormalized coefficients; the input coefficients will be re-ordered and normalized in line with the orthonormal standard. As Fringe is a 37 polynomial subset of the full set of Zernike polynomials any coefficients in the standard order missing a counterpart in the input vector will be set to zero.
 
 For the `standardize` subset method the tuples in `orders` must be of the form `(m, n)` associated with the respective coefficients at each index in `v_sub`.

src/Docstrings.jl

@@ -264,17 +264,6 @@ See also [`noll_to_j`](@ref), [`fringe_to_j`](@ref), [`j_to_fringe`](@ref), [`ge
 """
 j_to_noll
 
-"""
-    standardize!(noll::Vector)
-
-Re-order a Noll specified Zernike expansion coefficient vector according to the ANSI standard.
-
-This requires a full ordered vector up to `n_max`.
-
-See also [`standardize`](@ref), [`noll_to_j`](@ref), [`j_to_noll`](@ref), [`fringe_to_j`](@ref), [`j_to_fringe`](@ref), [`get_j`](@ref), [`get_mn`](@ref).
-"""
-standardize!
-
 """
     fringe_to_j(fringe::Int)
 
@@ -298,21 +287,31 @@ See also [`fringe_to_j`](@ref), [`noll_to_j`](@ref), [`j_to_noll`](@ref) [`stand
 j_to_fringe
 
 """
-    standardize(fringe::Vector)
+    standardize(noll::Noll)
+    standardize(fringe::Fringe)
+
+Format a `Noll` or `Fringe` specified Zernike expansion coefficient vector according to the ANSI standard.
 
-Format a Fringe specified Zernike expansion coefficient vector according to the ANSI standard.
+Floating-point coefficient vectors need to be wrapped in the index types (e.g. `standardize(Fringe(v))`).
 
-This function expects unnormalized coefficients; the input coefficients will be re-ordered and normalized in line with the orthonormal standard. As Fringe is a 37 polynomial subset of the full set of Zernike polynomials any coefficients in the standard order missing a counterpart in the input vector will be set to zero.
+The `Fringe` method expects unnormalized coefficients; the input coefficients will be re-ordered and normalized in line with the orthonormal standard. As Fringe is a 37 polynomial subset of the full set of Zernike polynomials any coefficients in the standard order missing a counterpart in the input vector will be set to zero.
 
 See also [`standardize!`](@ref), [`noll_to_j`](@ref), [`j_to_noll`](@ref), [`fringe_to_j`](@ref), [`j_to_fringe`](@ref), [`get_j`](@ref), [`get_mn`](@ref).
 
 ----
 
+    standardize(v_sub::FloatVec, j::AbstractVector{Int})
     standardize(v_sub::Vector, orders::Vector{Tuple{Int, Int}})
+    standardize(W::WavefrontError)
 
 Pad a subset Zernike expansion coefficient vector to the full standard length up to `n_max` (`1:j_max+1`).
 
 The tuples in `orders` must be of the form `(m, n)` associated with the respective coefficients at each index in `v_sub`.
+
+`j` is a vector of single-mode ordering indices associated with the coefficients.
+
+The `WavefrontError` method pads the `W.a` coefficient vector.
+
 """
 standardize
 

src/WavefrontError.jl

@@ -245,15 +245,18 @@ function sieve(v::Vector{Float64}, threshold::Float64)
 end
 
 # pads a subset Zernike expansion coefficient vector to standard length
-function standardize(v_sub::FloatVec, orders::Vector{Tuple{Int, Int}})
-    j = [get_j(mn...) for mn in orders]
+function standardize(v_sub::FloatVec, j::AbstractVector{Int})
     n_max = get_n(maximum(j))
     j_max = get_j(n_max)
     v_padded = zeros(eltype(v_sub), j_max + 1)
     v_padded[j.+1] .= v_sub
     return v_padded
 end
 
+function standardize(v_sub::FloatVec, orders::Vector{Tuple{Int, Int}})
+    standardize(v_sub, get_j.(orders))
+end
+
 standardize(W::WavefrontError) = standardize(W.a, [(i.m, i.n) for i ∈ W.recap])
 
 # methods

src/Zernike.jl

@@ -8,8 +8,8 @@
 module Zernike
 
 export zernike, wavefront, transform, Z, W, P, WavefrontError, get_j, get_mn,
-       noll_to_j, j_to_noll, fringe_to_j, j_to_fringe, standardize, standardize!,
-       Observable, plotconfig, zplot, reduce_wave
+       Noll, Fringe, noll_to_j, j_to_noll, fringe_to_j, j_to_fringe, standardize,
+       Standard, Observable, plotconfig, zplot, reduce_wave
 
 const public_names = "public \
     radial_coefficients, wavefront_coefficients, transform_coefficients, \
@@ -68,6 +68,46 @@ struct Output
     high_order::Bool
 end
 
+# normalization constant for Fringe coefficients
+N(j::Vector{Int}) = broadcast(x -> √radicand(x...), get_mn.(j))
+
+struct Standard
+    v::Vector{Float64}
+end
+
+struct Noll
+    v::Vector{Float64}
+end
+
+struct Fringe
+    v::Vector{Float64}
+end
+
+Standard(noll::Noll) = Standard(standardize(noll))
+Standard(fringe::Fringe) = Standard(standardize(fringe))
+
+Noll(fringe::Fringe) = (Noll ∘ Standard ∘ standardize)(fringe)
+
+Fringe(noll::Noll) = (Fringe ∘ Standard ∘ standardize)(noll)
+
+function Noll(s::Standard)
+    sv = s.v
+    j = eachindex(sv) .- 1
+    noll = j_to_noll.(j)
+    v = zeros(maximum(noll))
+    v[noll] = sv
+    return Noll(v)
+end
+
+function Fringe(s::Standard)
+    sv = s.v
+    j = [i for i ∈ 0:length(sv)-1 if i ∈ valid_fringes]
+    fringe = j_to_fringe.(j)
+    v = zeros(maximum(fringe))
+    v[fringe] .= N(j) .* sv[j.+1]
+    return Fringe(v)
+end
+
 macro domain(ex, msg::String, val)
     esc(:($ex ? nothing : throw(DomainError($val, $msg))))
 end
@@ -175,6 +215,8 @@ function get_j(m::Int, n::Int)
     return ((n + 2)n + m) ÷ 2
 end
 
+get_j((m, n)) = get_j(m, n)
+
 get_j(n_max::Int) = get_j(n_max, n_max)
 
 function get_mn(j::Int)
@@ -201,9 +243,12 @@ function j_to_noll(j::Int)
     (n + 1)n ÷ 2 + abs(m) + k
 end
 
-function standardize!(noll::Vector)
-    validate_length(noll)
-    invpermute!(noll, [noll_to_j(i) + 1 for i in eachindex(noll)])
+function standardize(noll::Noll)
+    (; v) = noll
+    p = [noll_to_j(i) + 1 for i in eachindex(v)]
+    u = zeros(maximum(p))
+    u[p] = v
+    return u
 end
 
 function fringe_to_j(fringe::Int)
@@ -229,21 +274,20 @@ function j_to_fringe(j::Int)
     @domain(fringe ∈ 1:36,
         """
         \nPolynomial does not have a Fringe representation.
-        Call fringe_to_j.(1:37) for valid polynomial indices.
+        See Zernike.valid_fringes for valid polynomial indices.
         """,
         j
     )
     return fringe
 end
 
-function standardize(fringe::Vector)
-    j = fringe_to_j.(eachindex(fringe))
+function standardize(fringe::Fringe)
+    (; v) = fringe
+    j = fringe_to_j.(eachindex(v))
     n_max = get_n(maximum(j))
     j_max = get_j(n_max)
-    a = zeros(eltype(fringe), j_max + 1)
-    # normalize
-    N = broadcast(x -> √radicand(x...), get_mn.(j))
-    a[j.+1] = fringe ./ N
+    a = zeros(eltype(v), j_max + 1)
+    a[j.+1] = v ./ N(j)
     return a
 end
 
@@ -361,6 +405,8 @@ Z(j::Int) = Z(get_mn(j)...)
 
 const piston = Z(0, 0)
 
+const valid_fringes = fringe_to_j.(1:37)
+
 # API namespace
 radial_coefficients(x...) = Φ(x...)[end]
 

test/runtests.jl

@@ -1,8 +1,8 @@
 using Test
 using Zernike
 using Zernike: radicand, Φ, get_i, canonical, coords, reconstruct, validate_length,
-               map_phase, format_strings, LaTeXString, latexstring, J, metrics,
-               polar, max_precision, Superposition, Product, sieve, (..)
+               valid_fringes, map_phase, format_strings, LaTeXString, latexstring, J,
+               metrics, polar, max_precision, Superposition, Product, sieve, (..)
 using StatsBase: mean, sample
 
 @testset "fringe" begin
@@ -11,12 +11,12 @@ using StatsBase: mean, sample
     @test fringe_to_j.(fringe) == 0:14
     v_fringe = collect(1.0:18.0)
     N_ = (√radicand(m, n) for n = 0:4 for m = -n:2:n)
-    compare = map(N_, standardize(v_fringe), v_fringe[fringe]) do N, v, o
+    compare = map(N_, standardize(Fringe(v_fringe)), v_fringe[fringe]) do N, v, o
                   N * v ≈ o
               end
     @test compare |> all
     @test_throws DomainError j_to_fringe(21) # Z(-6, 6) ∉ Fringe
-    @test j_to_fringe.(fringe_to_j.(1:37)) == 1:37
+    @test j_to_fringe.(valid_fringes) == 1:37
 end
 
 @testset "noll" begin
@@ -26,7 +26,7 @@ end
     v = zeros(15)
     v[noll] = 1:15
     @test !isequal(v, 1:15)
-    standardize!(v)
+    v = standardize(Noll(v))
     @test isequal(v, 1:15)
     @test noll_to_j.(j_to_noll.(0:1000)) == 0:1000
     j = get_j(101)
@@ -35,10 +35,20 @@ end
     w = copy(u)
     invpermute!(w, p)
     @test w != u
-    standardize!(w)
+    w = standardize(Noll(w))
     @test w == u
 end
 
+# more index conversions
+@testset "standard" begin
+    j = collect(0:101)
+    v = rand(102)
+    s = Standard(standardize(v, j))
+    i = valid_fringes .+ 1
+    @test standardize(Fringe(Noll(s)))[i] ≈ v[i]
+    @test standardize(Noll(Fringe(s)))[i] ≈ v[i]
+end
+
 function compare_coefficients(j_max)
     γ = Vector{Float64}[]
     zλ = Vector{Vector{Float64}}(undef, j_max + 1)