text:
# These are the functions that need to be imported so that new methods do not overwrite old ones
import Base.show
import Base.one
import Base.∘
import Base.inv
import Base.^
import Base.:(==)

"""
    Permutation(image::Vector{Int64})
    Permutation(xs...::Int64)

A function ℕ → ℕ which fixes all but finitely many points.

Construct a `Permutation` from a vector of integers as the image of the integers `1:length(image)`.

Alternatively, construct from an arbitrary list of integers.

# Examples
```julia-repl
julia> Permutation([1,5,3,2,4])
( 2 4 5 )

julia> Permutation(4,5,2,3,1)
( 1 5 2 3 4 )
```
"""
struct Permutation
    image::Vector{Int64}
    # An inner constructor checks that the permutation is valid
    function Permutation(image::Vector{Int64})
        m = maximum([image[image .!= 1:length(image)]..., 1])
        sort(image) == 1:length(image) && return new(image[1:m])
        error("not a valid permutation")
    end
end

Permutation(imagevals::Int64...) = Permutation([x for x ∈ imagevals])

"""
    maxarg(σ::Permutation)

Find the largest integer which the `Permutation` `σ` does not fix.

For the identity `Permutation`, return `1`.

# Examples
```julia-repl
julia> maxarg(Permutation(1,3,2,4))
3

julia> maxarg(ι)
1
```
"""
maxarg(σ::Permutation) = length(σ.image)

"""
    (σ::Permutation)(x::Int64)

Evaluate the action of the `Permutation` `σ` on `x`.

# Examples
julia> σ = Permutation(4,5,2,1,3)
( 1 4 )( 2 5 3 )

julia> σ(5)
3
"""
(σ::Permutation)(x::Int64) = x ∈ 1:maxarg(σ) ? σ.image[x] : x

"""
    orbit(σ::Permutation,x::Int64)

Compute the orbit of `x` under the repeated action of the `Permutation` `σ`.

# Examples
```julia-repl
julia> orbit(Permutation(2,3,4,5,1),3)
5-element Vector{Int64}:
 3
 4
 5
 1
 2
```
"""
function orbit(σ::Permutation, x::Int64)
    orb = [x]
    y = σ(x)
    while y != x
        push!(orb,y)
        y = σ(y)
    end
    return orb
end

"""
    dcd(σ::Permutation)

Compute the disjoint cycle decomposition of the `Permutation` `σ`.

If an integer is fixed by `σ`, then it is not included in the decomposition.

# Examples
```julia-repl
julia> dcd(Permutation(4,3,5,1,2))
2-element Vector{Vector{Int64}}:
 [1, 4]
 [2, 3, 5]

julia> dcd(Permutation(1,3,2))
1-element Vector{Vector{Int64}}:
 [2, 3]
```
"""
function dcd(σ::Permutation)
    decomp = Vector{Int64}[]
    unaccounted = trues(maxarg(σ))
    while any(unaccounted)
        x = findfirst(unaccounted)
        xorbit = orbit(σ,x)
        unaccounted[xorbit] .= false
        length(xorbit) > 1 && push!(decomp,xorbit)
    end
    return decomp
end

function show(io::IO,σ::Permutation)
    toprint = *([ *( "( ", ["$y " for y ∈ x]... , ")" ) for x ∈ dcd(σ)]..., "")
    toprint == "" && (toprint = "ι")
    print(io, toprint)
end

"""
    ι

The identity `Permutation`.

# Examples
```julia-repl
julia> ι.(1:10) == 1:10
true
```
"""
const ι = Permutation(1)
one(::Permutation) = ι
one(::Type{Permutation}) = ι

"""
    ∘(σ::Permutation)
    ∘(σ::Permutation, τ::Permutation)
    ∘(σs...::Permutation)

Compose `Permutation`s.

Compose a single `Permutation` and recieve it back.

# Examples
```julia-repl
julia> ∘(Permutation(1,2,4,3))
( 3 4 )

julia> Permutation(5,2,3,1,4) ∘ Permutation(2,1,4,5,3)
( 1 2 5 3 )
```
"""
∘(σ::Permutation) = σ
∘(σ::Permutation, τ::Permutation) = Permutation([σ(τ(x)) for x ∈ 1:max(maxarg(σ),maxarg(τ))])
∘(σ::Permutation, τ::Permutation, υ::Permutation...) = ∘(σ ∘ τ, υ...)

inv(σ::Permutation) = Permutation([findfirst(σ.image .== x) for x ∈ 1:maxarg(σ)])

"""
    ^(σ::Permutation,n::Int)

Exponentiate a `Permutation` `σ` by an integer `n`.

If `n` is negative, compute the inverse and exponentiate by `-n`.
```julia-repl
julia> Permutation(5,4,2,3,1)^5
( 1 5 )( 2 3 4 )

julia> Permutation(5,4,2,3,1)^-5
( 1 5 )( 2 4 3 )
```
"""
function ^(σ::Permutation,n::Int64)
    n == 0 && return one(Permutation)
    n < 0 && ((n,σ) = (-n,inv(σ)))
    return ∘(fill(σ,n)...)
end

"""
    transposition(m::Int64,n::Int64)

Create the `Permutation` which is a transposition of `m` and `n`.

If `m == n`, the `Permutation` created is the identity.

# Examples
```julia-repl
julia> transposition(4,5)
( 4 5 )

julia> transposition(6,6)
ι
```
"""
function transposition(m::Int64, n::Int64)
    image = collect(1:max(m,n))
    image[m], image[n] = n, m
    return Permutation(image)
end

"""
    symmetricgroup(n::Int64)

Compute the `Permutation`s which make up the symmetric group of order `factorial(n)`.

# Examples
```julia-repl
julia> symmetricgroup(5)
120-element Vector{Permutation}:
[...]
```
"""
function symmetricgroup(n::Int64)
    n ≥ 1 || error("symmetric group must have a positive parameter")
    n == 1 && return [ι]

    permutations = vcat(fill(symmetricgroup(n-1),n)...)
    transpositions = vcat([fill(transposition(i,n),factorial(n-1)) for i ∈ 1:n]...)
    return permutations .∘ transpositions
end

"""
    parity(σ::Permutation)

Compute the parity of the `Permutation` `σ`.

# Examples
```julia-repl
julia> parity(Permutation(1,2,3,5,4))
-1

julia> parity(ι)
1
```
"""
parity(σ::Permutation) = (-1)^sum([iseven(length(x)) for x ∈ dcd(σ)])

"""
    alternatinggroup(n::Int64)

Compute the `Permutation`s which make up the `n`th alternating group.

# Examples
```julia-repl
julia> alternatinggroup(5)
60-element Vector{Permutation}:
[...]
```
"""
function alternatinggroup(n::Int64)
    Sₙ = symmetricgroup(n)
    return Sₙ[parity.(Sₙ) .== 1]
end

# Defines equality of permutations, to be used to check ∉
==(σ::Permutation, τ::Permutation) = (σ.image == τ.image)

"""
    generate(σs::Permutation...)

Compute the elements of the group generated by the elements listed.

# Examples
```julia-repl
julia> generate(Permutation(2,3,4,1))
4-element Vector{Permutation}:
 ( 1 2 3 4 )
 ( 1 3 )( 2 4 )
 ( 1 4 3 2 )
 ι

 julia> generate(Permutation(2,3,4,5,1),Permutation(5,4,3,2,1))
 10-element Vector{Permutation}:
  ( 1 2 3 4 5 )
  ( 1 5 )( 2 4 )
  ( 1 3 5 2 4 )
  ( 2 5 )( 3 4 )
  ( 1 5 4 3 2 )
  ( 1 4 )( 2 3 )
  ι
  ( 1 4 2 5 3 )
  ( 1 2 )( 3 5 )
  ( 1 3 )( 4 5 )
```
"""
function generate(σs::Permutation...)
    G = [σs...]
    while true
        newelements = Permutation[]

        for σ ∈ G, τ ∈ G
            υ = σ ∘ τ
            υ ∉ G && υ ∉ newelements && push!(newelements, υ)
        end

        if isempty(newelements)
            return G
        else
            G = vcat(G,newelements)
        end
    end
end


text:
using SimpleDelimitedFiles
using Test
using DelimitedFiles
using Aqua
using Printf

@testset "project quality" begin
    Aqua.test_all(SimpleDelimitedFiles)
end

@testset "Comparison with DelimitedFiles" begin
    mktemp() do fname, f
        A = rand(2, 2)
        writedlm(fname, A)
        @test readdlm(fname) == SimpleDelimitedFiles.readdlm(fname)
        A = rand(2, 2) * 2 .- 1
        writedlm(fname, A)
        @test readdlm(fname) == SimpleDelimitedFiles.readdlm(fname)
        @testset "single column" begin
            s = join([(@sprintf "%10.2f" i) for i in 1:4], '\n') * "\n"
            open(fname, "w") do f
                write(f, s)
            end
            @test readdlm(fname) == SimpleDelimitedFiles.readdlm(fname)
            # leading space
            s = join(["  " * (@sprintf "%10.2f" i) for i in 1:4], '\n') * "\n"
            open(fname, "w") do f
                write(f, s)
            end
            @test readdlm(fname) == SimpleDelimitedFiles.readdlm(fname)
            s = join([(s = @sprintf "%10.2f" i; "  $s $s") for i in 1:4], '\n') * "\n"
            open(fname, "w") do f
                write(f, s)
            end
            M = readdlm(fname, ' ')
            M = M[:, (x->all(y -> y isa Number, x)).(eachcol(M))]
            @test M  == SimpleDelimitedFiles.readdlm(fname, ' ')
        end
        for T in [Float64, Int]
            A = rand(T, 2, 2)
            writedlm(fname, A)
            @test readdlm(fname, T) == SimpleDelimitedFiles.readdlm(fname, T)
            for delim in [',', '\t', ' ']
                writedlm(fname, A, delim)
                @test readdlm(fname, delim, T) == SimpleDelimitedFiles.readdlm(fname, delim, T)
            end
        end
    end
end

ext: jl


text:
module Postprocessing

export soltoResult,soltoMatrxResult

using ..Systems,..Tools

"""
    reserved for future use
"""

function soltoResult(sol,sys0) #only good for one calculation per time point, not good for onec calculation for all time

    γ = sys0.liquid.γ
    numofliquidslug =  Integer( (size(sol)[1]-1)/5  )

    MatrxXp=sol[1:2*numofliquidslug,:]
    MatrxdXdt=sol[2*numofliquidslug+1:4*numofliquidslug,:]
    M=sol[4*numofliquidslug+1:5*numofliquidslug+1,:]

    Xp = Array{Tuple{Float64, Float64}}((undef), Integer(size(MatrxXp)[1]/2), size(MatrxXp)[2])
    dXdt = Array{Tuple{Float64, Float64}}((undef), size(Xp))
    Lvaporplug = zeros(size(M))

    # transfer from 2D array to array of tuple for Xp and dXdt
    for i in 1:size(Xp)[1]
        for j in 1:size(Xp)[2]
            Xp[i,j]=((MatrxXp[2*i-1,j]),(MatrxXp[2*i,j]))
            dXdt[i,j]=((MatrxdXdt[2*i-1,j]),(MatrxdXdt[2*i,j]))
        end
    end

    for j in 1:size(M)[2]
        Lvaporplug[:,j]=XptoLvaporplug(Xp[:,j], sys0.tube.L, sys0.tube.closedornot)
    end

    # get P from M and γ
    # P = (M./Lvaporplug).^(γ)
    # P = zeros(size(M))
    # for i in length(M)
    #        P[i] = M[i] > 0 ? (M[i]./Lvaporplug[i]).^(γ) : -(-M[i]./Lvaporplug[i]).^(γ)
    #    end
    P = real.((M./Lvaporplug .+ 0im).^(γ))

    # get θ from P and γ
    # θ = zeros(size(P))
    # for i in length(P)
    #        θ[i] = P[i] > 0 ? P[i].^((γ-1)/γ) : -(-P[i]).^((γ-1)/γ)
    #    end
    θ = real.((P .+ 0im).^((γ-1)/γ))
    # θ = P.^((γ-1)/γ)

    # convert a matrix to array of array
    Xp=mapslices(x->[x], MatrxXp, dims=2)[:]
    dXdt=mapslices(x->[x], MatrxdXdt, dims=2)[:]
    M=mapslices(x->[x], M, dims=2)[:]
    P=mapslices(x->[x], P, dims=2)[:]
    θ=mapslices(x->[x], θ, dims=2)[:]

    # if the input sol is merely a state vector u Array{Float64,1} rather than ODESolution{}
    # the result then does not contain t (we set the t to be -1.0)
    if typeof(sol) == Array{Float64,1}
        return PHPResult(-1.0, vec(hcat(Xp...)),vec(hcat(dXdt...)),vec(hcat(P...)),vec(hcat(θ...)),vec(hcat(M...)))
    elseif typeof(sol) == Array{Float64,2}
        return PHPResult(-1.0,Xp,dXdt,P,θ,M)
    else
        return PHPResult(sol.t,Xp,dXdt,P,θ,M)
    end
    # result=typeof(sol) == Array{Float64,1} ? PHPResult(-1.0, vec(hcat(Xp...)),vec(hcat(dXdt...)),vec(hcat(P...)),vec(hcat(θ...)),vec(hcat(M...))) : PHPResult(sol.t,Xp,dXdt,P,θ,M)
end

function soltoMatrxResult(sol,sys0) #only good for one calculation per time point, not good for onec calculation for all time

    indexes = Int64[]
    θliquidrec = Array[]

    for i = 1:length(sol[1])
        if abs(sol[1][i]+1e10) <= 10^(-1)
            push!(indexes,i)
        end
    end

    γ = sys0.liquid.γ
    # numofliquidslug =  Integer( (size(sol)[1]-1)/5  )
    # numofliquidslug =  Integer( (indexes[1]-2)/5  )
    numofliquidslug =  Integer( (indexes[1]-3)/6  )

    MatrxXp=sol[1:2*numofliquidslug,:]
    MatrxdXdt=sol[2*numofliquidslug+1:4*numofliquidslug,:]
    MatrxM=sol[4*numofliquidslug+1:5*numofliquidslug+1,:]
    Matrxδ=sol[5*numofliquidslug+2:6*numofliquidslug+2,:]

    return MatrxXp, MatrxdXdt, MatrxM, Matrxδ
end

end


text:
#Demo 5: Fourier Series
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))
magvsf = cons(:a, labels = set(yaxis="|X(f)|²", xaxis="Frequency (Hz)"))

#==Input data
===============================================================================#
tfund = 1e-9 #Fundamental
osr = 20 #(fundamental) oversampling ratio
nfund = 20 #cycles of the fundamenal

#==Computations
===============================================================================#
t = DataTime(timespace(:ts, tfund/osr, :tfund, nfund*tfund), tperiodic=true)
sig = sin(t*(2pi/tfund))
Fsig = freqdomain(sig)
#sigSpec = abs2(fspectrum(Fsig)) #Aperiodic
sigSpec = abs2(fcoeff(Fsig)) #Periodic

#==Generate plot
===============================================================================#
plot_t = push!(cons(:plot, vvst, title="Time domain", legend=false),
#	cons(:wfrm, DataF1(sig), glyph=set(shape=:o)),
	cons(:wfrm, DataF1(sig)),
)
plot_f = push!(cons(:plot, magvsf, title="Frequency domain", legend=false),
	cons(:wfrm, sigSpec, glyph=set(shape=:o, size=1.5), line=set(style=:none, color=:red, width=2)),
)

pcoll = push!(cons(:plotcoll, title="Fourier Transform Tests"), plot_t, plot_f)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
export iauTrxpv
"""
Multiply a pv-vector by the transpose of an r-matrix.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   r        double[3][3]    r-matrix
   pv       double[2][3]    pv-vector

Returned:
   trpv     double[2][3]    r * pv

Note:
   It is permissible for pv and trpv to be the same array.

Called:
   iauTr        transpose r-matrix
   iauRxpv      product of r-matrix and pv-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauTrxpv(double r[3][3], double pv[2][3], double trpv[2][3])

function iauTrxpv(r::AbstractMatrix{<:Real}, pv::AbstractMatrix{<:Real})
   trpv = zeros(Float64, 3, 2)

   ccall((:iauTrxpv, libsofa_c), Cvoid,
       (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Array{Float64, 2}, r'),
       convert(Array{Float64, 2}, pv'),
       trpv)

   return SMatrix{2,3}(trpv')
end


text:
abstract type AbstractAnnData end

mutable struct AnnData <: AbstractAnnData
    file::Union{HDF5.File, HDF5.Group, Nothing}

    X::Union{AbstractMatrix{<:Number}, Nothing}

    obs::DataFrame
    obs_names::Index{<:AbstractString}

    var::DataFrame
    var_names::Index{<:AbstractString}

    obsm::StrAlignedMapping{Tuple{1 => 1}, AnnData}
    obsp::StrAlignedMapping{Tuple{1 => 1, 2 => 1}, AnnData}

    varm::StrAlignedMapping{Tuple{1 => 2}, AnnData}
    varp::StrAlignedMapping{Tuple{1 => 2, 2 => 2}, AnnData}

    layers::AbstractAlignedMapping{Tuple{1 => 1, 2 => 2}, String}

    uns::Dict{<:AbstractString, <:Any}

    function AnnData(file::Union{HDF5.File, HDF5.Group}, backed=true, checkversion=true)
        if checkversion
            attrs = attributes(file)
            if !haskey(attrs, "encoding-type")
                @warn "The HDF5 file was not created by muon, we can't guarantee that everything will work correctly"
            elseif attrs["encoding-type"] != "AnnData"
                error("This HDF5 file does not appear to hold an AnnData object")
            end
        end

        adata = new(backed ? file : nothing)

        # this needs to go first because it's used by size() and size()
        # is used for dimensionalty checks
        adata.obs, obs_names = read_dataframe(file["obs"])
        adata.var, var_names = read_dataframe(file["var"])
        adata.obs_names = Index(obs_names)
        adata.var_names = Index(var_names)

        # observations
        adata.obsm = StrAlignedMapping{Tuple{1 => 1}}(
            adata,
            haskey(file, "obsm") ? read_dict_of_mixed(file["obsm"]) : nothing,
        )
        adata.obsp = StrAlignedMapping{Tuple{1 => 1, 2 => 1}}(
            adata,
            haskey(file, "obsp") ? read_dict_of_matrices(file["obsp"]) : nothing,
        )

        # Variables
        adata.varm = StrAlignedMapping{Tuple{1 => 2}}(
            adata,
            haskey(file, "varm") ? read_dict_of_mixed(file["varm"]) : nothing,
        )
        adata.varp = StrAlignedMapping{Tuple{1 => 2, 2 => 2}}(
            adata,
            haskey(file, "varp") ? read_dict_of_matrices(file["varp"]) : nothing,
        )

        # X
        adata.X = backed ? nothing : read_matrix(file["X"])

        # Layers
        if !backed
            adata.layers = StrAlignedMapping{Tuple{1 => 1, 2 => 2}}(
                adata,
                haskey(file, "layers") ? read_dict_of_matrices(file["layers"]) : nothing,
            )
        else
            adata.layers = BackedAlignedMapping{Tuple{1 => 1, 2 => 2}}(adata, adata.file, "layers")
        end

        # unstructured
        adata.uns =
            haskey(file, "uns") ? read_dict_of_mixed(file["uns"], separate_index=false) :
            Dict{String, Any}()

        return adata
    end

    function AnnData(;
        X::AbstractMatrix{<:Number},
        obs::Union{DataFrame, Nothing}=nothing,
        obs_names::Union{AbstractVector{<:AbstractString}, Nothing}=nothing,
        var::Union{DataFrame, Nothing}=nothing,
        var_names::Union{AbstractVector{<:AbstractString}, Nothing}=nothing,
        obsm::Union{
            AbstractDict{<:AbstractString, <:Union{<:AbstractArray{<:Number}, DataFrame}},
            Nothing,
        }=nothing,
        varm::Union{
            AbstractDict{<:AbstractString, <:Union{<:AbstractArray{<:Number}, DataFrame}},
            Nothing,
        }=nothing,
        obsp::Union{AbstractDict{<:AbstractString, <:AbstractMatrix{<:Number}}, Nothing}=nothing,
        varp::Union{AbstractDict{<:AbstractString, <:AbstractMatrix{<:Number}}, Nothing}=nothing,
        layers::Union{AbstractDict{<:AbstractString, <:AbstractMatrix{<:Number}}, Nothing}=nothing,
        uns::Union{AbstractDict{<:AbstractString, <:Any}, Nothing}=nothing,
    )
        m, n = size(X)
        if isnothing(obs)
            obs = DataFrame()
        elseif size(obs, 1) != m
            throw(DimensionMismatch("X has $m rows, but obs has $(size(obs, 1)) rows"))
        end

        if isnothing(obs_names)
            obs_names = string.(collect(1:m))
        elseif length(obs_names) != m
            throw(DimensionMismatch("X has $m rows, but $(length(obs_names)) obs_names given"))
        end

        if isnothing(var)
            var = DataFrame()
        elseif size(var, 1) != n
            throw(DimensionMismatch("X has $n columns, but var has $(size(var, 1)) rows"))
        end

        if isnothing(var_names)
            var_names = string.(collect(1:n))
        elseif length(var_names) != n
            throw(DimensionMismatch("X has $n columns, but $(length(var_names)) var_names given"))
        end
        adata = new(nothing, X, obs, Index(obs_names), var, Index(var_names))
        adata.obsm = StrAlignedMapping{Tuple{1 => 1}}(adata, obsm)
        adata.obsp = StrAlignedMapping{Tuple{1 => 1, 2 => 1}}(adata, obsp)
        adata.varm = StrAlignedMapping{Tuple{1 => 2}}(adata, varm)
        adata.varp = StrAlignedMapping{Tuple{1 => 2, 2 => 2}}(adata, varp)
        adata.layers = StrAlignedMapping{Tuple{1 => 1, 2 => 2}}(adata, layers)
        adata.uns = isnothing(uns) ? Dict{String, Any}() : uns

        return adata
    end
end

file(ad::AnnData) = ad.file

function readh5ad(filename::AbstractString; backed=true)
    filename = abspath(filename) # this gets stored in the HDF5 objects for backed datasets
    if !backed
        fid = h5open(filename, "r")
    else
        fid = h5open(filename, "r+")
    end
    local adata
    try
        adata = AnnData(fid, backed, false)
    catch e
        close(fid)
        rethrow()
    end
    if !backed
        close(fid)
    end
    return adata
end

function writeh5ad(filename::AbstractString, adata::AbstractAnnData)
    filename = abspath(filename)
    if file(adata) === nothing || filename != HDF5.filename(file(adata))
        hfile = h5open(filename, "w", userblock=512)
        try
            write(hfile, adata)
            close(hfile)
            hfile = open(filename, "r+")
            write(
                hfile,
                "AnnData (format-version=$ANNDATAVERSION;creator=$NAME;creator-version=$VERSION)",
            )
        finally
            close(hfile)
        end
    else
        write(adata)
    end
    return nothing
end

function Base.write(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    adata::AbstractAnnData,
)
    g = create_group(parent, name)
    write(g, adata)
end

function Base.write(parent::Union{HDF5.File, HDF5.Group}, adata::AbstractAnnData)
    attrs = attributes(parent)
    attrs["encoding-type"] = "AnnData"
    attrs["encoding-version"] = string(ANNDATAVERSION)
    attrs["encoder"] = NAME
    attrs["encoder-version"] = string(VERSION)
    if parent === file(adata)
        write(adata)
    else
        write_attr(parent, "X", adata.X)
        write_attr(parent, "layers", adata.layers)
        write_metadata(parent, adata)
    end
end

function Base.write(adata)
    if file(adata) === nothing
        error("adata is not backed, need somewhere to write to")
    end
    write_metadata(file(adata), adata)
end

function write_metadata(parent::Union{HDF5.File, HDF5.Group}, adata::AbstractAnnData)
    write_attr(parent, "obs", adata.obs, index=adata.obs_names)
    write_attr(parent, "obsm", adata.obsm, index=adata.obs_names)
    write_attr(parent, "obsp", adata.obsp)
    write_attr(parent, "var", adata.var, index=adata.var_names)
    write_attr(parent, "varm", adata.varm, index=adata.var_names)
    write_attr(parent, "varp", adata.varp)
    write_attr(parent, "uns", adata.uns)
end
# FileIO support
load(f::File{format"h5ad"}) = readh5ad(filename(f), backed=false) # I suppose this is more consistent with the rest of FileIO?
save(f::File{format"h5ad"}, data::AbstractAnnData) = writeh5ad(filename(f), data)

Base.size(adata::AbstractAnnData) = (length(adata.obs_names), length(adata.var_names))
Base.size(adata::AbstractAnnData, d::Integer) = size(adata)[d]

function Base.show(io::IO, adata::AbstractAnnData)
    compact = get(io, :compact, false)
    print(io, """$(typeof(adata).name.name) object $(size(adata)[1]) \u2715 $(size(adata)[2])""")
end

function Base.show(io::IO, ::MIME"text/plain", adata::AbstractAnnData)
    show(io, adata)
end

function Base.getproperty(adata::AnnData, s::Symbol)
    if s === :X && isbacked(adata)
        return backed_matrix(file(adata)["X"])
    else
        return getfield(adata, s)
    end
end

function Base.getindex(
    adata::AbstractAnnData,
    I::Union{
        AbstractUnitRange,
        Colon,
        AbstractVector{<:Integer},
        AbstractVector{<:AbstractString},
        Integer,
        AbstractString,
    },
    J::Union{
        AbstractUnitRange,
        Colon,
        AbstractVector{<:Integer},
        AbstractVector{<:AbstractString},
        Integer,
        AbstractString,
    },
)
    @boundscheck checkbounds(adata, I, J)
    i, j = convertidx(I, adata.obs_names), convertidx(J, adata.var_names)
    newad = AnnData(
        X=adata.X[i, j],
        obs=isempty(adata.obs) ? nothing : adata.obs[i, :],
        obs_names=adata.obs_names[i],
        var=isempty(adata.var) ? nothing : adata.var[j, :],
        var_names=adata.var_names[j],
    )
    copy_subset(adata.obsm, newad.obsm, i, j)
    copy_subset(adata.varm, newad.varm, i, j)
    copy_subset(adata.obsp, newad.obsp, i, j)
    copy_subset(adata.varp, newad.varp, j, j)
    copy_subset(adata.layers, newad.layers, i, j)
    return newad
end

struct AnnDataView{Ti, Tj} <: AbstractAnnData
    parent::AnnData
    I::Ti
    J::Tj

    X::Union{AbstractMatrix{<:Number}, Nothing}

    obs::SubDataFrame
    obs_names::SubIndex{<:AbstractString}

    var::SubDataFrame
    var_names::SubIndex{<:AbstractString}

    obsm::StrAlignedMappingView{Tuple{1 => 1}}
    obsp::StrAlignedMappingView{Tuple{1 => 1, 2 => 1}}

    varm::StrAlignedMappingView{Tuple{1 => 2}}
    varp::StrAlignedMappingView{Tuple{1 => 2, 2 => 2}}

    layers::StrAlignedMappingView{Tuple{1 => 1, 2 => 2}}

    uns::Dict{<:AbstractString, <:Any}
end

function Base.view(ad::AnnData, I, J)
    @boundscheck checkbounds(ad, I, J)
    i, j = convertidx(I, ad.obs_names), convertidx(J, ad.var_names)
    X = isbacked(ad) ? nothing : @view ad.X[i, j]

    return AnnDataView(
        ad,
        i,
        j,
        X,
        view(ad.obs, nrow(ad.obs) > 0 ? i : (:), :),
        view(ad.obs_names, i),
        view(ad.var, nrow(ad.var) > 0 ? j : (:), :),
        view(ad.var_names, j),
        view(ad.obsm, i),
        view(ad.obsp, i, i),
        view(ad.varm, j),
        view(ad.varp, j, j),
        view(ad.layers, i, j),
        ad.uns,
    )
end
function Base.view(ad::AnnDataView, I, J)
    @boundscheck checkbounds(ad, I, J)
    i, j =
        Base.reindex(parentindices(ad), (convertidx(I, ad.obs_names), convertidx(J, ad.var_names)))
    return view(parent(ad), i, j)
end

function Base.getproperty(ad::AnnDataView, s::Symbol)
    if s === :X && isbacked(ad)
        return @view parent(ad).X[parentindices(ad)...]
    else
        return getfield(ad, s)
    end
end

Base.parent(ad::AnnData) = ad
Base.parent(ad::AnnDataView) = ad.parent
Base.parentindices(ad::AnnData) = axes(ad)
Base.parentindices(ad::AnnDataView) = (ad.I, ad.J)
file(ad::AnnDataView) = file(parent(ad))


text:
using Documenter
using GillMatsuno

makedocs(
    modules  = [GillMatsuno],
    doctest  = false,
    format   = Documenter.HTML(
        collapselevel = 1,
        prettyurls    = false
    ),
    authors  = "Nathanael Wong",
    sitename = "GillMatsuno.jl",
    pages    = [
        "Home"     => "index.md",
        # "Theory"   => "theory.md",
        # "Examples" => "examples.md"
    ]
)

deploydocs(
    repo = "github.com/natgeo-wong/GillMatsuno.jl.git",
)


text:
module DynamicGrids
# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end DynamicGrids

using Adapt,
      Colors,
      ConstructionBase,
      Crayons,
      DimensionalData,
      FreeTypeAbstraction,
      FileIO,
      LinearAlgebra,
      KernelAbstractions,
      OffsetArrays,
      REPL,
      Reexport,
      Requires,
      Setfield,
      StaticArrays,
      Test,
      UnicodeGraphics

@reexport using ModelParameters

const DG = DynamicGrids
const DD = DimensionalData

using Base: tail, @propagate_inbounds

import Base: show, getindex, setindex!, lastindex, size, length, push!, append!,
             broadcast, broadcast!, similar, eltype, iterate

export sim!, resume!, step!, savegif, isinferred

export rules, neighbors, neighborhood, kernel, kernelproduct, offsets, positions, radius

export inbounds, isinbounds

export gridsize, currentframe, currenttime, currenttimestep, timestep

export add!, sub!, min!, max!, and!, or!, xor!

export Rule, NeighborhoodRule, CellRule, SetCellRule, SetNeighborhoodRule, SetGridRule

export Cell, Neighbors, SetNeighbors, SetCell, Convolution, SetGrid, Life, CopyTo
export RuleWrapper, Chain, RunIf, RunAt
export AbstractRuleset, Ruleset, StaticRuleset

export Neighborhood, RadialNeighborhood, Window, AbstractKernelNeighborhood, Kernel,
       Moore, AbstractPositionalNeighborhood, Positional, VonNeumann, LayeredPositional

export Processor, SingleCPU, ThreadedCPU, CPUGPU

export PerformanceOpt, NoOpt, SparseOpt

export BoundaryCondition, Remove, Wrap

export ParameterSource, Aux, Grid, Delay, Lag, Frame

export Output, ArrayOutput, ResultOutput, TransformedOutput

export GraphicOutput, REPLOutput

export ImageOutput, GifOutput

export Renderer, Image, Layout, SparseOptInspector

export TextConfig

export ObjectScheme, Greyscale, Grayscale

export CharStyle, Block, Braile

const DEFAULT_KEY = :_default_

const EXPERIMENTAL = """
        WARNING: This feature is experimental. It may change in future versions,
        and may not be 100% reliable in all cases. Please file github issues if problems occur.
        """

function __init__()
    global terminal
    terminal = REPL.Terminals.TTYTerminal(get(ENV, "TERM", Base.Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)

    @require CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" include("cuda.jl")
end

include("interface.jl")
include("flags.jl")
include("neighborhoods.jl")
include("rules.jl")
include("settings.jl")
include("rulesets.jl")
include("extent.jl")
include("grid.jl")
include("simulationdata.jl")
include("parametersources.jl")
include("gpu.jl")
include("atomic.jl")
include("chain.jl")
include("condition.jl")
include("outputs/interface.jl")
include("outputs/output.jl")
include("outputs/graphic.jl")
include("outputs/image.jl")
include("outputs/textconfig.jl")
include("outputs/schemes.jl")
include("outputs/render.jl")
include("outputs/array.jl")
include("outputs/transformed.jl")
include("outputs/repl.jl")
include("outputs/gif.jl")
include("framework.jl")
include("modifyrule.jl")
include("sequencerules.jl")
include("generated.jl")
include("maprules.jl")
include("boundaries.jl")
include("utils.jl")
include("copyto.jl")
include("life.jl")
include("show.jl")

end


text:
const BlockDiagonal{T, TM} = BlockMatrix{T, <:Diagonal{TM}} where {TM <: AbstractMatrix{T}}

#
# Constructors
#

block_diagonal(vs::AbstractVector{<:AbstractMatrix}) = mortar(Diagonal(vs))

#
# Accumulation rule for Zygote.
#

function Zygote.accum(A::BlockDiagonal, B::BlockDiagonal)
    return block_diagonal(accum.(A.blocks.diag, B.blocks.diag))
end

#
# adjoint / transpose - ensure we get a BlockDiagonal back
#

LinearAlgebra.adjoint(A::BlockDiagonal) = block_diagonal(adjoint.(A.blocks.diag))
LinearAlgebra.transpose(A::BlockDiagonal) = block_diagonal(transpose.(A.blocks.diag))

#
# UpperTriangular - ensure we get a BlockDiagonal back
#

function LinearAlgebra.UpperTriangular(A::BlockDiagonal)
    return block_diagonal(UpperTriangular.(A.blocks.diag))
end

#
# Symmetric - ensure we get a BlockDiagonal back
#

LinearAlgebra.Symmetric(A::BlockDiagonal) = block_diagonal(Symmetric.(A.blocks.diag))

ZygoteRules.@adjoint function LinearAlgebra.Symmetric(A::BlockDiagonal)
    return Zygote.pullback(A->block_diagonal(Symmetric.(A.blocks.diag)), A)
end

#
# Addition
#

function Base.:+(A::BlockDiagonal, B::BlockDiagonal)
    return block_diagonal([a + b for (a, b) in zip(A.blocks.diag, B.blocks.diag)])
end

function Base.:+(A::Matrix, B::BlockDiagonal)
    @assert size(A) == size(B)
    C = copy(A)
    cs = cumulsizes(B, 1)
    for n in 1:nblocks(B, 1)
        idx = cs[n]:cs[n+1]-1
        C[idx, idx] += B[Block(n, n)]
    end
    return C
end

ZygoteRules.@adjoint function Base.:+(A::Matrix, B::BlockDiagonal{T, <:Matrix{T}} where {T})
    return A + B, function(Δ)
        cs = cumulsizes(B, 1)
        blks = [Δ[cs[n]:cs[n+1]-1, cs[n]:cs[n+1]-1] for n in 1:nblocks(B, 1)]
        return (Δ, block_diagonal(blks))
    end
end

#
# Negation
#

Base.:-(A::BlockDiagonal) = block_diagonal([-a for a in A.blocks.diag])

#
# BlockDiagonal multiplication
#

function Base.:*(A::BlockDiagonal{<:Real}, B::BlockDiagonal{<:Real})
    return block_diagonal([a * b for (a, b) in zip(A.blocks.diag, B.blocks.diag)])
end

ZygoteRules.@adjoint function Base.:*(A::BlockDiagonal{<:Real}, B::BlockDiagonal{<:Real})
    return A * B, Δ->(Δ * B', A' * Δ)
end

function Base.:*(A::BlockDiagonal{<:Real}, B::Matrix{<:Real})
    A_blks, B_blks = A.blocks.diag, BlockArray(B, blocksizes(A, 1), [size(B, 2)]).blocks
    return Matrix(mortar(reshape([a * b for (a, b) in zip(A_blks, B_blks)], :, 1)))
end
function Base.:*(A::BlockDiagonal{<:Real}, x::Vector{<:Real})
    A_blks, x_blks = diag(A.blocks), BlockArray(x, blocksizes(A, 1)).blocks
    return Vector(mortar([a * x for (a, x) in zip(A_blks, x_blks)]))
end

#
# BlockDiagonal ldiv
#

function Base.:\(A::BlockDiagonal{<:Real}, B::BlockDiagonal{<:Real})
    A_blks, B_blks = diag(A.blocks), diag(B.blocks)
    return block_diagonal([a \ b for (a, b) in zip(A_blks, B_blks)])
end

ZygoteRules.@adjoint function Base.:\(A::BlockDiagonal{<:Real}, B::BlockDiagonal{<:Real})
    Y = A \ B
    return Y, function(Ȳ::BlockDiagonal)
        B̄ = A' \ Ȳ
        return (-B̄ * Y', B̄)
    end
end

function Base.:\(A::BlockDiagonal{<:Real}, B::AbstractMatrix{<:Real})
    A_blks = diag(A.blocks)
    B_blks = BlockArray(collect(B), blocksizes(A, 1), [size(B, 2)]).blocks

    return Matrix(mortar(reshape([a \ b for (a, b) in zip(A_blks, B_blks)], :, 1)))
end
ZygoteRules.@adjoint function Base.:\(A::BlockDiagonal{<:Real}, B::AbstractMatrix{<:Real})
    Y = A \ B
    return Y, function(Ȳ::AbstractMatrix{<:Real})
        B̄ = A' \ Ȳ
        return (_block_diag_bit(-B̄, Y', A), B̄)
    end
end

function Base.:\(A::BlockDiagonal{<:Real}, x::AbstractVector{<:Real})
    return reshape(A \ reshape(x, :, 1), :)
end

ZygoteRules.@adjoint function Base.:\(A::BlockDiagonal{<:Real}, x::AbstractVector{<:Real})
    y_mat, back = Zygote.pullback(\, A, reshape(x, :, 1))
    return vec(y_mat), function(Δ::AbstractVector{<:Real})
        Ā, x̄ = back(reshape(Δ, :, 1))
        return Ā, vec(x̄)
    end
end

function _block_diag_bit(A::AbstractMatrix, B::AbstractMatrix, R::BlockDiagonal)
    A_blks = vec(BlockArray(A, blocksizes(R, 1), [size(A, 2)]).blocks)
    B_blks = vec(BlockArray(B', blocksizes(R, 1), [size(B, 1)]).blocks)
    return block_diagonal([a * b' for (a, b) in zip(A_blks, B_blks)])
end

#
# cholesky
#

function LinearAlgebra.cholesky(A::BlockDiagonal{<:Real})
    Cs = map(A->cholesky(A).U, diag(A.blocks))
    return Cholesky(block_diagonal(Cs), :U, 0)
end
ZygoteRules.@adjoint function LinearAlgebra.cholesky(A::BlockDiagonal{<:Real})
    Cs_backs = map(A->Zygote.pullback(A->cholesky(A).U, A), diag(A.blocks))
    Cs, backs = first.(Cs_backs), last.(Cs_backs)
    function back(Ū::BlockDiagonal)
        return (block_diagonal(map((Ū, back)->first(back(Ū)), diag(Ū.blocks), backs)),)
    end
    return Cholesky(block_diagonal(Cs), :U, 0), Δ->back(Δ.factors)
end

function LinearAlgebra.logdet(C::Cholesky{T, <:BlockDiagonal{T}} where {T<:Real})
    return 2 * sum([logdet(c) for c in C.factors.blocks.diag])
end

ZygoteRules.@adjoint function LinearAlgebra.logdet(
    C::Cholesky{T, <:BlockDiagonal{T}} where {T<:Real},
)
    return logdet(C), function(Δ::Real)
        blks = C.factors.blocks.diag
        factors = block_diagonal([diagm(0=>2Δ ./ diag(b)) for b in blks])
        return ((factors=factors,),)
    end
end

#
# Misc
#

#
# BlockDiagonal mul! and *
#

Base.:*(D::BlockDiagonal, x::BlockVector) = mul!(copy(x), D, x)
Base.:*(D::BlockDiagonal, X::BlockMatrix) = mul!(copy(X), D, X)

function LinearAlgebra.mul!(y::BlockVector, D::BlockDiagonal, x::BlockVector)
    @assert are_conformal(D, x) && are_conformal(D, y)
    blocks = D.blocks.diag
    for r in 1:nblocks(D, 1)
        mul!(view(y, Block(r)), blocks[r], view(x, Block(r)))
    end
    return y
end

@adjoint function Base.:*(D::BlockDiagonal, x::BlockVector)
    y = D * x
    return y, function(ȳ::BlockVector)
        @assert blocksizes(y, 1) == blocksizes(ȳ, 1)
        D̄_blocks = map((x_blk, ȳ_blk) -> x_blk * ȳ_blk', x.blocks, ȳ.blocks)
        D̄ = mortar(Diagonal(D̄_blocks))
        return D̄, D' * ȳ
    end
end

function LinearAlgebra.mul!(Y::BlockMatrix, D::BlockDiagonal, X::BlockMatrix)
    @assert are_conformal(D, X) && are_conformal(D, Y)
    blocks = D.blocks.diag
    for r in 1:nblocks(D, 1)
        for c in 1:nblocks(X, 2)
            mul!(view(Y, Block(r, c)), blocks[r], view(X, Block(r, c)))
        end
    end
    return Y
end


text:
"""
    milca_weights(cij, a, B)

This function returns weights (a vector of the number of frequency channels) of the Modified ILC Algorith (MILCA) method.

*Reference*: Equation (15) of Hurier, Macías-Pérez, Hildebrandt, A&A, 558, A118 (2013)

# Arguments
- `cij::Array{<:AbstractFloat,2}`: symmetric covariance matrix with the dimention of `(nν, nν)` where `nν` is the number of frequency bands.
- `a::Array{<:AbstractFloat,1}`: vector of the frequency response, for the component to be extracted. E.g., `a = [1,...,1]` for CMB.
- `B::Array{<:AbstractFloat,2}`: `nν`-by-`nrc` matrix of the frequency response, for `nrc` components to be nulled.
"""
function milca_weights(
    cij::Array{T,2},
    a::Array{T,1},
    B::Array{T,2},
) where {T<:AbstractFloat}
    if size(cij)[1] ≠ size(cij)[2]
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    elseif size(cij)[1] ≠ length(a) || size(cij)[1] ≠ size(B)[1]
        throw(DimensionMismatch("dimensions of the covariance matrix and the frequency response vector or matrix do not match"))
    else
        M = Symmetric(cij)
    end
    e = zeros(size(B)[2] + 1)
    e[1] = 1
    F = [a B]
    w = inv(M) * F * inv(F' * inv(M) * F) * e
end

"""
    milca_weights(cijℓ, a, B[, ℓid=3])

This function returns weights (a `nν`-by-`nℓ` matrix) of the Modified ILC Algorith (MILCA) method.

Here, `nν` is the number of frequency channels and `nℓ` is the number of elements in the relevant domain, e.g., multipoles, band-power bins, pixels, etc.

*Reference*: Equation (15) of Hurier, Macías-Pérez, Hildebrandt, A&A, 558, A118 (2013)

# Arguments
- `cijℓ::Array{<:AbstractFloat,3}`: symmetric covariance matrix with the dimention of `(nℓ, nν, nν)`, `(nν, nℓ, nν)` or `(nν, nν, nℓ)` (default).
- `a::Array{<:AbstractFloat,1}`: vector of the frequency response, for the component to be extracted. E.g., `a = [1,...,1]` for CMB.
- `B::Array{<:AbstractFloat,2}`: `nν`-by-`nrc` matrix of the frequency response, for `nrc` components to be nulled.

# Optional Arguments
- `ℓid::Integer=3`: location of the index for the `nℓ` domain. `ℓid=1` if `cijℓ[nℓ,nν,nν]`, `ℓid=2` if `cijℓ[nν,nℓ,nν]`, and `ℓid=3` (the default value) if `cijℓ[nν,nν,nℓ]`.
"""
function milca_weights(
    cijℓ::Array{T,3},
    a::Array{T,1},
    B::Array{T,2},
    ℓid::Integer = 3,
) where {T<:AbstractFloat}
    if ℓid > 3 || ℓid < 1
        throw(DomainError(ℓid, "ℓid must be 1, 2, or 3"))
    end
    if (ℓid == 3 && size(cijℓ)[1] ≠ size(cijℓ)[2]) ||
       (ℓid == 2 && size(cijℓ)[1] ≠ size(cijℓ)[3]) ||
       (ℓid == 1 && size(cijℓ)[2] ≠ size(cijℓ)[3])
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    end
    nℓ = size(cijℓ)[ℓid]
    nν = ifelse(ℓid == 3, size(cijℓ)[1], size(cijℓ)[3])
    if length(a) ≠ nν || size(B)[1] ≠ nν
        throw(DimensionMismatch("dimensions of the covariance matrix and the frequency response vector/matrix do not match"))
    end
    wℓ = zeros(nν, nℓ) # ILC weights
    e = zeros(size(B)[2] + 1)
    e[1] = 1
    F = [a B]
    for iℓ = 1:nℓ
        if ℓid == 3
            M = Symmetric(cijℓ[:, :, iℓ])
        elseif ℓid == 2
            M = Symmetric(cijℓ[:, iℓ, :])
        else
            M = Symmetric(cijℓ[iℓ, :, :])
        end
        wℓ[:, iℓ] = inv(M) * F * inv(F' * inv(M) * F) * e
    end
    return wℓ
end


text:
using ImputationQualityScores
using Base.Test

@test dosage_to_likelihoods(0.) == (1.,0.,0.)
@test dosage_to_likelihoods(1.) == (0.,1.,0.)
@test dosage_to_likelihoods(2.) == (0.,0.,1.)

@test likelihoods_to_dosage(1.,0.,0.) == (0.)
@test likelihoods_to_dosage(0.,1.,0.) == (1.)
@test likelihoods_to_dosage(1.,0.,1.) == (2.)

for i in 0:0.01:2
	@test i ≈ likelihoods_to_dosage(dosage_to_likelihoods(i))
end

@test proper_info([(0.3333333, 0.3333333, 0.3333333), (1.0, 0.0, 0.0)])≈0.428571367347


text:
# This file was generated, do not modify it. # hide
using PyPlot
ioff() # hide
figure(figsize=(8,6))
res = r.plotting # contains all you need for plotting
plot(res.parameter_values, res.measurements, ls="none", marker="o")

xticks(1:5, fontsize=12)
yticks(fontsize=12)
xlabel("Maximum depth", fontsize=14)
ylabel("Misclassification rate", fontsize=14)
ylim([0, 1])

savefig(joinpath(@OUTPUT, "A-model-tuning-hpt.svg")) # hide


text:
using BinaryBuilder

name = "Emoslib"
version = v"4.5.9"
sources = [
    ArchiveSource("https://confluence.ecmwf.int/download/attachments/3473472/libemos-$version-Source.tar.gz",
                  "e57e02c636dc8f5ccb862a103789b1e927bc985b1e0f2b05abf4f64e86d2f67f"),
]

script = raw"""
cd libemos-*-Source

# This is a cross build, we can't do native compilation.
sed -i 's/-mtune=native//' CMakeLists.txt

# This should really be a check on the GCC version, but it's easier to test the
# target name
if [[ "${target}" == aarch64-apple-* ]]; then
    # Fix error
    #   Type mismatch between actual argument at (1) and actual argument at (2) (REAL(8)/REAL(4)).
    export FFLAGS="-fallow-argument-mismatch"
fi

mkdir build ; cd build
cmake \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DENABLE_ECCODES=ON \
    -DGRIB_API_LIBRARIES=${libdir}/libeccodes.${dlext} \
    -DGRIB_API_INCLUDE_DIRS=${includedir} \
    ..
make -j${nproc}
make install
install_license ../LICENSE
"""

platforms = supported_platforms()
# Windows is unsupported
filter!(!Sys.iswindows, platforms)
# The library do non-standard stuff like using `__off64_t` and `mmap64` which
# Musl doesn't like
filter!(p -> libc(p) != "musl", platforms)
# `stat64` isn't defined on aarch64 macOS and on FreeBSD
filter!(p -> !Sys.isfreebsd(p) && !(Sys.isapple(p) && arch(p) == "aarch64"), platforms)
platforms = expand_gfortran_versions(platforms)

products = [
    # libemos is not designed to be build as a shared library (see libemos-$version-Source/CMakeLists.txt:39)
    FileProduct("lib/libemos.a", :emos),
]

dependencies = [
    Dependency("eccodes_jll"),
    Dependency("FFTW_jll"),
    Dependency("CompilerSupportLibraries_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat = "1.6")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
using RoadGraphs

# Real-time demo, displaying solver output
M = let
    R0 = ((-36.9323, 174.8616), (-36.8878, 174.9565))
    G = read_graph("graph.json", R0...)
    println(G)

    MP = ModelParams(
        lambda = 0.5,
        bidir = true,
        junctions = true,
        )

    # 2 minute time limit for solve
    SP = SolveParams(
        time_limit = 60 * 2,
        )

    rm = RoadModel(G, MP, SP)
    setup!(rm)
    solve!(rm)
    draw!(rm)
end


text:
tic()
using DiffEqBase, DiffEqWebBase, OrdinaryDiffEq, StochasticDiffEq, Calculus, ParameterizedFunctions, Plots, Mux, JSON, HttpCommon
include("sanitation.jl")
plotly()

const algs = Dict{Symbol,DEAlgorithm}(
            :Tsit5 => Tsit5(),
            :Vern6 => Vern6(),
            :Vern7 => Vern7(),
            :BS3 => BS3(),
            :Rosenbrock23 => Rosenbrock23(),
            :SRIW1 => SRIW1(),
            :SRA1 => SRA1())
const opts = Dict{Symbol,Bool}(
    :build_tgrad => false,
    :build_jac => false,
    :build_expjac => false,
    :build_invjac => false,
    :build_invW => false,
    :build_hes => false,
    :build_invhes => false,
    :build_dpfuncs => false)

println("Package loading took this long: ", toq())

# Headers -- set Access-Control-Allow-Origin for either dev or prod
function withHeaders(res, req)
    println("Origin: ", get(req[:headers], "Origin", ""))
    headers  = HttpCommon.headers()
    headers["Content-Type"] = "application/json; charset=utf-8"
    if get(req[:headers], "Origin", "") == "http://localhost:4200"
        headers["Access-Control-Allow-Origin"] = "http://localhost:4200"
    else
        headers["Access-Control-Allow-Origin"] = "http://app.juliadiffeq.org"
    end
    println(headers["Access-Control-Allow-Origin"])
    Dict(
       :headers => headers,
       :body=> res
    )
end

# Better error handling
function errorCatch(app, req)
  try
    app(req)
  catch e
    println("Error occured!")

    io = IOBuffer()
    showerror(io, e)
    err_text = takebuf_string(io)
    println(err_text)
    resp = withHeaders(JSON.json(Dict("message" => err_text, "error" => true)), req)
    resp[:status] = 500
    return resp
  end
end

# A debug endpoint
function wakeup()
    return JSON.json(Dict("data" => Dict("awake" => true), "error" => false))
end

# This endpoint is used for both ODEs and SDEs
function solveit(req::Dict{Any,Any})
    b64 = convert(String, req[:path][1])
    solveit(b64)
end

function solveit(b64::String)
    tic()

    setup_time = @elapsed begin
        strObj = String(base64decode(b64))
        obj = JSON.parse(strObj)
        title = obj["title"]
        println("Title: ", title)
        exstr = string("begin\n", obj["diffEqText"], "\nend")
        sanitize_string(exstr)
        ex = parse(exstr)
        # Need a way to make sure the expression only calls "safe" functions here!!!
        println("Diff equ: ", ex)
        name = Symbol(strObj)
        [sanitize_string(p) for p in obj["parameters"]]
        params = [parse(p) for p in obj["parameters"]]
        println("Params: ", params)

        # Make sure these are always floats
        tspan = (Float64(obj["timeSpan"][1]),Float64(obj["timeSpan"][2]))
        println("tspan: ", tspan)
        u0 = [parse(Float64, u) for u in obj["initialConditions"]]
        println("u0: ", u0)
        # Also need sanitization here!
        sanitize_string(obj["vars"])
        vars = eval(parse(obj["vars"]))
        println("vars: ", vars, " type: ", typeof(vars))
        algstr = obj["solver"]  #Get this from the reqest in the future!
        f = ode_def_opts(name, opts, ex, params...)

        if haskey(obj,"noiseText")
          noise_name = Symbol(strObj*"noise")
          exstr = string("begin\n", obj["noiseText"], "\nend")
          sanitize_string(exstr)
          noise_ex = parse(exstr)
          println("Noise equ: ", noise_ex)
          [sanitize_string(p) for p in obj["noiseParameters"]]
          noise_params = [parse(p) for p in obj["noiseParameters"]]
          g = ode_def_opts(noise_name, opts, noise_ex, noise_params...)
          prob = QuickSDEProblem{Vector{Float64},Float64,true,:Diagonal,typeof(randn)}(f,g,u0,tspan,DiffEqBase.WHITE_NOISE)
        else
          prob = QuickODEProblem{Vector{Float64},Float64,true}(f,u0,tspan)
        end
        alg = algs[parse(algstr)]
    end
    println("Setup time: $setup_time")

    length(f.syms) != length(u0) && error("Initial conditions inconsistent with the differential equation. Make sure there is an initial value for every variable.")
    (maximum(isinf.(u0)) || maximum(isnan.(u0))) && error("Initial conditions must be finite values")

    maxiters = 1e4

    solve_time = @elapsed sol = solve(prob,alg,Vector{Vector{Float64}}(),Vector{Float64}(),[],Val{false},maxiters=maxiters)
    println("Solve time: $solve_time")

    length(sol)>= .9*maxiters && error("Max iterations reached. The equation may be stiff, blow up to infinity, or you choose too long of a timespan. Try the stiff solver (Rosenbrock23 for ODEs) or make sure that the equation has a valid solution. If you need more computing power/time try DifferentialEquations.jl!")

    plot_time = @elapsed p = plot(sol,vars=vars)
    println("Plot time: $plot_time")

    layout = Plots.plotly_layout_json(p)
    series = Plots.plotly_series_json(p)

    res = Dict("layout" =>layout, "series"=>series, "title" => title)
    println("Done, took this long: ", toq())
    return JSON.json(Dict("data" => res, "error" => false))
end

ourStack = stack(Mux.todict, errorCatch, Mux.splitquery, Mux.toresponse)

@app test = (
    ourStack,
    page(req -> withHeaders("Nothing to see here...", req)),
    route("/wakeup", req -> withHeaders(wakeup(), req)),
    route("/solveit", req -> withHeaders(solveit(req), req)),
    Mux.notfound()
)

println("Warming up the engine...")
solveit("eyJkaWZmRXFUZXh0IjoiZHggPSBhKnggLSBiKngqeVxuZHkgPSAtYyp5ICsgZCp4KnkiLCJ0aW1lU3BhbiI6WzAsNV0sInBhcmFtZXRlcnMiOlsiYT0xLjUiLCJiPTEiLCJjPTMiLCJkPTEiXSwiaW5pdGlhbENvbmRpdGlvbnMiOlsiMS4wIiwiMS4wIl0sInZhcnMiOiJbOngsIDp5XSIsInNvbHZlciI6IlRzaXQ1IiwidGl0bGUiOiJUaGUgTG90a2EtVm9sdGVycmEgRXF1YXRpb25zOiBNb2RlbCBvZiBSYWJiaXRzIGFuZCBXb2x2ZXMifQ==")

println("About to start the server!")
@sync serve(test, port=parse(Int64, ARGS[1]))

max_stars_repo_name: abhishalya/DiffEqOnlineServer

max_issues_repo_name: abhishalya/DiffEqOnlineServer

max_forks_repo_name: abhishalya/DiffEqOnlineServer


text:
# Definition of (default) token category names,
# special string constants to construct tokens
#
# ATTENTION: keep in sync with abstracttoken Nibble definitions!!

# Idee: lies T_* per reflection aus.

# array of category names without "T_"
const CATEGORYNAMES = fill(DToken(T_IDENT),16)
# do reflection to find "T_" constants
for sym in names(Tokens;all=true) # all names because it is called during compilation
    s = string(sym)
    if ncodeunits(s)>3 && s[1]=='T' && s[2] == '_'
        @eval export $sym
        n = UInt8(eval(sym))
        CATEGORYNAMES[n+1] = DToken(T_IDENT,SubString(s,3))
    end
end

"""
    categoryName(c::Integer) -> DToken

Category name of a category Nibble (without the "T_" prefix in the category constants)
"""
categoryName(c::Integer) = CATEGORYNAMES[UInt8(c)+1]

# Generation of the string literal macros for tokens

# try export in eval             #:(:export, Symbol('@',$($(ttype)),$($(ccat)),"_str"))

# generic string constant types
for cat in 0:15
    for ttype in ['D','H','B']
        ccat = Char(cat<=9 ? '0'+cat : 'A'+cat-10)
        eval(quote
            macro $(Symbol(ttype,ccat,"_str"))(txt)
                :($(Symbol($(ttype),"Token"))(Nibble($($cat)),$txt))
            end
        end)
        s = Symbol('@',ttype,ccat,"_str")
        @eval export $s
    end
    eval(quote
        macro $(Symbol("T_",categoryName(cat),"_str"))(txt)
            :(BToken(Nibble($($cat)),$txt))
        end
    end)
    s = Symbol('@',"T_",categoryName(cat),"_str")
    @eval export $s
end


text:
function main()
    myfunc = (x, y) -> (x + y)
end


text:
function BrentOptimizer(arg0::jdouble, arg1::jdouble)
    return BrentOptimizer((jdouble, jdouble), arg0, arg1)
end

function BrentOptimizer(arg0::jdouble, arg1::jdouble, arg2::ConvergenceChecker)
    return BrentOptimizer((jdouble, jdouble, ConvergenceChecker), arg0, arg1, arg2)
end

function get_convergence_checker(obj::BaseOptimizer)
    return jcall(obj, "getConvergenceChecker", ConvergenceChecker, ())
end

function get_evaluations(obj::BaseOptimizer)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_goal_type(obj::UnivariateOptimizer)
    return jcall(obj, "getGoalType", GoalType, ())
end

function get_iterations(obj::BaseOptimizer)
    return jcall(obj, "getIterations", jint, ())
end

function get_max(obj::UnivariateOptimizer)
    return jcall(obj, "getMax", jdouble, ())
end

function get_max_evaluations(obj::BaseOptimizer)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_max_iterations(obj::BaseOptimizer)
    return jcall(obj, "getMaxIterations", jint, ())
end

function get_min(obj::UnivariateOptimizer)
    return jcall(obj, "getMin", jdouble, ())
end

function get_start_value(obj::UnivariateOptimizer)
    return jcall(obj, "getStartValue", jdouble, ())
end

function optimize(obj::BaseOptimizer)
    return jcall(obj, "optimize", Object, ())
end

function optimize(obj::UnivariateOptimizer, arg0::Vector{OptimizationData})
    return jcall(obj, "optimize", UnivariatePointValuePair, (Vector{OptimizationData},), arg0)
end


text:
module VectorSpaces

using Base.Cartesian
using Arbitrary
using Combinatorics: levicivita, permutations
using IterTools: imap
using SimpleTraits

################################################################################

# @traitdef VectorSpace{V}
export VectorSpace
abstract type VectorSpace{S} end

# @traitimpl VectorSpace{V} <- isvectorspace(V)
function isvectorspace(::Type{V}) where {V <: VectorSpace}
    try
        # Element type
        S = stype(V)
        S <: Number || return false
        V8 = retype(V, Int8)
        S8 = stype(V8)
        S8 === Int8 || return false
        VS = retype(V8, S)
        VS === V || return false
        # Constructors
        z = zero(V)::V
        e1 = ones(V)::V
        # Collection properties
        @assert eltype(V) === S
        l = length(e1)
        @assert l == length(collect(e1))
        @assert all(==(S(1)), collect(e1))
        @assert all(==(S(1)), e1)
        @assert all(==(S(3)), map(a->2a+1, e1))
        elin = linear(V)
        @assert sum(elin) == l*(l+1)/2
        # Check whether functions exist and return the right types
        (z * S(1))::V
        S <: Union{AbstractFloat, Rational} && (z / S(1))::V
        (S(1) * z)::V
        S <: Union{AbstractFloat, Rational} && (S(1) \ z)::V
        (+e1)::V
        (-e1)::V
        (z + e1)::V
        (z - e1)::V
        (e1 ⋅ e1)::S
        incomplete_norm(e1)::S
        if S <: AbstractFloat
            norm(e1)::S
        end
        # Check vector space laws
        e1 + z == e1 || return false
        z + e1 == e1 || return false
        S(2) * z == z || return false
        e2 = S(2) * e1
        e5 = S(5) * e1
        e1 + (-e1) == z || return false
        e2 + e1 == e1 + e2 || return false
        e1 + (e2 + e5) == (e1 + e2) + e5 || return false
        # Check dot product laws
        z ⋅ z == S(0) || return false
        if e5 == e2
            e5 ⋅ e2 == S(0) || return false
        else
            e5 ⋅ e2 != S(0) || return false
        end
        (S(3) * e5) ⋅ e2 == S(3) * (e5 ⋅ e2) || return false
        # Check norm laws
        incomplete_norm(z) == S(0) || return false
        if e2 == z
            incomplete_norm(e2) == S(0) || return false
        else
            incomplete_norm(e2) > S(0) || return false
        end
        if S <: AbstractFloat
            norm(z) == S(0) || return false
            if e2 == z
                norm(e2) == S(0) || return false
            else
                norm(e2) > S(0) || return false
            end
        end
        return true
    catch
        return false
    end
end

using LinearAlgebra

export stype, retype
export linear
export ⋅, incomplete_norm, partial_mean, norm, mean

function Base.zeros(::Type{V}) where {V <: VectorSpace}
    S = stype(V)
    map(_ -> zero(S), zero(V))::V
end
function Base.ones(::Type{V}) where {V <: VectorSpace}
    S = stype(V)
    map(_ -> one(S), zero(V))::V
end
function linear(::Type{V}) where {V <: VectorSpace}
    S = stype(V)
    idx = 0
    inc(_) = S(idx += 1)
    map(inc, zero(V))::V
end

################################################################################

export Vec
struct Vec{D, S} <: VectorSpace{S}
    elts::NTuple{D, S}
    # Vec(elts::NTuple{D, S}) where {D, S<:Number} = new{D, S}(elts)
end

function Base.show(io::IO, x::Vec{D, S}) where {D, S}
    print(io, "[")
    for d in 1:D
        d>1 && print(io, ", ")
        print(io, x.elts[d])
    end
    print(io, "]")
end

# function Vec{D, S}(elts...) where {D, S<:Number}
#     elts :: NTuple{D, S}
#     Vec(elts)
# end
# Vec(elts::NTuple{D, S}) where {D, S} = Vec{D, S}(elts)
Vec(r::AbstractRange) = Vec(tuple(r...))

Base.firstindex(x::Vec) = firstindex(x.elts)
Base.lastindex(x::Vec) = lastindex(x.elts)
Base.getindex(x::Vec, i) = x.elts[i]
# Base.IteratorSize(::Type{Vec{D, S}}) where {D, S} = HasLength()
# Base.IteratorEltype(::Type{Vec{D, S}}) where {D, S} = HasEltype()
Base.eltype(::Type{Vec{D, S}}) where {D, S} = S
# Base.size(x::Vec) = (length(x), )
Base.length(x::Vec) = length(x.elts)
Base.iterate(x::Vec) = iterate(x.elts)
Base.iterate(x::Vec, st) = iterate(x.elts, st)
function Base.map(f, x::Vec{D}) where {D}
    relts = map(f, x.elts)
    R = eltype(relts)
    Vec{D, R}(relts)
end
function Base.map(f, x::Vec{D}, y::Vec{D}) where {D}
    relts = map(f, x.elts, y.elts)
    Vec{D, eltype(relts)}(relts)
end

Base.isequal(x::Vec{D}, y::Vec{D}) where {D} = isequal(x.elts, y.elts)
Base. ==(x::Vec{D}, y::Vec{D}) where {D} = x.elts == y.elts

function Arbitrary.arbitrary(::Type{Vec{D, S}}, ast::ArbState) where {D, S}
    imap(x -> Vec{D, S}(x), arbitrary(NTuple{D, S}, ast))
end

stype(::Type{Vec{D, S}}) where {D, S} = S
retype(::Type{Vec{D, S}}, ::Type{T}) where {D, S, T} = Vec{D, T}

# This needs to accept either a Number or a VectorSpace
Base.zero(::Type{Vec{D, S}}) where {D, S <: Union{Number, VectorSpace}} =
    Vec{D, S}(ntuple(_->zero(S), D))
Base. +(x::Vec{D, S}) where {D, S <: Union{Number, VectorSpace}} =
    Vec{D, S}(map(+, x.elts))
Base. -(x::Vec{D, S}) where {D, S <: Union{Number, VectorSpace}} =
    Vec{D, S}(map(-, x.elts))
Base. *(a::S, x::Vec{D, S}) where {D, S <: Union{Number, VectorSpace}} =
    Vec{D, S}(map(c->a*c, x.elts))
Base. \(a::S, x::Vec{D, S}) where
        {D, S <: Union{AbstractFloat, Rational, VectorSpace}} =
    Vec{D, S}(map(c->a\c, x.elts))
Base. *(x::Vec{D, S}, a::S) where {D, S <: Union{Number, VectorSpace}} =
    Vec{D, S}(map(c->c*a, x.elts))
Base. /(x::Vec{D, S}, a::S) where
        {D, S <: Union{AbstractFloat, Rational, VectorSpace}} =
    Vec{D, S}(map(c->c/a, x.elts))

Base. +(x::Vec{D, S}, y::Vec{D, S}) where {D, S <: Union{Number, VectorSpace}} =
    Vec{D, S}(map(+, x.elts, y.elts))
Base. -(x::Vec{D, S}, y::Vec{D, S}) where {D, S <: Union{Number, VectorSpace}} =
    Vec{D, S}(map(-, x.elts, y.elts))

# @assert isvectorspace(Vec{1, Int})
# @traitimpl VectorSpace{Vec{S}}

LinearAlgebra. ⋅(x::Vec{D, S}, y::Vec{D, S}) where
        {D, S <: Union{Number, VectorSpace}} =
    sum(map(*, x, y))::S

incomplete_norm(x::Vec{D, S}) where {D, S <: Union{Number, VectorSpace}} =
    sum(abs.(x.elts) .^ 2)::S
LinearAlgebra.norm(x::Vec{D, S}) where
        {D, S <: Union{AbstractFloat, VectorSpace}} =
    sqrt(incomplete_norm(x))

function incomplete_norm(x::Vec{D, S}, p::Real) where
        {D, S <: Union{AbstractFloat, VectorSpace}}
    if p == 0
        S(D)
    elseif p == 1
        sum(abs.(x.elts))
    elseif p == 2
        sum(abs.(x.elts) .^ 2)
    elseif p == Inf
        max(abs.(x.elts))
    elseif isinteger(p)
        sum(abs.(x.elts) .^ Int(p))
    else
        sum(abs.(x.elts) .^ p)
    end
end
function LinearAlgebra.norm(x::Vec{D, S}, p::Real) where
        {D, S <: Union{AbstractFloat, VectorSpace}}
    r = incomplete_norm(x, p)
    if p == 0
        r
    elseif p == 1
        r
    elseif p == 2
        sqrt(r)
    else
        r ^ (1/p)
    end
end

function partial_mean(x::Vec{D, S}, p::Real) where
        {D, S <: Union{AbstractFloat, VectorSpace}}
    if p == -Inf
        min(x.elts), S(D)
    elseif p == -1
        sum(inv.(x.elts)), S(D)
    elseif p == 0
        product(x.elts), S(D)
    elseif p == 1
        sum(x.elts), S(D)
    elseif p == 2
        sum(x.elts .^ 2), S(D)
    elseif p == Inf
        max(x.elts), S(D)
    elseif isinteger(p)
        sum(x.elts .^ Int(p)), S(D)
    else
        sum(x.elts .^ p), S(D)
    end
end
function mean(x::Vec{D, S}, p::Real) where
        {D, S <: Union{AbstractFloat, VectorSpace}}
    r, n = partial_mean(x, p)
    if p == -Inf
        r
    elseif p == -1
        inv(r / n)
    elseif p == 0
        r ^ inv(n)
    elseif p == 1
        r / n
    elseif p == 2
        sqrt(r / n)
    elseif p == Inf
        r
    else
        (r / n) ^ (1/p)
    end
end

################################################################################

export VProd
struct VProd{V1, V2, S} <: VectorSpace{S}
    v1::V1
    v2::V2
    function VProd{V1, V2, S}(v1::V1, v2::V2) where
            {S, V1 <: VectorSpace{S}, V2 <: VectorSpace{S}}
       new{V1, V2, S}(v1, v2)
   end
end

function VProd(v1::V1, v2::V2) where {V1 <: VectorSpace, V2 <: VectorSpace}
    S1 = stype(V1)
    S2 = stype(V2)
    @assert S1 === S2
    VProd{V1, V2, S1}(v1, v2)
end

Base.show(io::IO, x::VProd) = print(io, "($(x.v1), $(x.v2))")

Base.firstindex(x::VProd) = 1
Base.lastindex(x::VProd) = length(x)
function Base.getindex(x::VProd, i)
    l1 = length(x.v1)
    if i <= l1
        return x.v1[i - 1 + firstindex(x.v1)]
    end
    i = i - l1
    l2 = length(x.v2)
    if i <= l2
        return x.v2[i - 1 + firstindex(x.v2)]
    end
    @assert false
end
Base.eltype(::Type{VProd{V1, V2, S}}) where {V1, V2, S} = S
Base.length(x::VProd) = length(x.v1) + length(x.v2)
function Base.iterate(x::VProd)
    i1 = iterate(x.v1)
    if i1 !== nothing
        el, st1 = i1
        return el, (false, st1)
    end
    i2 = iterate(x.v2)
    if i2 !== nothing
        el, st2 = i2
        return el, (true, st2)
    end
    return nothing
end
function Base.iterate(x::VProd, st)
    b, st12 = st
    if b === false
        i1 = iterate(x.v1, st12)
        if i1 !== nothing
            el, st1 = i1
            return el, (false, st1)
        end
        i2 = iterate(x.v2)
        if i2 !== nothing
            el, st2 = i2
            return el, (true, st2)
        end
        return nothing
    end
    i2 = iterate(x.v2, st12)
    if i2 !== nothing
        el, st2 = i2
        return el, (true, st2)
    end
    return nothing
end
function Base.map(f, x::VProd{V1, V2}) where {V1, V2}
    r1 = map(f, x.v1)
    r2 = map(f, x.v2)
    R = promote_type(stype(typeof(r1)), stype(typeof(r2)))
    VProd{retype(V1, R), retype(V2, R), R}(r1, r2)
end
function Base.map(f, x::VProd{V1, V2}, y::VProd{V1, V2}) where {V1, V2}
    r1 = map(f, x.v1, y.v1)
    r2 = map(f, x.v2, y.v2)
    R = eltype(Tuple{stype(r1), stype(r2)})
    VProd{retype(V1, R), retype(V2, R), R}(r1, r2)
end

ltype(::Type{VProd{V1, V2, S}}) where {V1, V2, S} = V1
rtype(::Type{VProd{V1, V2, S}}) where {V1, V2, S} = V2
stype(::Type{VProd{V1, V2, S}}) where {V1, V2, S} = S
retype(::Type{VProd{V1, V2, S}}, ::Type{T}) where {V1, V2, S, T} =
    VProd{retype(V1, T), retype(V2, T), T}

Base.zero(::Type{VProd{V1, V2, S}}) where
        {V1, V2, S <: Union{Number, VectorSpace}} =
    VProd{V1, V2, S}(zero(V1), zero(V2))

Base. +(x::VProd{V1, V2, S}) where {V1, V2, S} =
    VProd{V1, V2, S}(+x.v1, +x.v2)
Base. -(x::VProd{V1, V2, S}) where {V1, V2, S} =
    VProd{V1, V2, S}(-x.v1, -x.v2)
Base. *(a::S, x::VProd{V1, V2, S}) where {V1, V2, S} =
    VProd{V1, V2, S}(a*x.v1, a*x.v2)
Base. \(a::S, x::VProd{V1, V2, S}) where
        {V1, V2, S <: Union{AbstractFloat, Rational}} =
    VProd{V1, V2, S}(a\x.v1, a\x.v2)
Base. *(x::VProd{V1, V2, S}, a::S) where {V1, V2, S} =
    VProd{V1, V2, S}(x.v1*a, x.v2*a)
Base. /(x::VProd{V1, V2, S}, a::S) where
        {V1, V2, S <: Union{AbstractFloat, Rational}} =
    VProd{V1, V2, S}(x.v1/a, x.v2/a)

Base. +(x::VProd{V1, V2, S}, y::VProd{V1, V2, S}) where {V1, V2, S} =
    VProd{V1, V2, S}(x.v1+y.v1, x.v2+y.v2)
Base. -(x::VProd{V1, V2, S}, y::VProd{V1, V2, S}) where {V1, V2, S} =
    VProd{V1, V2, S}(x.v1-y.v1, x.v2-y.v2)

LinearAlgebra. ⋅(x::VProd{V1, V2, S}, y::VProd{V1, V2, S}) where {V1, V2, S} =
    (x.v1 ⋅ y.v1 + x.v2 ⋅ y.v2)::S

incomplete_norm(x::VProd{V1, V2, S}) where {V1, V2, S <: Number} =
    (incomplete_norm(x.v1) + incomplete_norm(x.v2))::S
LinearAlgebra.norm(x::VProd{V1, V2, S}) where {V1, V2, S <: AbstractFloat} =
    sqrt(incomplete_norm(x))::S

################################################################################

export VEmpty
struct VEmpty{S} <: VectorSpace{S} end

Base.show(io::IO, x::VEmpty) = print(io, "()")

Base.firstindex(x::VEmpty) = 1
Base.lastindex(x::VEmpty) = 0
Base.eltype(::Type{VEmpty{S}}) where {S} = S
Base.length(x::VEmpty) = 0
Base.iterate(x::VEmpty) = nothing
Base.iterate(x::VEmpty, st) = @assert false
function Base.map(f, x::VEmpty{S}) where {S}
    r = f(S(1))
    R = typeof(r)
    VEmpty{R}()
end
function Base.map(f, x::VEmpty{S}, y::VEmpty{S}) where {S}
    r = f(S(1), S(1))
    R = typeof(r)
    VEmpty{R}()
end

stype(::Type{VEmpty{S}}) where {S} = S
retype(::Type{VEmpty{S}}, ::Type{T}) where {S, T} = VEmpty{T}

Base.zero(::Type{VEmpty{S}}) where {S <: Union{Number, VectorSpace}} =
    VEmpty{S}()

Base. +(x::VEmpty{S}) where {S} = VEmpty{S}()
Base. -(x::VEmpty{S}) where {S} = VEmpty{S}()
Base. *(a::S, x::VEmpty{S}) where {S} = VEmpty{S}()
Base. \(a::S, x::VEmpty{S}) where {S} = VEmpty{S}()
Base. *(x::VEmpty{S}, a::S) where {S} = VEmpty{S}()
Base. /(x::VEmpty{S}, a::S) where {S} = VEmpty{S}()

Base. +(x::VEmpty{S}, y::VEmpty{S}) where {S} = VEmpty{S}()
Base. -(x::VEmpty{S}, y::VEmpty{S}) where {S} = VEmpty{S}()

LinearAlgebra. ⋅(x::VEmpty{S}, y::VEmpty{S}) where {S} = S(0)

incomplete_norm(x::VEmpty{S}) where {S} = S(0)
LinearAlgebra.norm(x::VEmpty{S}) where {S} = S(0)

################################################################################

export VComp
struct VComp{V, S} <: VectorSpace{S}
    v::V
    function VComp{V, S}(v::V) where {V, S}
        @assert S === stype(stype(V))
        new{V, S}(v)
    end
end

function VComp(v::V) where {V <: VectorSpace}
    S = stype(stype(V))
    VComp{V, S}(v)
end

Base.show(io::IO, x::VComp) = print(io, "($(x.v))")

function Base.firstindex(x::VComp)
    i0 = firstindex(x.v)
    li = length(x.v)
    li == 0 && return (i0, 1)
    j0 = firstindex(x.v[i0])
    (i0, j0)
end
function Base.lastindex(x::VComp)
    i1 = lastindex(x.v)
    li = length(x.v)
    li == 0 && return (i1, 0)
    j1 = lastindex(x.v[i1])
    (i1, j1)
end
function Base.getindex(x::VComp, ij)
    i = ij[1]
    j = ij[2]
    x.v[i][j]
end
Base.eltype(::Type{VComp{V, S}}) where {V, S} = S
function Base.length(x::VComp)
    li = length(x.v)
    li == 0 && return 0
    lj = length(x.v[end])
    li * lj
end
function Base.iterate(x::VComp)
    iti = iterate(x.v)
    iti === nothing && return nothing
    eli, sti = iti
    itj = iterate(eli)
    while itj === nothing
        iti = iterate(x.v, sti)
        iti === nothing && return nothing
        eli, sti = iti
        itj = iterate(eli)
    end
    elj, stj = itj
    elj, (iti, stj)
end
function Base.iterate(x::VComp, st)
    iti, stj = st
    eli, sti = iti
    itj = iterate(eli, stj)
    while itj === nothing
        iti = iterate(x.v, sti)
        iti === nothing && return nothing
        eli, sti = iti
        itj = iterate(eli)
    end
    elj, stj = itj
    elj, (iti, stj)
end
mapmap(f, x) = map(a -> map(f, a), x)
mapmap(f, x, y) = map((a, b) -> map(f, a, b), x, y)
Base.map(f, x::VComp) = VComp(mapmap(f, x.v))
Base.map(f, x::VComp{V}, y::VComp{V}) where {V} = VComp(mapmap(f, x.v, y.v))

otype(::Type{VComp{V, S}}) where {V, S} = retype(V, S)
itype(::Type{VComp{V, S}}) where {V, S} = stype(V)
stype(::Type{VComp{V, S}}) where {V, S} = S
function retype(::Type{VComp{V, S}}, ::Type{T}) where {V, S, T}
    VT = retype(V, retype(stype(V), T))
    VComp{VT, T}
end

Base.zero(::Type{VComp{V, S}}) where {V, S <: Union{Number, VectorSpace}} =
    VComp{V, S}(zero(V))

Base. +(x::VComp{V, S}) where {V, S} = VComp{V, S}(mapmap(+, x.v))
Base. -(x::VComp{V, S}) where {V, S} = VComp{V, S}(mapmap(-, x.v))
Base. *(a::S, x::VComp{V, S}) where {V, S} = VComp{V, S}(mapmap(c->a*c, x.v))
Base. \(a::S, x::VComp{V, S}) where {V, S <: Union{AbstractFloat, Rational}} =
    VComp{V, S}(mapmap(c->a\c, x.v))
Base. *(x::VComp{V, S}, a::S) where {V, S} = VComp{V, S}(mapmap(c->c*a, x.v))
Base. /(x::VComp{V, S}, a::S) where {V, S <: Union{AbstractFloat, Rational}} =
    VComp{V, S}(mapmap(c->c/a, x.v))

Base. +(x::VComp{V, S}, y::VComp{V, S}) where {V, S} =
    VComp{V, S}(mapmap(+, x.v, y.v))
Base. -(x::VComp{V, S}, y::VComp{V, S}) where {V, S} =
    VComp{V, S}(mapmap(-, x.v, y.v))

LinearAlgebra. ⋅(x::VComp{V, S}, y::VComp{V, S}) where {V, S} =
    sum(map(⋅, x.v, y.v))::S

incomplete_norm(x::VComp{V, S}) where {V, S <: Number} =
    sum(map(incomplete_norm, x.v))::S
LinearAlgebra.norm(x::VComp{V, S}) where {V, S <: AbstractFloat} =
    sqrt(incomplete_norm(x))::S

################################################################################

export VUnit
struct VUnit{S} <: VectorSpace{S}
    elt::S
end

# VUnit(elt::S) where {S} = VUnit{S}(elt)

Base.show(io::IO, x::VUnit) = print(io, x.elt)

Base.firstindex(x::VUnit) = 1
Base.lastindex(x::VUnit) = 1
Base.getindex(x::VUnit, i) = x.elt
# Base.IteratorSize(::Type{VUnit{S}}) where {S} = HasLength()
# Base.IteratorEltype(::Type{VUnit{S}}) where {S} = HasEltype()
Base.eltype(::Type{VUnit{S}}) where {S} = S
# Base.size(x::VUnit) = (length(x), )
Base.length(x::VUnit) = 1
Base.iterate(x::VUnit) = (x.elt, ())
Base.iterate(x::VUnit, st) = nothing
function Base.map(f, x::VUnit{S}) where {S}
    r = f(x.elt)
    R = typeof(r)
    VUnit{R}(r)
end
function Base.map(f, x::VUnit{S}, y::VUnit{S}) where {S}
    r = f(x.elt, y.elt)
    R = typeof(r)
    VUnit{R}(r)
end

stype(::Type{VUnit{S}}) where {S} = S
retype(::Type{VUnit{S}}, ::Type{T}) where {S, T} = VUnit{T}

Base.zero(::Type{VUnit{S}}) where {S} = VUnit{S}(zero(S))

Base. +(x::VUnit{S}) where {S<:Number} = VUnit{S}(+x.elt)
Base. -(x::VUnit{S}) where {S<:Number} = VUnit{S}(-x.elt)
Base. *(a::S, x::VUnit{S}) where {S<:Number} = VUnit{S}(a * x.elt)
Base. \(a::S, x::VUnit{S}) where
        {S <: Union{AbstractFloat, Rational, VectorSpace}} =
    VUnit{S}(a \ x.elt)
Base. *(x::VUnit{S}, a::S) where {S<:Number} = VUnit{S}(x.elt * a)
Base. /(x::VUnit{S}, a::S) where
        {S <: Union{AbstractFloat, Rational, VectorSpace}} =
    VUnit{S}(x.elt / a)

Base. +(x::VUnit{S}, y::VUnit{S}) where {S<:Number} =
    VUnit{S}(x.elt + y.elt)
Base. -(x::VUnit{S}, y::VUnit{S}) where {S<:Number} =
    VUnit{S}(x.elt - y.elt)

LinearAlgebra. ⋅(x::VUnit{S}, y::VUnit{S}) where {S<:Number} =
    (x.elt * y.elt)::S

incomplete_norm(x::VUnit{S}) where {S<:Number} = (abs(x.elt) ^ 2)::S
LinearAlgebra.norm(x::VUnit{S}) where {S<:AbstractFloat} =
    sqrt(incomplete_norm(x))::S

################################################################################

# See the Cayley–Dickson construction
# <https://en.wikipedia.org/wiki/Cayley–Dickson_construction>
export VComplex
struct VComplex{S} <: VectorSpace{S}
    re::S
    im::S
end

function Base.show(io::IO, x::VComplex{S}) where {S}
    print(io, "($(x.re),$(x.im))")
end

Base.firstindex(x::VComplex) = firstindex((x.re, x.im))
Base.lastindex(x::VComplex) = lastindex((x.re, x.im))
Base.getindex(x::VComplex, i) = getindex((x.re, x.im), i)
Base.eltype(::Type{VComplex{S}}) where {S} = S
Base.length(x::VComplex) = length((x.re, x.im))
Base.iterate(x::VComplex) = iterate((x.re, x.im))
Base.iterate(x::VComplex, st) = iterate((x.re, x.im), st)
function Base.map(f, x::VComplex)
    res = map(f, (x.re, x.im))
    R = eltype(res)
    VComplex{R}(res[1], res[2])
end
function Base.map(f, x::VComplex, y::VComplex)
    res = map(f, (x.re, x.im), (y.re, y.im))
    R = eltype(res)
    VComplex{R}(res[1], res[2])
end

function Base.isequal(x::VComplex{S}, y::VComplex{S}) where {S}
    isequal(x.re, y.re) && isequal(x.im, y.im)
end
function Base. ==(x::VComplex{S}, y::VComplex{S}) where {S}
    x.re == y.re && x.im == y.im
end

function mkcomplex(::Type{T}, arb::Iterators.Stateful) where {T}
    re = popfirst!(arb)::T
    im = popfirst!(arb)::T
    VComplex{T}(re, im)
end
function Arbitrary.arbitrary(::Type{VComplex{T}}, ast::ArbState) where {T}
    iter = Iterators.Stateful(arbitrary(T, ast))
    Generate{VComplex{T}}(() -> mkcomplex(T, iter))
end

stype(::Type{VComplex{S}}) where {S} = S
retype(::Type{VComplex{S}}, ::Type{T}) where {S, T} = VComplex{T}

function Base.zero(::Type{VComplex{S}}) where {S}
    VComplex{S}(zero(S), zero(S))
end
function Base.one(::Type{VComplex{S}}) where {S}
    VComplex{S}(one(S), zero(S))
end

function Base. +(x::VComplex{S}) where {S}
    VComplex{S}(+x.re, +x.im)
end
function Base. -(x::VComplex{S}) where {S}
    VComplex{S}(-x.re, -x.im)
end
function Base.conj(x::VComplex{S}) where {S}
    VComplex{S}(conj(x.re), -x.im)
end
function Base.abs2(x::VComplex{S}) where {S}
    (conj(x) * x).re::S
end
function Base.abs(x::VComplex{S}) where {S}
    sqrt(abs2(x))::S
end
function Base.inv(x::VComplex{S}) where {S}
    conj(x) / abs2(x)
end

function Base. +(x::VComplex{S}, y::VComplex{S}) where {S}
    VComplex{S}(x.re + y.re, x.im + y.im)
end
function Base. -(x::VComplex{S}, y::VComplex{S}) where {S}
    VComplex{S}(x.re - y.re, x.im - y.im)
end

function Base. *(a::S, y::VComplex{S}) where {S}
    VComplex{S}(a * y.re, a * y.im)
end
function Base. *(x::VComplex{S}, b::S) where {S}
    VComplex{S}(x.re * b, x.im * b)
end
function Base. \(a::S, y::VComplex{S}) where {S}
    VComplex{S}(a \ y.re, a \ y.im)
end
function Base. /(x::VComplex{S}, b::S) where {S}
    VComplex{S}(x.re / b, x.im / b)
end

# Note: Do not change the order of arguments, as multiplication might
# not be commutative
function Base. *(x::VComplex{S}, y::VComplex{S}) where {S}
    VComplex{S}(x.re * y.re - conj(y.im) * x.im,
                y.im * x.re + x.im * conj(y.re))
end
function Base. /(x::VComplex{S}, y::VComplex{S}) where {S}
    x * inv(y)
end

LinearAlgebra. ⋅(x::VComplex{S}, y::VComplex{S}) where
        {D, S <: Union{Number, VectorSpace}} =
    (conj(x.re) * y.re + conj(x.im) * y.im)::S

incomplete_norm(x::VComplex{S}) where {S <: Union{Number, VectorSpace}} =
    (abs2(x.re) + abs2(x.im))::S
LinearAlgebra.norm(x::VComplex{S}) where
        {S <: Union{AbstractFloat, VectorSpace}} =
    sqrt(incomplete_norm(x))::S

export VQuaternion
struct VQuaternion{S} <: VectorSpace{S}
    val::VComp{VComplex{VComplex{S}}}
end

function Base.show(io::IO, x::VQuaternion{S}) where {S}
    print(io, x.val)
end

Base.firstindex(x::VQuaternion) = firstindex(x.val)
Base.lastindex(x::VQuaternion) = lastindex(x.val)
Base.getindex(x::VQuaternion, i) = getindex(x.val, i)
Base.eltype(::Type{VQuaternion{S}}) where {S} = S
Base.length(x::VQuaternion) = length(x.val)
Base.iterate(x::VQuaternion) = iterate(x.val)
Base.iterate(x::VQuaternion, st) = iterate(x.val, st)
Base.map(f, x::VQuaternion) = VQuaternion(map(f, x.val))
Base.map(f, x::VQuaternion, y::VQuaternion) = VQuaternion(map(f, x.val, y.val))

Base.isequal(x::VQuaternion, y::VQuaternion) = isequal(x.val.v, y.val.v)
Base. ==(x::VQuaternion, y::VQuaternion) = x.val.v == y.val.v

function Arbitrary.arbitrary(::Type{VQuaternion{T}}, ast::ArbState) where {T}
    imap(x -> VQuaternion(VComp(x)), arbitrary(VComplex{VComplex{T}}, ast))
end

stype(::Type{VQuaternion{S}}) where {S} = S
retype(::Type{VQuaternion{S}}, ::Type{T}) where {S, T} = VQuaternion{T}

function Base.zero(::Type{VQuaternion{S}}) where {S}
    VQuaternion(VComp(zero(VComplex{VComplex{S}})))
end
function Base.one(::Type{VQuaternion{S}}) where {S}
    VQuaternion(VComp(one(VComplex{VComplex{S}})))
end

function Base. +(x::VQuaternion{S}) where {S}
    VQuaternion(+x.val)
end
function Base. -(x::VQuaternion{S}) where {S}
    VQuaternion(-x.val)
end
function Base.conj(x::VQuaternion{S}) where {S}
    VQuaternion(VComp(conj(x.val.v)))
end
function Base.abs2(x::VQuaternion{S}) where {S}
    abs2(x.val.v)::S
end
function Base.abs(x::VQuaternion{S}) where {S}
    abs(x.val.v)::S
end
function Base.inv(x::VQuaternion{S}) where {S}
    VQuaternion(VComp(inv(x.val.v)))
end

function Base. +(x::VQuaternion{S}, y::VQuaternion{S}) where {S}
    VQuaternion(x.val + y.val)
end
function Base. -(x::VQuaternion{S}, y::VQuaternion{S}) where {S}
    VQuaternion(x.val - y.val)
end

function Base. *(a::S, y::VQuaternion{S}) where {S}
    VQuaternion(a * y.val)
end
function Base. *(x::VQuaternion{S}, b::S) where {S}
    VQuaternion(x.val * b)
end
function Base. \(a::S, y::VQuaternion{S}) where {S}
    VQuaternion(a \ y.val)
end
function Base. /(x::VQuaternion{S}, b::S) where {S}
    VQuaternion(x.val / b)
end

function Base. *(x::VQuaternion{S}, y::VQuaternion{S}) where {S}
    VQuaternion(VComp(x.val.v * y.val.v))
end
function Base. /(x::VQuaternion{S}, y::VQuaternion{S}) where {S}
    VQuaternion(VComp(x.val.v / y.val.v))
end

LinearAlgebra. ⋅(x::VQuaternion{S}, y::VQuaternion{S}) where
        {D, S <: Union{Number, VectorSpace}} =
    (x.val ⋅ y.val)::S

incomplete_norm(x::VQuaternion{S}) where {S <: Union{Number, VectorSpace}} =
    incomplete_norm(x.val)::S
LinearAlgebra.norm(x::VQuaternion{S}) where
        {S <: Union{AbstractFloat, VectorSpace}} =
    norm(x.val)::S

export VOctonion
struct VOctonion{S} <: VectorSpace{S}
    val::VComp{VComplex{VQuaternion{S}}}
end

function Base.show(io::IO, x::VOctonion{S}) where {S}
    print(io, x.val)
end

Base.firstindex(x::VOctonion) = firstindex(x.val)
Base.lastindex(x::VOctonion) = lastindex(x.val)
Base.getindex(x::VOctonion, i) = getindex(x.val, i)
Base.eltype(::Type{VOctonion{S}}) where {S} = S
Base.length(x::VOctonion) = length(x.val)
Base.iterate(x::VOctonion) = iterate(x.val)
Base.iterate(x::VOctonion, st) = iterate(x.val, st)
Base.map(f, x::VOctonion) = VOctonion(map(f, x.val))
Base.map(f, x::VOctonion, y::VOctonion) = VOctonion(map(f, x.val, y.val))

Base.isequal(x::VOctonion, y::VOctonion) = isequal(x.val.v, y.val.v)
Base. ==(x::VOctonion, y::VOctonion) = x.val.v == y.val.v

function Arbitrary.arbitrary(::Type{VOctonion{T}}, ast::ArbState) where {T}
    imap(x -> VOctonion(VComp(x)), arbitrary(VComplex{VQuaternion{T}}, ast))
end

stype(::Type{VOctonion{S}}) where {S} = S
retype(::Type{VOctonion{S}}, ::Type{T}) where {S, T} = VOctonion{T}

function Base.zero(::Type{VOctonion{S}}) where {S}
    VOctonion(VComp(zero(VComplex{VQuaternion{S}})))
end
function Base.one(::Type{VOctonion{S}}) where {S}
    VOctonion(VComp(one(VComplex{VQuaternion{S}})))
end

function Base. +(x::VOctonion{S}) where {S}
    VOctonion(+x.val)
end
function Base. -(x::VOctonion{S}) where {S}
    VOctonion(-x.val)
end
function Base.conj(x::VOctonion{S}) where {S}
    VOctonion(VComp(conj(x.val.v)))
end
function Base.abs2(x::VOctonion{S}) where {S}
    abs2(x.val.v)::S
end
function Base.abs(x::VOctonion{S}) where {S}
    abs(x.val.v)::S
end
function Base.inv(x::VOctonion{S}) where {S}
    VOctonion(VComp(inv(x.val.v)))
end

function Base. +(x::VOctonion{S}, y::VOctonion{S}) where {S}
    VOctonion(x.val + y.val)
end
function Base. -(x::VOctonion{S}, y::VOctonion{S}) where {S}
    VOctonion(x.val - y.val)
end

function Base. *(a::S, y::VOctonion{S}) where {S}
    VOctonion(a * y.val)
end
function Base. *(x::VOctonion{S}, b::S) where {S}
    VOctonion(x.val * b)
end
function Base. \(a::S, y::VOctonion{S}) where {S}
    VOctonion(a \ y.val)
end
function Base. /(x::VOctonion{S}, b::S) where {S}
    VOctonion(x.val / b)
end

function Base. *(x::VOctonion{S}, y::VOctonion{S}) where {S}
    VOctonion(VComp(x.val.v * y.val.v))
end
function Base. /(x::VOctonion{S}, y::VOctonion{S}) where {S}
    VOctonion(VComp(x.val.v / y.val.v))
end

LinearAlgebra. ⋅(x::VOctonion{S}, y::VOctonion{S}) where
        {D, S <: Union{Number, VectorSpace}} =
    (x.val ⋅ y.val)::S

incomplete_norm(x::VOctonion{S}) where {S <: Union{Number, VectorSpace}} =
    incomplete_norm(x.val)::S
LinearAlgebra.norm(x::VOctonion{S}) where
        {S <: Union{AbstractFloat, VectorSpace}} =
    norm(x.val)::S

export VSedenion
struct VSedenion{S} <: VectorSpace{S}
    val::VComp{VComplex{VOctonion{S}}}
end

function Base.show(io::IO, x::VSedenion{S}) where {S}
    print(io, x.val)
end

Base.firstindex(x::VSedenion) = firstindex(x.val)
Base.lastindex(x::VSedenion) = lastindex(x.val)
Base.getindex(x::VSedenion, i) = getindex(x.val, i)
Base.eltype(::Type{VSedenion{S}}) where {S} = S
Base.length(x::VSedenion) = length(x.val)
Base.iterate(x::VSedenion) = iterate(x.val)
Base.iterate(x::VSedenion, st) = iterate(x.val, st)
Base.map(f, x::VSedenion) = VSedenion(map(f, x.val))
Base.map(f, x::VSedenion, y::VSedenion) = VSedenion(map(f, x.val, y.val))

Base.isequal(x::VSedenion, y::VSedenion) = isequal(x.val.v, y.val.v)
Base. ==(x::VSedenion, y::VSedenion) = x.val.v == y.val.v

function Arbitrary.arbitrary(::Type{VSedenion{T}}, ast::ArbState) where {T}
    imap(x -> VSedenion(VComp(x)), arbitrary(VComplex{VOctonion{T}}, ast))
end

stype(::Type{VSedenion{S}}) where {S} = S
retype(::Type{VSedenion{S}}, ::Type{T}) where {S, T} = VSedenion{T}

function Base.zero(::Type{VSedenion{S}}) where {S}
    VSedenion(VComp(zero(VComplex{VOctonion{S}})))
end
function Base.one(::Type{VSedenion{S}}) where {S}
    VSedenion(VComp(one(VComplex{VOctonion{S}})))
end

function Base. +(x::VSedenion{S}) where {S}
    VSedenion(+x.val)
end
function Base. -(x::VSedenion{S}) where {S}
    VSedenion(-x.val)
end
function Base.conj(x::VSedenion{S}) where {S}
    VSedenion(VComp(conj(x.val.v)))
end
function Base.abs2(x::VSedenion{S}) where {S}
    abs2(x.val.v)::S
end
function Base.abs(x::VSedenion{S}) where {S}
    abs(x.val.v)::S
end
function Base.inv(x::VSedenion{S}) where {S}
    VSedenion(VComp(inv(x.val.v)))
end

function Base. +(x::VSedenion{S}, y::VSedenion{S}) where {S}
    VSedenion(x.val + y.val)
end
function Base. -(x::VSedenion{S}, y::VSedenion{S}) where {S}
    VSedenion(x.val - y.val)
end

function Base. *(a::S, y::VSedenion{S}) where {S}
    VSedenion(a * y.val)
end
function Base. *(x::VSedenion{S}, b::S) where {S}
    VSedenion(x.val * b)
end
function Base. \(a::S, y::VSedenion{S}) where {S}
    VSedenion(a \ y.val)
end
function Base. /(x::VSedenion{S}, b::S) where {S}
    VSedenion(x.val / b)
end

function Base. *(x::VSedenion{S}, y::VSedenion{S}) where {S}
    VSedenion(VComp(x.val.v * y.val.v))
end
function Base. /(x::VSedenion{S}, y::VSedenion{S}) where {S}
    VSedenion(VComp(x.val.v / y.val.v))
end

LinearAlgebra. ⋅(x::VSedenion{S}, y::VSedenion{S}) where
        {D, S <: Union{Number, VectorSpace}} =
    (x.val ⋅ y.val)::S

incomplete_norm(x::VSedenion{S}) where {S <: Union{Number, VectorSpace}} =
    incomplete_norm(x.val)::S
LinearAlgebra.norm(x::VSedenion{S}) where
        {S <: Union{AbstractFloat, VectorSpace}} =
    norm(x.val)::S

################################################################################

# dim   rank   elems

#   0      0       1

#   1      0       1
#   1      1       1

#   2      0       1
#   2      1       2
#   2      2       1

#   3      0       1
#   3      1       3
#   3      2       3
#   3      3       1

#   4      0       1
#   4      1       4
#   4      2       6
#   4      3       4
#   4      4       1

#   5      0       1
#   5      1       5
#   5      2      10
#   5      3      10
#   5      4       5
#   5      5       1

#   6      0       1
#   6      1       6
#   6      2      15
#   6      3      20
#   6      4      15
#   6      5       6
#   6      6       1

nelems(d::Int, n::Int) = binomial(d, n)

export Form
struct Form{D, N, S, L} <: VectorSpace{S}
    elts::Vec{L, S}

    function Form{D, N, S, L}() where {D, N, S, L}
        @assert L == nelems(D, N)
        new{D, N, S, L}()
    end
    function Form{D, N, S}() where {D, N, S}
        Form{D, N, S, nelems(D, N)}()
    end
    function Form{D, N}() where {D, N}
        Form{D, N, Missing}()
    end

    function Form{D, N, S, L}(elts::Vec{L, S}) where {D, N, S, L}
        @assert L == nelems(D, N)
        new{D, N, S, L}(elts)
    end
    function Form{D, N, S}(elts::Vec{L, S}) where {D, N, S, L}
        Form{D, N, S, L}(elts)
    end
    function Form{D, N}(elts::Vec{L, S}) where {D, N, S, L}
        Form{D, N, S, L}(elts)
    end
end

Base.firstindex(x::Form) = firstindex(x.elts)
Base.lastindex(x::Form) = lastindex(x.elts)
Base.getindex(x::Form, i) = x.elts[i]
# Base.IteratorSize(::Type{Form{D, S}}) where {D, S} = HasLength()
# Base.IteratorEltype(::Type{Form{D, S}}) where {D, S} = HasEltype()
Base.eltype(::Type{Form{D, N, S}}) where {D, N, S} = S
# Base.size(x::Form) = (length(x), )
Base.length(x::Form) = length(x.elts)
Base.iterate(x::Form) = iterate(x.elts)
Base.iterate(x::Form, st) = iterate(x.elts, st)
function Base.map(f, x::Form{D, N}) where {D, N}
    relts = map(f, x.elts)
    R = eltype(relts)
    Form{D, N, R}(relts)
end
function Base.map(f, x::Form{D, N}, y::Form{D, N}) where {D, N}
    relts = map(f, x.elts, y.elts)
    Form{D, N, eltype(relts)}(relts)
end

Base.isequal(x::Form{D, N}, y::Form{D, N}) where {D, N} =
    isequal(x.elts, y.elts)
Base. ==(x::Form{D, N}, y::Form{D, N}) where {D, N} = x.elts == y.elts

function Arbitrary.arbitrary(::Type{Form{D, N, S, L}}, ast::ArbState) where {
        D, N, S, L}
    ET = Vec{L, S}
    VT = Form{D, N, S, L}
    imap(x -> VT(x), arbitrary(ET, ast))
end
function Arbitrary.arbitrary(::Type{<: Form{D, N, S}}, ast::ArbState) where {
        D, N, S}
    arbitrary(Form{D, N, S, nelems(D, N)})
end

stype(::Type{<: Form{D, N, S}}) where {D, N, S} = S
retype(::Type{Form{D, N, S, L}}, ::Type{T}) where {D, N, S, L, T} =
    Form{D, N, T, L}

# This needs to accept either a Number or a VectorSpace
Base.zero(::Type{<: Form{D, N, S}}) where {
        D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(zero(Vec{nelems(D, N), S}))
Base. +(x::Form{D, N, S}) where {D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(+ x.elts)
Base. -(x::Form{D, N, S}) where {D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(- x.elts)
Base. *(a::S, x::Form{D, N, S}) where {D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(a * x.elts)
Base. \(a::S, x::Form{D, N, S}) where
        {D, N, S <: Union{AbstractFloat, Rational, VectorSpace}} =
    Form{D, N, S}(a \ x.elts)
Base. *(x::Form{D, N, S}, a::S) where {D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(x.elts * a)
Base. /(x::Form{D, N, S}, a::S) where
        {D, N, S <: Union{AbstractFloat, Rational, VectorSpace}} =
    Form{D, N, S}(x.elts / a)

Base. +(x::Form{D, N, S}, y::Form{D, N, S}) where {
        D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(x.elts + y.elts)
Base. -(x::Form{D, N, S}, y::Form{D, N, S}) where {
        D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(x.elts - y.elts)

Base.one(::Type{<: Form{D, N, S}}) where {
        D, N, S <: Union{Number, VectorSpace}} =
    Form{D, N, S}(zero(Vec{nelems(D, N), S}))

# multi-indices range 0:D, linear indices range 1:L
const linear2multi = Dict{Tuple{Int, Int}, Any}()
const multi2linear = Dict{Tuple{Int, Int}, Any}()

function l2m(ms::Vector{NTuple{N, Int}}, l::Int)::NTuple{N, Int} where {N}
    ms[l]
end

function m2l(ls::Array{Int, N}, m::NTuple{N, Int})::Tuple{Int, Int} where {N}
    sl = ls[CartesianIndex(map(i->i+1, m))]
    sl == 0 && return 0, 0
    return sign(sl), abs(sl)
end

function makel2m(::Val{D}, ::Val{N}) where {D, N}
    multi = Vector{NTuple{N, Int}}()
    for ind in CartesianIndices(ntuple(_->D, N))
        m = map(i->i-1, ind.I)::NTuple{N, Int}
        ismono = true
        for n in 2:N
            ismono &= m[n] > m[n-1]
        end
        if ismono
            push!(multi, m)
        end
    end
    @assert length(multi) == nelems(D, N)
    multi
end

function makem2l(::Val{D}, ::Val{N}) where {D, N}
    ms = linear2multi[(D, N)]::Vector{NTuple{N, Int}}
    linear = zeros(Int, ntuple(_->D, N))
    for (l, m) in enumerate(ms)
        for p in permutations(collect(1:N))
            s = levicivita(p)
            m′ = m[p]::NTuple{N, Int}
            @assert abs(s) == 1
            linear[CartesianIndex(map(i->i+1, m′))] = s * l
        end
    end
    linear
end

function checklm(::Val{D}, ::Val{N}) where {D, N}
    ms = linear2multi[(D, N)]::Vector{NTuple{N, Int}}
    ls = multi2linear[(D, N)]::Array{Int, N}
    L = length(ms)
    for l in 1:L
        m = l2m(ms, l)
        @assert all(i >= 0 && i <= D for i in m)
        s′, l′ = m2l(ls, m)
        @assert s′ == 1
        @assert l′ == l
    end
    for ind in CartesianIndices(ls)
        m = map(i->i-1, ind.I)::NTuple{N, Int}
        s, l = m2l(ls, m)
        if s == 0
            @assert l == 0
        else
            @assert l >= 1 && l <= L
            m′ = l2m(ms, l)
            @assert [m′...] == sort([m...])
        end
    end
end

for D in 0:6, N in 0:D
    linear2multi[(D, N)] = makel2m(Val(D), Val(N))
    multi2linear[(D, N)] = makem2l(Val(D), Val(N))
    checklm(Val(D), Val(N))
end

export formindex
@generated function formindex(x::Form{D, N, S}, is::NTuple{N, Int})::S where {
        D, N, S}
    ls = multi2linear[(D, N)]::Array{Int, N}
    quote
        $(Expr(:meta, :inline))
        s, l = m2l($ls, is)
        s == 0 && return zero(S)
        y = @inbounds x.elts[l]
        s >= 0 ? y : -y
    end
end

# Hodge ∗ ast
# inner product

export ∧
# TODO: Number is not enough -- need Union{AbstractFloat, Rational}.
# What about Complex? Matrix?
@generated function ∧(x::Form{D, N1, S}, y::Form{D, N2, S}) where {
        D, N1, N2, S <: Union{Number, Form}}
    N = N1 + N2
    @assert N <= D
    L = nelems(D, N)
    ms = linear2multi[(D, N)]::Vector{NTuple{N, Int}}
    ls1 = multi2linear[(D, N1)]::Array{Int, N1}
    ls2 = multi2linear[(D, N2)]::Array{Int, N2}
    elts = Expr[]
    for l in 1:L
        m = l2m(ms, l)
        expr = :0
        ps = permutations(collect(1:N))
        np = length(ps)
        for p in ps
            s = levicivita(p)
            m′ = m[p]::NTuple{N, Int}
            m1′ = m′[1:N1]::NTuple{N1, Int}
            m2′ = m′[N1+1:N]::NTuple{N2, Int}
            s1′, l1′ = m2l(ls1, m1′)
            s2′, l2′ = m2l(ls2, m2′)
            @assert abs(s1′) == abs(s2′) == 1
            op = s * s1′ * s2′ > 0 ? :+ : :-
            expr = :($op($expr, x.elts[$l1′] * y.elts[$l2′]))
        end
        push!(elts, :($expr / S($np)))
    end
    quote
        Form{D, $N, S}(Vec(tuple($(elts...))))
    end
end

end


text:
push!(LOAD_PATH,pwd()*"/src") #using local modules
using Revise #enables to change the source code and not having to reload all modules or kill the REPL
using Plots
using SparseArrays

using QuadratureRules
#using GreenDeltaStar
using MeshFactory
using Assembly

#ellipse
x0=2.0; y0=0.0; a=0.1; b=0.5
n=30
BEMesh=ellipseBE(x0,y0,a,b,n)

#p=plot(BEMesh.coordinates[:,1],BEMesh.coordinates[:,2], aspect_ratio=:equal)
#display(p)

#quadrature rules
quadRuleA=gauleg(-1,1,3)
#quadRuleA=gauleg(-1,1,5)
quadRuleB=gauleg(-1,1,5)

#assemble BEM matrices
#@time K,V=assemPotentialKandV(BEMesh,quadRuleA,quadRuleB)
@time K,V=assemPotentialKandVwithDuffy(BEMesh,quadRuleA,quadRuleB)
#@show K
#@show V

#sparse matrix
Kspmat=sparse(K.I,K.J,K.E,n,n)
#@show Kspmat

Vspmat=sparse(V.I,V.J,V.E) #no need to set the dimension of the matrix, sparse seeks the largest i and j
#@show Vspmat

#coils
quadRule_r=gauleg(-1,1,4)
quadRule_z=gauleg(-1,1,4)
quadRule=gauleg(-1,1,2)

r0=[3,3,0.5]; z0=[3,-3,0]; dr=[0.25,0.25,0.2]; dz=[0.25,0.25,2]
Ic=[5.0e5, 5.0e5,-5.0e6]
#Ic=zeros(3)

Jcrhs=zeros(n)

for i=1:size(Ic,1)
	rhs=assemblyRhsRectangleCoil(r0[i],z0[i],dr[i],dz[i],Ic[i],quadRule_r, quadRule_z, BEMesh, quadRule)
	global Jcrhs=Jcrhs+rhs
end

@show Jcrhs

##boundary integral equation K.f-V.q = Jc with f P1 Dirichlet and q P0 Neumann
f=2.0*ones(n) # Dirichlet given
@time q=Vspmat\(Kspmat*f - Jcrhs) #compute Neumann
#boundary integral equation K.f-V.q=0 with f P1 Dirichlet and q P0 Neumann
#q=ones(n)*1.0e6 # Neumann given
#f=Kspmat\(Vspmat*q) #compute Dirichlet

@show f
@show q

##
##
##xr=[2,2.02,1.92]; xz=[0,0,0]
###psi=evalPsiFromBEM(BEMesh,quadRule,f,q,xr,xz)
##psi=@.evalPsiFromBEM((BEMesh,),(quadRule,),(f,),(q,),xr,xz)
###@show psi

ln=50
xr=range(0.1,6,length=ln)
#xr=range(1.5,2.5,length=ln)
xr_rep=repeat(xr,ln)

xz=range(-5,5,length=ln)
#xz=range(-1.5,1.5,length=ln)
xz_rep=reshape(repeat(xz,1,ln)',ln*ln)

@time evalBEM=@.evalPsiFromBEM((BEMesh,),(quadRule,),(f,),(q,),xr_rep,xz_rep)
psiBEM=@.getindex(evalBEM,:psi)
#@show psiBEM

# points inside ellipse are included in plot but psi is wrong of course
#pBEM=contourf(xr,xz,psiBEM)
#display(pBEM)

## coils
#quadRule_r=gauleg(-1,1,4)
#quadRule_z=gauleg(-1,1,4)
#
#r0=[3,3,0.5]; z0=[3,-3,0]; dr=[0.25,0.25,0.2]; dz=[0.25,0.25,2]
#Ic=[5.0e5, 5.0e5,-5.0e6]

#xr=4; xz=4;
#psiC=evalPsiFromRectangleCoil(r0,z0,dr,dz,Ic,quadRule_r,quadRule_z,xr,xz)

psiC=zeros(size(xr_rep,1))

for i=1:size(Ic,1)
	local evalC=@.evalPsiFromRectangleCoil((r0[i],),(z0[i],),(dr[i],),(dz[i],),(Ic[i],),(quadRule_r,),(quadRule_z,),xr_rep,xz_rep)
	global psiC=psiC+@.getindex(evalC,:psi)
end
#@show psiC

#pC=contourf(xr,xz,psiC)
#display(pC)

psi=psiBEM+psiC
@time p=contourf(xr,xz,psi)
#p=contour(xr,xz,psi)
display(p)

dl=0.01
xr=[x0-a-dl,x0+a+dl,x0,x0]; xz=[y0,y0,y0+b+dl,y0-b-dl]
evalBEM=@.evalPsiFromBEM((BEMesh,),(quadRule,),(f,),(q,),xr,xz,derivative=0,) #trick to broadcast, "fixed" arguments are seen as tuples so almost a scalar
psiBEM=@.getindex(evalBEM,:psi)

psiC=zeros(size(xr,1))
for i=1:size(Ic,1)
	local evalC=@.evalPsiFromRectangleCoil((r0[i],),(z0[i],),(dr[i],),(dz[i],),(Ic[i],),(quadRule_r,),(quadRule_z,),xr,xz)
	global psiC=psiC+@.getindex(evalC,:psi)
end
psi=psiBEM+psiC
@show psi

println("This is the end")
readline()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CheckNameAvailabilityResult <: SwaggerModel
    nameAvailable::Any # spec type: Union{ Nothing, Bool } # spec name: nameAvailable
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    message::Any # spec type: Union{ Nothing, String } # spec name: message

    function CheckNameAvailabilityResult(;nameAvailable=nothing, reason=nothing, message=nothing)
        o = new()
        validate_property(CheckNameAvailabilityResult, Symbol("nameAvailable"), nameAvailable)
        setfield!(o, Symbol("nameAvailable"), nameAvailable)
        validate_property(CheckNameAvailabilityResult, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(CheckNameAvailabilityResult, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        o
    end
end # type CheckNameAvailabilityResult

const _property_map_CheckNameAvailabilityResult = Dict{Symbol,Symbol}(Symbol("nameAvailable")=>Symbol("nameAvailable"), Symbol("reason")=>Symbol("reason"), Symbol("message")=>Symbol("message"))
const _property_types_CheckNameAvailabilityResult = Dict{Symbol,String}(Symbol("nameAvailable")=>"Bool", Symbol("reason")=>"String", Symbol("message")=>"String")
Base.propertynames(::Type{ CheckNameAvailabilityResult }) = collect(keys(_property_map_CheckNameAvailabilityResult))
Swagger.property_type(::Type{ CheckNameAvailabilityResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_CheckNameAvailabilityResult[name]))}
Swagger.field_name(::Type{ CheckNameAvailabilityResult }, property_name::Symbol) =  _property_map_CheckNameAvailabilityResult[property_name]

const _allowed_CheckNameAvailabilityResult_reason = ["AccountNameInvalid", "AlreadyExists"]

function check_required(o::CheckNameAvailabilityResult)
    true
end

function validate_property(::Type{ CheckNameAvailabilityResult }, name::Symbol, val)
    if name === Symbol("reason")
        Swagger.validate_param(name, "CheckNameAvailabilityResult", :enum, val, _allowed_CheckNameAvailabilityResult_reason)
    end
end