text:
# B-spline space

abstract type AbstractBSplineSpace end

@doc raw"""
Construct B-spline space from given polynominal degree and knot vector.
```math
\mathcal{P}[p,k]
```
This type `BSplineSpace` is slower than `FastBSplineSpace`, but this type is not limited with degree.
"""
struct BSplineSpace <: AbstractBSplineSpace
    degree::Int
    knots::Knots
    function BSplineSpace(degree::Integer, knots::Knots)
        if degree < 0
            throw(DomainError(p, "degree of polynominal must be non-negative"))
        end
        new(degree, knots)
    end
end

"""
convert AbstractBSplineSpace to BSplineSpace
"""
function BSplineSpace(P::AbstractBSplineSpace)
    return BSplineSpace(degree(P), knots(P))
end

function degree(P::BSplineSpace)
    return P.degree
end

function knots(P::BSplineSpace)
    return P.knots
end

function bsplineunity(P::AbstractBSplineSpace)
    p = degree(P)
    k = knots(P)
    return k[1+p]..k[end-p]
end

@doc raw"""
Return dimention of a B-spline space.
```math
\dim(\mathcal{P}[p,k])
=\sharp k - p -1
```
"""
function dim(bsplinespace::AbstractBSplineSpace)
    p = degree(bsplinespace)
    k = knots(bsplinespace)
    return length(k) - p - 1
end

@doc raw"""
Check inclusive relationship between B-spline spaces.
```math
\mathcal{P}[p,k]
\subseteq\mathcal{P}[p',k']
```
"""
function Base.:⊆(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
    p = degree(P)
    k = knots(P)
    p′ = degree(P′)
    k′ = knots(P′)
    p₊ = p′ - p

    return p₊ ≥ 0 && (k + p₊ * unique(k) ⊆ k′)
end

@doc raw"""
Check inclusive relationship between B-spline spaces.
```math
\mathcal{P}[p,k]
\sqsubseteq\mathcal{P}[p',k']
\Leftrightarrow
\mathcal{P}[p,k]|_{[k_{p+1},k_{l-p}]}
\subseteq\mathcal{P}[p',k']|_{[\sharp k'_{p'+1},k'_{\sharp k'-p'}]}
```
"""
function issqsubset(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
    p = degree(P)
    k = knots(P)
    p′ = degree(P′)
    k′ = knots(P′)
    p₊ = p′ - p

    if p₊ < 0
        return false
    elseif bsplineunity(P) ≠ bsplineunity(P′)
        return false
    end

    inner_knots = k[p+2:end-p-1]
    inner_knots′ = k′[p′+2:end-p′-1]

    return inner_knots + p₊ * unique(inner_knots) ⊆ inner_knots′
end

const ⊑ = issqsubset
⊒(l, r) = r ⊑ l
⋢(l, r) = !⊑(l, r)
⋣(l, r) = r ⋢ l

≃(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace) = (P1 ⊑ P2)&(P2 ⊑ P1)

function iszeros(P::AbstractBSplineSpace)
    p = degree(P)
    k = knots(P)
    n = dim(P)
    return [k[i] == k[i+p+1] for i in 1:n]
end

function isproper(P::AbstractBSplineSpace)
    return !|(iszeros(P)...)
end

function properdim(P::AbstractBSplineSpace)
    return dim(P) - sum(iszeros(P))
end

function _knotindex(P,t)
    p = degree(P)
    k = knots(P)
    l = length(k)
    return _knotindex(view(k.vector, 1+p:l-p), t) + p
end


text:
using StrFs, Test

randstr(pieces::AbstractVector, K) = string((rand(pieces) for _ in 1:K)...)
randstr(pieces::String, K) = randstr(collect(pieces), K)

@testset "literal and show" begin
    @test strf"λ the ultimate" ≡ StrF("λ the ultimate")
    str = strf"λ the ultimate"
    @test repr(str) == "\"λ the ultimate\""
end

@testset "conversions" begin
    for _ in 1:100
        str = randstr("aα∃", rand(3:5))
        len = sizeof(str)
        for _ in 1:100
            # to string and back
            S = len + rand(0:2)
            strf = StrF{S}(str)
            @test sizeof(strf) == len
            @test length(strf) == length(str)
            str2 = String(strf)
            @test str2 == str

            # to other strfs types
            ssame = StrF{S}(strf)
            slong = StrF{S + 1}(strf)
            @test ssame ≡ strf
            @test slong.bytes[1:(len+1)] == vcat(codeunits(str), [0x00])
            @test_throws InexactError StrF{len-1}(strf)

            # io
            io = IOBuffer()
            write(io, strf)
            seekstart(io)
            @test strf == read(io, StrF{len})
        end
    end
end

@testset "promotion" begin
    p1 = [strf"a", strf"bb", "ccc"]
    @test p1 isa Vector{String}
    p2 = [strf"a", strf"bb"]
    @test p2 isa Vector{StrF{2}}
end

"Convert to StrF, padding size with `Δs`."
function StrFmultilen(str, Δs)
    S = sizeof(str)
    Any[StrF{S + Δ}(str) for Δ in Δs]
end

@testset "concatenation, comparisons, and hashing" begin
    for _ in 1:1000
        str = randstr("abcηβπ", rand(2:8))
        stra = str * "a"
        strλ = str * "λ"
        fstr = StrFmultilen(str, 0:2)
        fstra = StrFmultilen(stra, 0:2)
        fstrλ = StrFmultilen(strλ, 0:2)
        @test all(str .== fstr)
        @test all(fstr .== permutedims(fstr))
        @test all(fstr .== permutedims(fstr))
        @test all(str .< fstra)
        @test all(fstr .< permutedims(fstra))
        @test all(fstra .< permutedims(fstrλ))
        @test all(hash(str) .== hash.(fstr))
    end
end

@testset "shortening conversions" begin
    str = StrF{6}("foo")
    str3 = StrF{3}(str)
    str4 = StrF{4}(str)
    str5 = StrF{5}(str)
    @test str == str3 == str4 == str5
    @test_throws InexactError StrF{2}(str)
end

@testset "zero length strings" begin
    s0 = StrF{0}("")
    s = StrF("")
    @test length(s) == length(s0) == 0
    @test s == s0 == StrF{1}("")
end

@testset "iteration" begin
    for s in ["abc", "ηβπ", "∇κ≠0 ↔ ∃ζ>0"]
        @test collect(StrF(s)) == collect(s)
    end
    for _ in 1:1000
        s = randstr("aα∃bc", rand(3:5))
        @test collect(StrF(s)) == collect(s)
    end
    S = StrF{9}
    @test eltype(S) ≡ eltype(String)
    @test Base.IteratorSize(S) ≡ Base.IteratorSize(String)
    @test Base.IteratorEltype(S) ≡ Base.IteratorEltype(S)
end

ext: jl


text:
# Realtime Database

using HTTP, JSON

BASE_URL = nothing

"""
In the examples on this page,
you would replace [PROJECT_ID] with the identifier of your Firebase project.
"""

"""
realdb_init(base_url; query = Dict())

Initialize the realtimedb with baseurl to make things easier

# Example

```julia
realdb_init("https://[PROJECT_ID].asia-southeast1.firebasedatabase.app")
```
"""
function realdb_init(base_url; query = Dict())
    global BASE_URL = base_url
    println("BASE_URL set:", BASE_URL)
end

"""
realdb_get(url; query = Dict())

GET request to an endpoint.

# Example
```julia
realdb_init("https://[PROJECT_ID].asia-southeast1.firebasedatabase.app")
realdb_get("/users/jack/name")
```

"""
function realdb_get(url; query = Dict())
    pagesize = 300
    pagetoken = ""
    final_url = "$BASE_URL$url.json"
    println("FINAL URL:", final_url)
    query = Dict{String,Any}("pageSize" => pagesize, "pageToken" => pagetoken)
    res = HTTP.get(final_url; query = query)
    if res.status == 200
        println("GET successful")
    else
        println("GET errored")
    end
    JSON.parse(String(res.body))
end

"""
realdb_post(url, body = "{"name": "real_db_test"}"; query = Dict())

POST request to an endpoint.

# Example
```julia
realdb_init("https://[PROJECT_ID].asia-southeast1.firebasedatabase.app")
# realdb_post("/message_list")
body =Dict("user_id" => "jack", "text" => "Ahoy!")
realdb_post("/message_list",body)
```

## Notes:

According to the REST API documentation for Realtime Database,
a POST is the equivalent of a "push" operation when using the client SDKs.
This push operation always involves adding data under a random ID.
There is no avoiding that for a push.

If you know the name of the node where you want to add data,
you should use a PUT instead. This is the equivalent of using "set"
operation with the client SDKs.

"""
function realdb_post(url, body = Dict("name" => "real_db_test"); query = Dict())
    pagesize = 300
    pagetoken = ""
    final_url = "$BASE_URL$url.json"
    println("FINAL URL:", final_url)
    body = JSON.json(body)
    println("Body:", body)
    query = Dict{String,Any}("pageSize" => pagesize, "pageToken" => pagetoken)
    res = HTTP.post(final_url, "", body; query = query)
    if res.status == 200
        println("POST successful")
    else
        println("POST errored")
    end
    JSON.parse(String(res.body))
end

"""
realdb_patch(url, body = Dict("name"=> "real_db_test"); query = Dict())

`PATCH` request to an endpoint.

# Example
```julia
realdb_init("https://[PROJECT_ID].asia-southeast1.firebasedatabase.app")
body =Dict("last"=>"Jones")
realdb_patch("/users/jack/name/",body)
```
"""
function realdb_patch(url, body = Dict("name" => "real_db_test"); query = Dict())
    pagesize = 300
    pagetoken = ""
    final_url = "$BASE_URL$url.json"
    println("FINAL URL:", final_url)
    query = Dict{String,Any}("pageSize" => pagesize, "pageToken" => pagetoken)
    body = JSON.json(body)
    println("Body:", body)
    res = HTTP.patch(final_url, "", body; query = query)
    if res.status == 200
        println("PATCH successful")
    else
        println("PATCH errored")
    end
    JSON.parse(String(res.body))
end

"""
realdb_delete(url, body = Dict("name"=> "real_db_test"); query = Dict())

`DELETE` request to an endpoint

# Example
```julia
realdb_init("https://[PROJECT_ID].asia-southeast1.firebasedatabase.app")
realdb_delete("/users/jack/name/last")
```
"""
function realdb_delete(url, body = Dict("name" => "real_db_test"); query = Dict())
    pagesize = 300
    pagetoken = ""
    final_url = "$BASE_URL$url.json"
    println("FINAL URL:", final_url)
    query = Dict{String,Any}("pageSize" => pagesize, "pageToken" => pagetoken)
    body = JSON.json(body)
    println("Body:", body)
    res = HTTP.delete(final_url, "", body; query = query)
    if res.status == 200
        println("DELETE successful")
    else
        println("DELETE errored")
    end
    JSON.parse(String(res.body))
end

"""
realdb_put(url, body = Dict("name"=> "real_db_test"); query = Dict())

`PUT` request to a endpoint.

# Example

```julia
realdb_init("https://[PROJECT_ID].asia-southeast1.firebasedatabase.app")
body = Dict("first"=>"Ash", "last"=>"Sparrow")
realdb_put("/users/jack/name",body)
```
"""
function realdb_put(url, body = Dict("name" => "real_db_test"); query = Dict())
    pagesize = 300
    pagetoken = ""
    final_url = "$BASE_URL$url.json"
    println("FINAL URL:", final_url)
    query = Dict{String,Any}("pageSize" => pagesize, "pageToken" => pagetoken)
    body = JSON.json(body)
    println("Body:", body)
    res = HTTP.put(final_url, "", body; query = query)
    if res.status == 200
        println("PUT successful")
    else
        println("PUT errored")
    end
    JSON.parse(String(res.body))
end

"""
readdb_download(url, filename = "test"; query = Dict())

Download request

# Example

```julia

```
"""
function readdb_download(url, filename = "test"; query = Dict())
    pagesize = 300
    pagetoken = ""
    final_url = "$BASE_URL$url.json?download=$filename.txt"
    println("FINAL URL:", final_url)
    query = Dict{String,Any}("pageSize" => pagesize, "pageToken" => pagetoken)
    res = HTTP.get(final_url; query = query)
    if res.status == 200
        println("GET successful")
    else
        println("GET errored")
    end
    JSON.parse(String(res.body))
end


text:
using ControlAllocationSwitching
using Test

@testset "ControlAllocationSwitching.jl" begin
    # Write your tests here.
end


text:
println("Test metropolis")

@testset "metropolis" begin
  @testset "gaussian_random_walk_metropolis" begin
    function sample_with_metropolis(n::Integer, rwm::MCMC.Metropolis)
      x = 0.0
      xs = Float64[]
      for _ in 1:n
        x = MCMC.update(rwm, x, normlogpdf)
        append!(xs, x)
      end
      return xs
    end

    Random.seed!(0)
    xs = sample_with_metropolis(1_000_000, MCMC.StaticRWM(Normal()))
    @test isapprox(mean(xs), 0, atol=1e-2)
    @test isapprox(std(xs), 1, atol=1e-2)
    @test isapprox(quantile(xs, 0.025), -1.96, atol=1e-2)
    @test isapprox(quantile(xs, 0.975), 1.96, atol=1e-2)
  end

  @testset "RWM univariate" begin
    mu1, sd1 = (0, 1)
    mu2, sd2 = (5, .5)
    d1 = Normal(mu1, sd1)
    d2 = Normal(mu2, sd2)

    spl = Gibbs(M3(),
                RWM(:mu, (m, s, x) -> normlogpdf(mu1, sd1, x), Normal(0, 1)),
                RWM(:sigma, (m, s, x) -> normlogpdf(mu2, sd2, x), Normal(0, 1)),
                Conditional(:eta, (m, s) -> s.eta .- 1))

    nsamps = 10000
    Random.seed!(0)
    chain = mcmc(spl, nsamps, discard=500, thin=2).chain
    mu = getindex.(chain, :mu)
    sigma = getindex.(chain, :sigma)

    for (param, d) in zip([mu, sigma], [d1, d2])
      @test isapprox(mean(param), mean(d), atol=0.05)
      @test isapprox(std(param), std(d), atol=0.05)
      @test isapprox(quantile(param, 0.025), quantile(d, .025), atol=0.05)
      @test isapprox(quantile(param, 0.975), quantile(d, .975), atol=0.05)
    end
  end

  @testset "RWM multivariate" begin
    mu1, sd1 = (0, 1)
    mu2, sd2 = (5, .5)
    d1 = Normal(mu1, sd1)
    d2 = Normal(mu2, sd2)
    d3 = MvNormal(0:2, 1)

    spl = Gibbs(M3(),
                RWM(:mu, (m, s, x) -> normlogpdf(mu1, sd1, x), Normal(0, 1)),
                RWM(:sigma, (m, s, x) -> normlogpdf(mu2, sd2, x), Normal(0, 1)),
                RWM(:eta, (m, s, x) -> logpdf(d3, x), MvNormal(3, 0.1/3)))

    nsamps = 50000
    Random.seed!(1)
    chain = mcmc(spl, nsamps, discard=10000, thin=2).chain
    mu = getindex.(chain, :mu)
    sigma = getindex.(chain, :sigma)

    for (param, d) in zip([mu, sigma], [d1, d2])
      @test isapprox(mean(param), mean(d), atol=0.05)
      @test isapprox(std(param), std(d), atol=0.05)
      @test isapprox(quantile(param, 0.025), quantile(d, .025), atol=0.05)
      @test isapprox(quantile(param, 0.975), quantile(d, .975), atol=0.05)
    end

    eta = getindex.(chain, :eta)
    @test isapprox(mean(eta), mean(d3), atol=0.2)
    @test isapprox(std(eta), LinearAlgebra.diag(d3.Σ), atol=0.2)
  end

  @testset "RWM bijector" begin
    d = Beta(2, 8)

    spl = Gibbs(M4(),
                RWM(:theta, (m, s, x) -> logpdf(d, x), Normal(0, 1),
                    bijector=Bijectors.Logit(0, 1)),
                RWM(:mu, (m, s, x) -> logpdf(Normal(), x), Normal(0, 1)))

    nsamps = 50000
    Random.seed!(1)
    chain = mcmc(spl, nsamps, discard=10000, thin=2).chain
    theta = getindex.(chain, :theta)

    for (param, d) in zip([theta], [d])
      @test isapprox(mean(param), mean(d), atol=1e-2)
      @test isapprox(std(param), std(d), atol=1e-2)
      @test isapprox(quantile(param, 0.025), quantile(d, .025), atol=1e-2)
      @test isapprox(quantile(param, 0.975), quantile(d, .975), atol=1e-2)
    end
  end

  @testset "MvAdaptiveRWM" begin
    Random.seed!(1)
    K = 10
    d1 = let
      M = randn(K, K)
      MvNormal(randn(K), M*M')
    end
    d2 = Normal()

    spl = Gibbs(M5(K),
                RWM(:theta, (m, s, x) -> logpdf(d1, x), mvarwm(randn(K))),
                RWM(:mu, (m, s, x) -> logpdf(d2, x), Normal(0, 1)))

    nsamps = 50000
    chain = mcmc(spl, nsamps, discard=10000, thin=2).chain
    theta = getindex.(chain, :theta)
    mu = getindex.(chain, :mu)

    @test isapprox(mean(theta), mean(d1), atol=0.2)
    @test isapprox(std(theta), sqrt.(LinearAlgebra.diag(cov(d1))), atol=0.2)
    @test isapprox(mean(mu), mean(d2), atol=1e-2)
    @test isapprox(std(mu), std(d2), atol=1e-2)
  end
end


text:
#=--------------------------------------------------------------------

geometry - Functions for some basic geometric operations

Part of the ImageProjectiveGeometry Module

Copyright (c) 2016 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK June 2016

---------------------------------------------------------------------=#
using LinearAlgebra

export ray2raydist, circleintersectionpts
export rectintersect, pointinconvexpoly, isconvexpoly
export convexpolyintersect

#--------------------------------------------------------------------------
"""
ray2raydist -  Minimum distance between two 3D rays

```
Usage: d = ray2raydist(p1, v1, p2, v2)

Arguments:
      p1, p2 - 3D points that lie on rays 1 and 2.
      v1, v2 - 3D vectors defining the direction of each ray.

Returns:
           d - The minimum distance between the rays.
```

Each ray is defined by a point on the ray and a vector giving the direction
of the ray.  Thus a point on ray 1 will be given by  p1 + alpha*v1  where
alpha is some scalar.
"""
function ray2raydist(p1::Vector, v1::Vector, p2::Vector, v2::Vector)

    # Get vector perpendicular to both rays
    n = cross(v1, v2)
    nnorm = sqrt(dot(n,n))

    # Check if lines are parallel. If so, form a vector perpendicular to v1
    # that is within the plane formed by the parallel rays.
    if nnorm < eps()
        n = cross(v1, p1 .- p2) # Vector perpendicular to plane formed by pair
                                # of rays.
        n = cross(v1, n)        # Vector perpendicular to v1 within the plane
                                # formed by the 2 rays.
        nnorm = sqrt(dot(n,n))

        if nnorm < eps()      # p1 and p2 must be coincident
            return  0.0
        end
    end

    n /= nnorm         # Unit vector

    # The vector joining the two closest points on the rays is:
    #    d*n = p2 + beta*v2 - (p1 + alpha*v1)
    # for some unknown values of alpha and beta.
    #
    # Take dot product of n on both sides
    #    d*n.n = p2.n + beta*v2.n - p1.n - alpha*v1.n
    #
    # as n is prependicular to v1 and v2 this reduces to
    #   d = (p2 - p1).n
    return abs(dot((p2 .- p1), n))
end

#--------------------------------------------------------------------------
"""
circleintersectionpts - Finds intersections of two circles.

Function computes the intersection points between two circles given
their centres and radii.

```
Usage: (i1, i2) = circleintersectionpts(c1, r1, c2, r2, lr)
             i1 = circleintersectionpts(c1, r1, c2, r2, lr)

Arguments:
        c1, c2 - 2-vectors specifying the centres of the two circles.
        r1, r2 - Radii of the two circles.
            lr - Optional string specifying what solutions are wanted:
                "l" for the solution to the left of the line from c1 to c2
                "r" for the solution to the right of the line from c1 to c2
                "lr" if both solutions are wanted (default).

Returns:
     (i1, i2) - Tuple of 2D intersection points if both solutions requested.
          i1  - Single 2D intersection point if a single solution requested.
```

If no solution exists empty vectors are returned for i1, i2.  If the
two circles are identical there are an infinite number of solutions.
In this case we simply choose to subtract/add r1 in the x direction to
c1 to obtain the 'left' and 'right' solutions.
"""
function circleintersectionpts(c1i::Vector, r1::Real, c2i::Vector, r2::Real, lr="lr")

    c1 = float(c1i) # Ensure c1 and c2 are floats for type stability
    c2 = float(c2i)

    maxmag = maximum(abs.([c1; c2; r1; r2]))  # Maximum value in data input
    tol = (maxmag+1)*eps()              # Tolerance used for deciding whether values are equal

    bv = c2 .- c1        # Baseline vector from c1 to c2
    b = norm(bv)         # The distance between the centres

    # Trap case of baseline of zero length.  If r1 == r2 we have a
    # valid geometric situation, but there are an infinite number of
    # solutions.  Here we simply choose to subtract/add r1 in the x
    # direction to c1 to obtain the 'left' and 'right' solutions.
    if b < tol && abs(r1-r2) < tol
        i1 = c1 .- [r1, 0]
        i2 = c1 .+ [r1, 0]
    else
        bv = bv/b             # Normalised baseline vector.
        bvp = [-bv[2], bv[1]] # Vector perpendicular to baseline

        # Trap the degenerate cases where one of the radii are zero, or nearly zero
        if r1 < tol && abs(b - r2) < tol
            i1 = c1
            i2 = c1
        elseif r2 < tol && abs(b - r1) < tol
            i1 = c2
            i2 = c2

        # Check triangle inequality
        elseif b > (r1+r2) || r1 > (b+r2) || r2 > (b+r1)
            i1 = zeros(eltype(c1), 0)
            i2 = zeros(eltype(c1), 0)

        else   # Normal solution
            cosR2 = (b^2 + r1^2 - r2^2)/(2*b*r1)
            sinR2 = sqrt(1-cosR2^2)

            i1 = c1 .+ r1*cosR2*bv .+ r1*sinR2*bvp   # 'left' solution
            i2 = c1 .+ r1*cosR2*bv .- r1*sinR2*bvp   # and 'right solution
        end
    end

    # Return requested solution(s)
    if lr == "lr"
        return (i1, i2)

    elseif lr == "l"
        return i1

    elseif lr == "r"
        return i2

    else
        error("illegal left/right solution request")
    end

end

#--------------------------------------------------------------------------
"""
rectintersect - Test if rectangles intersect

```
Usage:  intersect = rectintersect(r1, r2)

Arguments:  r1, r2 - The two rectangles to be tested where the rectangles
                     are defined using 4-vectors with the elements
                     defining  [left, bottom, width, height]

Returns:             Boolean result
```

See also: convexpolyintersect()
"""
function rectintersect(r1::Vector, r2::Vector)

    intersect = (r1[1] <= r2[1]+r2[3] &&
                 r2[1] <= r1[1]+r1[3] &&
                 r1[2] <= r2[2]+r2[4] &&
                 r2[2] <= r1[2]+r1[4])
end

#--------------------------------------------------------------------------
"""
pointinconvexpoly - Determine if a 2D point is within a convex polygon

```
Usage:  v = pointinconvexpoly(p, poly)

Arguments:   p - 2-vector specifying the point.
          poly - Convex polygon defined as a series of vertices in
                 sequence, clockwise or anticlockwise around the polygon as
                 a 2 x N array.

Returns:    v - +1 if within the polygon
                -1 if outside
                 0 if on the boundary
```

Note:  There is no check to see if the polygon is indeed convex, use
isconvexpoly() if needed

See also: isconvexpoly()
"""
function pointinconvexpoly(p::Vector, poly::Array)

    (dim, N) = size(poly)

    if length(p) != 2 || dim != 2
        error("Data must be 2D")
    end

    # Strategy: Determine whether, in traveling from p to a vertex and
    # then to the following vertex, we turn clockwise or anticlockwise.
    c = zeros(Int, N)
    for n = 1:N-1
        c[n] = clockwise(p, view(poly,:,n), view(poly,:,n+1))
    end
    c[N] = clockwise(p, view(poly,:,N), view(poly,:,1))

    # If for every vertex we turn consistently clockwise, or
    # consistently anticlockwise we are inside the polygon.  If for
    # one of the vertices we did not turn clockwise or anticlockwise
    # then we must be on the boundary.
    if all(c .>= 0) || all(c .<= 0)
        if any(c .== 0)  # We are on the boundary
            v = 0
        else             # We are inside
            v = 1
        end
    else                 # We are outside
        v = -1
    end

    return v
end

#----------------------------------------------------------------------
# Determine whether, in traveling from p1 to p2 to p3 we turn clockwise or
# anticlockwise.  Returns +1 for clockwise, -1 for anticlockwise, and 0 for
# p1, p2, p3 in a straight line.

function clockwise(p1::AbstractArray, p2::AbstractArray, p3::AbstractArray)

    # Form vectors p1->p2 and p2->p3 with z component = 0, form cross product
    # if the resulting z value is -ve the vectors turn clockwise, if +ve
    # anticlockwise, and if 0 the points are in a line.

    return -sign((p2[1]-p1[1])*(p3[2]-p2[2]) - (p2[2]-p1[2])*(p3[1]-p2[1]))

end

#----------------------------------------------------------------------
"""
convexpolyintersect - Returns true if convex polygons intersect.

```
Usage:  intersect = convexpolyintersect(p1, p2)

Arguments:  p1, p2 - The two convex polygons to be tested where the
                     polygons are defined as 2xN arrays of vertex
                     coordinates in sequence around the polygon.

Returns:             Boolean result.
```

Note: There is no check to see if the polygons are indeed convex, use
isconvexpoly() if needed

See also: isconvexpoly(), rectintersect()
"""
function convexpolyintersect(p1::Array, p2::Array)

    # Function to test if two 1D segments overlap.
    # Each segment is defined by a 2 element vector with 1st element
    # being the minimum x value and 2nd element the maximum x value.
    function  overlap(s1::Vector, s2::Vector)
        return  s2[2] >= s1[1] && s1[2] >= s2[1]
    end

    if size(p1,1) != 2 || size(p2,1) != 2
        error("Polygon data must be in 2xN arrays")
    end

    # Build a sequence of vertices that trace the edges of both
    # polygons.  We prepend the last vertex to each polygon list to
    # form a complete cycle.  We also concatenate the two polygon
    # arrays together, this adds an extra 'edge' that connects the
    # last vertex of the 1st polygon with the last vertex of the 2nd
    # polygon.  This allows us to run a single loop of tests below.
    p = [p1[:,end] p1  p2[:,end] p2]

    for n = 2:size(p,2)
        # Generate a 2x2 rotation matrix describing a frame with the y axis
        # aligned with the polygon side.
        y = p[:,n] - p[:,n-1]
        normy = sqrt(y[1]^2 + y[2]^2)
        if normy < eps()  # repeated vertex
            continue
        end
        y = y./normy
        x = [y[2], -y[1]]
        R = [x y]'

        # Rotate all vertices to frame defined by this side
        p1r = R*p1
        p2r = R*p2

        # Check whether x coordinates of rotated p1 vertices overlap
        # with x coordinates of rotated p2 vertices.  If they do not
        # overlap we have found a separating axis indicating that the
        # polygons do not intersect.
        if !overlap([minimum(p1r[1,:]), maximum(p1r[1,:])],
                    [minimum(p2r[1,:]), maximum(p2r[1,:])])
            return false
        end
    end

    return true        # No separating axis found
end

#--------------------------------------------------------------------------
"""
isconvexpoly - Test if 2D polygon is convex.

```
Usage: v = isconvexpoly(p)

Argument:  p - 2D polygon are defined as 2xN array of vertices.

Returns:   true if polygon is convex.

```

See also: convexpolyintersect(), pointinconvexpoly()
"""
function isconvexpoly(poly::Array)

    (dim, N) = size(poly)

    if dim != 2
        error("Polygon data must be a 2xN array")
    end

    # Strategy: Determine turning direction as one travels through
    # each set of three successive vertices.  If, for all vertices,
    # the turning direction is always clockwise or straight, or always
    # anti-clockwise or straight then the polygon is convex.

    c = zeros(Int, N)
    for n = 1:N-2
        c[n] = clockwise(view(poly,:,n), view(poly,:,n+1), view(poly,:,n+2))
    end
    c[N-1] = clockwise(view(poly,:,N-1), view(poly,:,N), view(poly,:,1))
    c[N] = clockwise(view(poly,:,N), view(poly,:,1), view(poly,:,2))

    if all(c .>= 0) || all(c .<= 0)
        return true
    else
        return false
    end
end


text:
struct NestedSimple
    x::Int64
end

struct Simple
    x::String
    y::NestedSimple
end

struct WithOptional
    x::Union{Nothing, Int64}
end

struct StructType
    x::Int64
    y::Float64
end

StructTypes.StructType(::Type{StructType}) = StructTypes.Struct()

struct Evolved1
    x::Int64
end

LightBSON.bson_schema_version(::Type{Evolved1}) = Int32(1)

struct Evolved2
    x::Int64
    y::Int64
end

Evolved2(x::Evolved1) = Evolved2(x.x, 0)

LightBSON.bson_schema_version(::Type{Evolved2}) = Int32(2)

struct Evolved3
    x::Int64
    y::Float64
end

Evolved3(x::Evolved2) = Evolved3(x.x, x.y)

LightBSON.bson_schema_version(::Type{Evolved3}) = Int32(3)

function LightBSON.bson_read_versioned(::Type{Evolved3}, v::Int32, reader::AbstractBSONReader)
    if v == 1
        Evolved3(Evolved2(bson_read_unversioned(Evolved1, reader)))
    elseif v == 2
        Evolved3(bson_read_unversioned(Evolved2, reader))
    elseif v == 3
        bson_read_unversioned(Evolved3, reader)
    else
        error("Unsupported schema version $v")
    end
end

struct Evolved4
    y::Float64
end

LightBSON.bson_schema_version(::Type{Evolved4}) = Int32(4)

const Evolved = Evolved3

@testset "struct" begin

@testset "simple" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    x = Simple("test", NestedSimple(123))
    writer[] = x
    close(writer)
    reader = BSONReader(buf)
    @test reader[Simple] == x
end

@testset "optional empty" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    x = WithOptional(nothing)
    writer[] = x
    close(writer)
    reader = BSONReader(buf)
    @test reader[WithOptional] == x
end

@testset "optional set" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    x = WithOptional(123)
    writer[] = x
    close(writer)
    reader = BSONReader(buf)
    @test reader[WithOptional] == x
end

@testset "StructType" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    x = StructType(123, 1.25)
    writer[] = x
    close(writer)
    reader = BSONReader(buf)
    @test reader[StructType] == x
end

@testset "evolved prev" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer[] = Evolved2(123, 456)
    close(writer)
    reader = BSONReader(buf)
    @test reader[Evolved] == Evolved3(123, 456)
end

@testset "evolved prev prev" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer[] = Evolved1(123)
    close(writer)
    reader = BSONReader(buf)
    @test reader[Evolved] == Evolved3(123, 0)
end

@testset "evolved current" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer[] = Evolved3(123, 1.25)
    close(writer)
    reader = BSONReader(buf)
    @test reader[Evolved] == Evolved3(123, 1.25)
end

@testset "evolved next" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer[] = Evolved4(1.25)
    close(writer)
    reader = BSONReader(buf)
    @test_throws ErrorException reader[Evolved]
end

@testset "NamedTuple" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    x = (; x = 123, y = 456, z = (; a = 1.25, b = "test"))
    writer[] = x
    close(writer)
    reader = BSONReader(buf)
    reader[typeof(x)] == x
end

end


text:
"""
# module LaminartParameterFunc

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Functions to create named tuple containing parameters, and input image in approprate type and dimensions for NNlib.conv/imfilter and CPU/GPU implementation.
"""
module LaminartInitFunc

include("./LaminartEqConv.jl")
include("./LaminartEqImfilter.jl")
include("./LaminartKernels.jl")

using NNlib, ImageFiltering, Images, OffsetArrays, CUDA, Noise

"""
Returns named tuple with input image, retina, kernels and parameters.
For use with equations with NNlib conv on GPU.
"""
function parameterInit_conv_gpu(imgLoc::String, p::NamedTuple)
    img = convert(Array{Float32,2}, load(imgLoc))
    img = reshape2d_4d(img)
    img = cu(img)

    r = similar(img)

    parameters = kernels_conv_gpu(img, p)

    LaminartEqConv.I_u!(r, img, parameters)
    temp_out = (I = img, r = r)
    parameters = merge(parameters, temp_out)
    return parameters
end

"""
Adds multiply Gaussian noise to input image

Returns named tuple for GPU with input image, retina, kernels and parameters.
For use with equations with NNlib conv on GPU.
"""
function parameterInit_conv_gpu_noise(
    imgLoc::String,
    p::NamedTuple,
    noise::Real,
)
    img = convert(Array{Float32,2}, load(imgLoc))
    img = mult_gauss(img, noise)
    img = reshape2d_4d(img)
    img = cu(img)

    r = similar(img)

    parameters = kernels_conv_gpu(img, p)

    LaminartEqConv.I_u!(r, img, parameters)
    temp_out = (I = img, r = r)
    parameters = merge(parameters, temp_out)
    return parameters
end

"""
Returns named tuple  with input image, retina, kernels and parameters.
For use with equations with NNlib conv on CPU.
"""
function parameterInit_conv_cpu(imgLoc::String, p::NamedTuple)
    img = convert(Array{Float32,2}, load(imgLoc))
    img = reshape2d_4d(img)

    r = similar(img)

    parameters = kernels_conv_cpu(img, p)

    LaminartEqConv.I_u!(r, img, parameters)
    temp_out = (I = img, r = r)
    parameters = merge(parameters, temp_out)
    return parameters
end

"""
Returns named tuple  with input image, retina, kernels and parameters.
For use with equations with imfilter on CPU.
"""
function parameterInit_imfil_cpu(imgLoc::String, p::NamedTuple)
    img = convert(Array{Float32,2}, load(imgLoc))

    r = similar(img)

    parameters = kernels_imfil_cpu(img, p)

    LaminartEqImfilter.I_u!(r, img, parameters)
    temp_out = (I = img, r = r)
    parameters = merge(parameters, temp_out)
    return parameters
end

"""
Reshapes ixj array to ixjx1x1 array.
"""
function reshape2d_4d(img::AbstractArray)
    reshape(img, size(img)[1], size(img)[2], 1, 1)
end

function add_I_u_p(I::AbstractArray, p::NamedTuple)
    # 	todo fix
    # 	I_4d = cu(reshape2d_4d(I))
    I_4d = reshape2d_4d(I)
    r = similar(I_4d)
    I_u!(r, I_4d, p)
    temp_out = (I = I_4d, r = r)
    return merge(p, temp_out)
end

"""
Generates kernels for use with NNLib conv and GPU.
"""
function kernels_conv_gpu(img::AbstractArray, p::NamedTuple)
    C_Q_temp = reshape(
        Array{eltype(img)}(undef, p.C_AB_l, p.C_AB_l * p.K),
        p.C_AB_l,
        p.C_AB_l,
        1,
        p.K,
    )
    C_P_temp = similar(C_Q_temp)
    H_temp = reshape(
        zeros(eltype(img), p.H_l, p.H_l * p.K * p.K),
        p.H_l,
        p.H_l,
        p.K,
        p.K,
    )
    T_temp = reshape(Array{eltype(img)}(undef, p.K * p.K), 1, 1, p.K, p.K)
    W_p_temp = reshape(
        Array{eltype(img)}(undef, p.W_l, p.W_l * p.K * p.K),
        p.W_l,
        p.W_l,
        p.K,
        p.K,
    )
    W_m_temp = similar(W_p_temp)

    for k ∈ 1:p.K
        θ = π * (k - 1.0f0) / p.K
        C_Q_temp[:, :, 1, k] = LaminartKernels.kern_d(p.σ_2, θ)
        C_P_temp[:, :, 1, k] = LaminartKernels.kern_b(p.σ_2, θ)
        H_temp[:, :, k, k] =
            p.H_fact .* LaminartKernels.gaussian_rot(p.H_σ_x, p.H_σ_y, θ, p.H_l)
        # 		todo make T kernel more general for higher K
        T_temp[1, 1, k, 1] = p.T_fact[k]
        T_temp[1, 1, 2, 2] = p.T_fact[1]
        T_temp[1, 1, 1, 2] = p.T_fact[2]
        #todo: generalise T and W for higher K
        #         T_temp[:,:,k] = KernelFactors.gaussian(p.T_σ, p.K)
        #         for l ∈ 1:p.K
        #             W_temp[:,:,l,k] =
        #         end
    end

    W_p_temp[:, :, 1, 1] =
        p.W_p_same_fact .* LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_a,
            p.W_p_σ_y_same_a,
            0f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_b,
            p.W_p_σ_y_same_b,
            0f0,
            p.W_l,
        )

    W_p_temp[:, :, 2, 2] =
        p.W_p_same_fact .* LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_a,
            p.W_p_σ_y_same_a,
            π / 2f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_b,
            p.W_p_σ_y_same_b,
            π / 2f0,
            p.W_l,
        )

    W_p_temp[:, :, 1, 2] =
        relu.(
            p.W_p_opp_fact_a .*
            (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_p_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_b,
                    p.W_p_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_c,
                    p.W_p_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    W_p_temp[:, :, 2, 1] =
        relu.(
            p.W_p_opp_fact_a .*
            (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_p_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_b,
                    p.W_p_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_c,
                    p.W_p_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    W_m_temp[:, :, 1, 1] =
        p.W_m_same_fact .* LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_a,
            p.W_m_σ_y_same_a,
            0f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_b,
            p.W_m_σ_y_same_b,
            0f0,
            p.W_l,
        )

    W_m_temp[:, :, 2, 2] =
        p.W_m_same_fact .* LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_a,
            p.W_m_σ_y_same_a,
            π / 2f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_b,
            p.W_m_σ_y_same_b,
            π / 2f0,
            p.W_l,
        )

    W_m_temp[:, :, 1, 2] =
        relu.(
            p.W_m_opp_fact_a .*
            (Kernel.gaussian((p.W_m_σ_opp_a, p.W_m_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_m_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_m_σ_x_opp_b,
                    p.W_m_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_m_σ_x_opp_c,
                    p.W_m_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    W_m_temp[:, :, 2, 1] =
        relu.(
            p.W_p_opp_fact_a .*
            (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_p_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_b,
                    p.W_p_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_c,
                    p.W_p_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    temp_out = (
        k_gauss_1 = cu(reshape2d_4d(Kernel.gaussian(p.σ_1))),
        k_gauss_2 = cu(reshape2d_4d(Kernel.gaussian(p.σ_2))),
        k_C_d = cu(C_Q_temp),
        k_C_b = cu(C_P_temp),

        # 		todo use mean of x_lgn?
        k_x_lgn = cu(reshape(ones(Float32, 1, p.K), 1, 1, p.K, 1)),
        # 		k_x_lgn = cu(reshape(ones(Float32,1,p.K),1,1,p.K,1) ./ p.K),
        k_W_p = cu(W_p_temp),
        k_W_m = cu(W_m_temp),
        k_H = cu(H_temp),
        k_T_p = cu(T_temp),
        k_T_m = cu(p.T_p_m .* T_temp),
        k_T_p_v2 = cu(p.T_v2_fact .* T_temp),
        k_T_m_v2 = cu(p.T_v2_fact .* p.T_p_m .* T_temp),
        dim_i = size(img)[1],
        dim_j = size(img)[2],
        x_V2 = cu(reshape(
            zeros(Float32, size(img)[1], size(img)[2] * p.K),
            size(img)[1],
            size(img)[2],
            p.K,
            1,
        )),
        ν_pw_n = p.ν^p.n,
    )

    merge(p, temp_out)
end

"""
Generates kernels for use with NNLib conv and CPU.
"""
function kernels_conv_cpu(img::AbstractArray, p::NamedTuple)
    C_Q_temp = reshape(
        Array{eltype(img)}(undef, p.C_AB_l, p.C_AB_l * p.K),
        p.C_AB_l,
        p.C_AB_l,
        1,
        p.K,
    )
    C_P_temp = similar(C_Q_temp)
    H_temp = reshape(
        zeros(eltype(img), p.H_l, p.H_l * p.K * p.K),
        p.H_l,
        p.H_l,
        p.K,
        p.K,
    )
    T_temp = reshape(Array{eltype(img)}(undef, p.K * p.K), 1, 1, p.K, p.K)
    W_p_temp = reshape(
        Array{eltype(img)}(undef, p.W_l, p.W_l * p.K * p.K),
        p.W_l,
        p.W_l,
        p.K,
        p.K,
    )
	W_m_temp = reshape(
        Array{eltype(img)}(undef, p.W_l, p.W_l * p.K * p.K),
        p.W_l,
        p.W_l,
        p.K,
        p.K,
    )
    for k ∈ 1:p.K
        θ = π * (k - 1.0f0) / p.K
        C_Q_temp[:, :, 1, k] = LaminartKernels.kern_d(p.σ_2, θ)
        C_P_temp[:, :, 1, k] = LaminartKernels.kern_b(p.σ_2, θ)
        H_temp[:, :, k, k] =
            p.H_fact .* LaminartKernels.gaussian_rot(p.H_σ_x, p.H_σ_y, θ, p.H_l)
        # 		todo make T kernel more general for higher K
        T_temp[1, 1, k, 1] = p.T_fact[k]
        T_temp[1, 1, 2, 2] = p.T_fact[1]
        T_temp[1, 1, 1, 2] = p.T_fact[2]
        #todo: generalise T and W for higher K
        #         T_temp[:,:,k] = KernelFactors.gaussian(p.T_σ, p.K)
        #         for l ∈ 1:p.K
        #             W_temp[:,:,l,k] =
        #         end
    end

    W_p_temp[:, :, 1, 1] =
        p.W_p_same_fact .* LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_a,
            p.W_p_σ_y_same_a,
            0f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_b,
            p.W_p_σ_y_same_b,
            0f0,
            p.W_l,
        )

    W_p_temp[:, :, 2, 2] =
        p.W_p_same_fact .* LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_a,
            p.W_p_σ_y_same_a,
            π / 2f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_b,
            p.W_p_σ_y_same_b,
            π / 2f0,
            p.W_l,
        )

    W_p_temp[:, :, 1, 2] =
        relu.(
            p.W_p_opp_fact_a .*
            (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_p_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_b,
                    p.W_p_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_c,
                    p.W_p_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    W_p_temp[:, :, 2, 1] =
        relu.(
            p.W_p_opp_fact_a .*
            (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_p_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_b,
                    p.W_p_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_c,
                    p.W_p_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    W_m_temp[:, :, 1, 1] =
        p.W_m_same_fact .* LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_a,
            p.W_m_σ_y_same_a,
            0f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_b,
            p.W_m_σ_y_same_b,
            0f0,
            p.W_l,
        )

    W_m_temp[:, :, 2, 2] =
        p.W_m_same_fact .* LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_a,
            p.W_m_σ_y_same_a,
            π / 2f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_b,
            p.W_m_σ_y_same_b,
            π / 2f0,
            p.W_l,
        )

    W_m_temp[:, :, 1, 2] =
        relu.(
            p.W_m_opp_fact_a .*
            (Kernel.gaussian((p.W_m_σ_opp_a, p.W_m_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_m_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_m_σ_x_opp_b,
                    p.W_m_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_m_σ_x_opp_c,
                    p.W_m_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    W_m_temp[:, :, 2, 1] =
        relu.(
            p.W_p_opp_fact_a .*
            (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
            p.W_p_opp_fact_b .* (
                LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_b,
                    p.W_p_σ_y_opp_b,
                    0f0,
                    p.W_l,
                ) .+ LaminartKernels.gaussian_rot(
                    p.W_p_σ_x_opp_c,
                    p.W_p_σ_y_opp_c,
                    0f0,
                    p.W_l,
                )
            ),
        )

    temp_out = (
        k_gauss_1 = reshape2d_4d(Kernel.gaussian(p.σ_1)),
        k_gauss_2 = reshape2d_4d(Kernel.gaussian(p.σ_2)),
        k_C_d = C_Q_temp,
        k_C_b = C_P_temp,

        # 		todo use mean of x_lgn?
        k_x_lgn = reshape(ones(Float32, 1, p.K), 1, 1, p.K, 1),
        # 		k_x_lgn = reshape(ones(Float32,1,p.K),1,1,p.K,1)./p.K,
        k_W_p = W_p_temp,
        k_W_m = W_m_temp,
        k_H = H_temp,
        k_T_p = T_temp,
        k_T_m = p.T_p_m .* T_temp,
        k_T_p_v2 = p.T_v2_fact .* T_temp,
        k_T_m_v2 = p.T_v2_fact .* p.T_p_m .* T_temp,
        dim_i = size(img)[1],
        dim_j = size(img)[2],
        x_V2 = reshape(
            zeros(Float32, size(img)[1], size(img)[2] * p.K),
            size(img)[1],
            size(img)[2],
            p.K,
            1,
        ),
        ν_pw_n = p.ν^p.n,
    )
    merge(p, temp_out)
end

"""
Generates kernels for use with imfilter equations.
"""
function kernels_imfil_cpu(img::AbstractArray, p::NamedTuple)
    C_Q_temp = reshape(
        Array{eltype(img)}(undef, p.C_AB_l, p.C_AB_l * p.K),
        p.C_AB_l,
        p.C_AB_l,
        p.K,
    )
    C_P_temp = copy(C_Q_temp)
    H_temp = reshape(
        Array{eltype(img)}(undef, p.H_l, p.H_l * p.K),
        p.H_l,
        p.H_l,
        p.K,
    )
    T_temp = reshape(Array{eltype(img)}(undef, 1, 1 * p.K), 1, 1, p.K)     #ijk,  1x1xk,   ijk
    W_p_temp = reshape(
        Array{eltype(img)}(undef, p.W_l, p.W_l * p.K * p.K),
        p.W_l,
        p.W_l,
        p.K,
        p.K,
    )
	W_m_temp = reshape(
        Array{eltype(img)}(undef, p.W_l, p.W_l * p.K * p.K),
        p.W_l,
        p.W_l,
        p.K,
        p.K,
    )    #ijk,  1x1xk,   ijk
    for k ∈ 1:p.K
        θ = π * (k - 1) / p.K
        C_Q_temp[:, :, k] = reflect(centered(LaminartKernels.kern_d(p.σ_2, θ)))           #ij ijk ijk
        C_P_temp[:, :, k] = reflect(centered(LaminartKernels.kern_b(p.σ_2, θ)))               #ij ijk ijk
        H_temp[:, :, k] = reflect(
            p.H_fact .*
            LaminartKernels.gaussian_rot(p.H_σ_x, p.H_σ_y, θ, p.H_l),
        )  #ijk, ij for each k; ijk
        T_temp[:, :, k] = reshape([p.T_fact[k]], 1, 1)
        #todo: generalise T and W for higher K
        #         T_temp[:,:,k] = KernelFactors.gaussian(p.T_σ, p.K)
        #         for l ∈ 1:p.K
        #             W_temp[:,:,l,k] =
        #         end
    end
    W_p_temp[:, :, 1, 1] = reflect(
        p.W_p_same_fact .* LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_a,
            p.W_p_σ_y_same_a,
            0f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_b,
            p.W_p_σ_y_same_b,
            0f0,
            p.W_l,
        ),
    )

    W_p_temp[:, :, 2, 2] = reflect(
        p.W_p_same_fact .* LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_a,
            p.W_p_σ_y_same_a,
            π / 2f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_p_σ_x_same_b,
            p.W_p_σ_y_same_b,
            π / 2f0,
            p.W_l,
        ),
    )

    W_p_temp[:, :, 1, 2] = reflect(relu.(
        p.W_p_opp_fact_a .*
        (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
        p.W_p_opp_fact_b .* (
            LaminartKernels.gaussian_rot(
                p.W_p_σ_x_opp_b,
                p.W_p_σ_y_opp_b,
                0f0,
                p.W_l,
            ) .+ LaminartKernels.gaussian_rot(
                p.W_p_σ_x_opp_c,
                p.W_p_σ_y_opp_c,
                0f0,
                p.W_l,
            )
        ),
    ))

    W_p_temp[:, :, 2, 1] = reflect(relu.(
        p.W_p_opp_fact_a .*
        (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
        p.W_p_opp_fact_b .* (
            LaminartKernels.gaussian_rot(
                p.W_p_σ_x_opp_b,
                p.W_p_σ_y_opp_b,
                0f0,
                p.W_l,
            ) .+ LaminartKernels.gaussian_rot(
                p.W_p_σ_x_opp_c,
                p.W_p_σ_y_opp_c,
                0f0,
                p.W_l,
            )
        ),
    ))

    W_m_temp[:, :, 1, 1] = reflect(
        p.W_m_same_fact .* LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_a,
            p.W_m_σ_y_same_a,
            0f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_b,
            p.W_m_σ_y_same_b,
            0f0,
            p.W_l,
        ),
    )

    W_m_temp[:, :, 2, 2] = reflect(
        p.W_m_same_fact .* LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_a,
            p.W_m_σ_y_same_a,
            π / 2f0,
            p.W_l,
        ) .+ LaminartKernels.gaussian_rot(
            p.W_m_σ_x_same_b,
            p.W_m_σ_y_same_b,
            π / 2f0,
            p.W_l,
        ),
    )

    W_m_temp[:, :, 1, 2] = reflect(relu.(
        p.W_m_opp_fact_a .*
        (Kernel.gaussian((p.W_m_σ_opp_a, p.W_m_σ_opp_a), (p.W_l, p.W_l))) .-
        p.W_m_opp_fact_b .* (
            LaminartKernels.gaussian_rot(
                p.W_m_σ_x_opp_b,
                p.W_m_σ_y_opp_b,
                0f0,
                p.W_l,
            ) .+ LaminartKernels.gaussian_rot(
                p.W_m_σ_x_opp_c,
                p.W_m_σ_y_opp_c,
                0f0,
                p.W_l,
            )
        ),
    ))

    W_m_temp[:, :, 2, 1] = reflect(relu.(
        p.W_p_opp_fact_a .*
        (Kernel.gaussian((p.W_p_σ_opp_a, p.W_p_σ_opp_a), (p.W_l, p.W_l))) .-
        p.W_p_opp_fact_b .* (
            LaminartKernels.gaussian_rot(
                p.W_p_σ_x_opp_b,
                p.W_p_σ_y_opp_b,
                0f0,
                p.W_l,
            ) .+ LaminartKernels.gaussian_rot(
                p.W_p_σ_x_opp_c,
                p.W_p_σ_y_opp_c,
                0f0,
                p.W_l,
            )
        ),
    ))

    # todo fix W kernel
    #  W_temp[:,:,1,1] = reflect(LaminartKernels.gaussian_rot(3,0.8,0,19))
    #     W_temp[:,:,2,2] = reflect(LaminartKernels.gaussian_rot(3,0.8,0,19))
    #     W_temp[:,:,1,2] = reflect(LaminartKernels.gaussian_rot(3,0.8,0,19))
    #     W_temp[:,:,2,1] = reflect(LaminartKernels.gaussian_rot(3,0.8,0,19))

    # todo: fix range of W H
    #     W_range = -(p.W_size-1)/2:(p.W_size-1)/2
    #     H_range = -(p.H_size-1)/2:(p.H_size-1)/2
    W_range = -9:9
    H_range = -9:9

    temp_out = (
        k_gauss_1 = reflect(Kernel.gaussian(p.σ_1)),
        k_gauss_2 = reflect(Kernel.gaussian(p.σ_2)),
        k_C_d = C_Q_temp,
        k_C_b = C_P_temp,
        k_W_p = W_p_temp,
        k_W_m = W_m_temp,
        # k_W_m = OffsetArray(W_temp, W_range, W_range, 1:p.K, 1:p.K),
        k_H = H_temp,
        # k_H = OffsetArray(H_temp, H_range, H_range, 1:p.K),
        k_T_p = T_temp,
        k_T_m = (p.T_p_m .* T_temp),
        k_T_p_v2 = (p.T_v2_fact .* T_temp),
        k_T_m_v2 = (p.T_v2_fact .* p.T_p_m .* T_temp),
        dim_i = size(img)[1],
        dim_j = size(img)[2],
        x_V2 = reshape(
            zeros(typeof(img[1, 1]), size(img)[1], size(img)[2] * p.K),
            size(img)[1],
            size(img)[2],
            p.K,
        ),
    )
    return merge(p, temp_out)
end

end


text:
export iplot, iplot!, iscatter, iscatter!, iheatmap

const bgcolor = RGB(1,1,1)

iplot(y::AbstractVecOrMat, args...; kwargs...) = iplot(1:size(y,1), y, args...; kwargs...)
iplot!(y::AbstractVecOrMat; kwargs...) = iplot!(1:size(y,1), y; kwargs...)

function iplot(x::AbstractVector, y::AbstractMatrix, args...; axes=true, xlabel=missing, ylabel=missing, axescolor=:black, grid=false, overlay=false, cursor=false, kwargs...)
  colors = distinguishable_colors(size(y,2)+1, [bgcolor])
  iplot(x, y[:,1], args...; color=colors[2], axes=axes, xlabel=xlabel, ylabel=ylabel, axescolor=axescolor, grid=grid, overlay=overlay, cursor=cursor, kwargs...)
  for k ∈ 2:size(y,2)
    iplot!(x, y[:,k]; color=colors[1+k], kwargs...)
  end
end

function iplot!(x::AbstractVector, y::AbstractMatrix; kwargs...)
  colors = distinguishable_colors(size(y,2)+1, [bgcolor])
  for k ∈ 1:size(y,2)
    iplot!(x, y[:,k]; color=colors[1+k], kwargs...)
  end
end

function iplot(x::AbstractVector, y::AbstractVector, x1=missing, x2=missing, y1=missing, y2=missing; axes=true, xlabel=missing, ylabel=missing, axescolor=:black, grid=false, overlay=false, cursor=false, kwargs...)
  length(x) != size(y,1) && error("x and y must be of equal length")
  if x1 === missing || x2 === missing
    x1a, x2a = autorange(x)
    x1 === missing && (x1 = x1a)
    x2 === missing && (x2 = x2a)
  end
  if y1 === missing || y2 === missing
    y1a, y2a = autorange(y, 0.1)
    y1 === missing && (y1 = y1a)
    y2 === missing && (y2 = y2a)
  end
  x2 <= x1 && error("Bad x range")
  y2 <= y1 && error("Bad y range")
  datasrc = datasource(;
    generate! = x isa AbstractRange ? linepool! : linecrop!,
    xyrect = ℛ{Float64}(x1, y1, x2-x1, y2-y1),
    xmin = minimum(x),
    xmax = maximum(x),
    data = hcat(x, y),
    kwargs...
  )
  viz = ifigure()
  ijrect = axes && !overlay ? inset(viz.ijrect, 100, 50, 100, 50) : viz.ijrect
  c = addcanvas!(LineCanvas, viz, datasrc; rect=ijrect, kwargs...)
  axes && addaxes!(c; xlabel=xlabel, ylabel=ylabel, color=axescolor, grid=grid, inset = overlay ? 100 : 0, border = overlay ? 0 : 150)
  cursor && addcursor!(c, color=axescolor)
  c
end

function iplot!(x::AbstractVector, y::AbstractVector; kwargs...)
  viz = ifigure(hold=true)
  length(viz.children) == 0 && error("No previous canvas available to plot over")
  prev = viz.children[end]
  datasrc = datasource(;
    generate! = x isa AbstractRange ? linepool! : linecrop!,
    xyrect = prev.datasrc.xyrect,
    xylock = prev.datasrc.xylock,
    data = hcat(x, y),
    kwargs...
  )
  c = addcanvas!(LineCanvas, viz, datasrc; rect=prev.ijrect, kwargs...)
  c.xyrect[] = prev.xyrect[]
  c
end

function iplot(f::Function, x1=0.0, x2=1.0, y1=missing, y2=missing; axes=true, xlabel=missing, ylabel=missing, axescolor=:black, grid=false, overlay=false, cursor=false, kwargs...)
  x2 <= x1 && error("Bad x range")
  if y1 === missing || y2 === missing
    y1a, y2a = autorange(f.(x1 .+ rand(1000)*(x2-x1)), 0.1)
    y1 === missing && (y1 = y1a)
    y2 === missing && (y2 = y2a)
  end
  y2 <= y1 && error("Bad y range")
  datasrc = datasource(
    generate! = (b, c) -> apply!(f, b, c),
    xyrect = ℛ{Float64}(x1, y1, x2-x1, y2-y1),
    kwargs...
  )
  viz = ifigure()
  ijrect = axes && !overlay ? inset(viz.ijrect, 100, 50, 100, 50) : viz.ijrect
  c = addcanvas!(LineCanvas, viz, datasrc; rect=ijrect, kwargs...)
  axes && addaxes!(c; xlabel=xlabel, ylabel=ylabel, color=axescolor, grid=grid, inset = overlay ? 100 : 0, border = overlay ? 0 : 150)
  cursor && addcursor!(c, color=axescolor)
  c
end

function iplot!(f::Function; kwargs...)
  viz = ifigure(hold=true)
  length(viz.children) == 0 && error("No previous canvas available to plot over")
  prev = viz.children[end]
  datasrc = datasource(;
    generate! = (b, c) -> apply!(f, b, c),
    xyrect = prev.datasrc.xyrect,
    xylock = prev.datasrc.xylock,
    kwargs...
  )
  c = addcanvas!(LineCanvas, viz, datasrc; rect=prev.ijrect, kwargs...)
  c.xyrect[] = prev.xyrect[]
  c
end

function iscatter(x::AbstractVector, y::AbstractVector, x1=missing, x2=missing, y1=missing, y2=missing; aggregate=false, xylock=true, axes=true, xlabel=missing, ylabel=missing, axescolor=:black, grid=false, overlay=false, cursor=false, kwargs...)
  length(x) != length(y) && error("x and y must be of equal length")
  if x1 === missing || x2 === missing
    x1a, x2a = autorange(x)
    x1 === missing && (x1 = x1a)
    x2 === missing && (x2 = x2a)
  end
  if y1 === missing || y2 === missing
    y1a, y2a = autorange(y)
    y1 === missing && (y1 = y1a)
    y2 === missing && (y2 = y2a)
  end
  x2 <= x1 && error("Bad x range")
  y2 <= y1 && error("Bad y range")
  datasrc = datasource(;
    generate! = aggregate ? aggregate! : pointcrop!,
    xyrect = ℛ{Float64}(x1, y1, x2-x1, y2-y1),
    xylock = xylock,
    xmin = minimum(x),
    xmax = maximum(x),
    ymin = minimum(y),
    ymax = maximum(y),
    data = hcat(x, y),
    kwargs...
  )
  viz = ifigure()
  ijrect = axes && !overlay ? inset(viz.ijrect, 100, 50, 100, 50) : viz.ijrect
  if aggregate
    haskey(kwargs, :colormap) || (kwargs = merge(Dict((:colormap => :Greys)), kwargs))
    c = addcanvas!(HeatmapCanvas, viz, datasrc; rect=ijrect, kwargs...)
  else
    c = addcanvas!(ScatterCanvas, viz, datasrc; rect=ijrect, kwargs...)
  end
  axes && addaxes!(c; xlabel=xlabel, ylabel=ylabel, color=axescolor, grid=grid, inset = overlay ? 100 : 0, border = overlay ? 0 : 150)
  cursor && addcursor!(c, color=axescolor)
  c
end

function iscatter!(x::AbstractVector, y::AbstractVector; kwargs...)
  viz = ifigure(hold=true)
  length(viz.children) == 0 && error("No previous canvas available to plot over")
  prev = viz.children[end]
  datasrc = datasource(;
    generate! = pointcrop!,
    xyrect = prev.datasrc.xyrect,
    xylock = prev.datasrc.xylock,
    data = hcat(x, y),
    kwargs...
  )
  c = addcanvas!(ScatterCanvas, viz, datasrc; rect=prev.ijrect, kwargs...)
  c.xyrect[] = prev.xyrect[]
  c
end

# TODO: allow starting with a smaller view than the full z matrix
function iheatmap(z::AbstractMatrix, x1=0.0, x2=1.0, y1=0.0, y2=1.0; clim=missing, pooling=mean, axes=true, xlabel=missing, ylabel=missing, axescolor=:black, grid=false, overlay=false, cursor=false, kwargs...)
  x2 <= x1 && error("Bad x range")
  y2 <= y1 && error("Bad y range")
  clim === missing && (clim = extrema(z))
  viz = ifigure()
  datasrc = datasource(;
    generate! = (b, c) -> heatmappool!(b, c; pooling=pooling),
    xyrect = ℛ{Float64}(x1, y1, x2-x1, y2-y1),
    data = z,
    clim = clim,
    xmin = x1,
    xmax = x2,
    ymin = y1,
    ymax = y2,
    kwargs...
  )
  ijrect = axes && !overlay ? inset(viz.ijrect, 100, 50, 100, 50) : viz.ijrect
  c = addcanvas!(HeatmapCanvas, viz, datasrc; rect=ijrect, kwargs...)
  axes && addaxes!(c; xlabel=xlabel, ylabel=ylabel, grid=grid, color=axescolor, inset = overlay ? 100 : 0, border = overlay ? 0 : 150)
  cursor && addcursor!(c, color=axescolor)
  c
end

function iheatmap(f::Function, x1=0.0, x2=1.3, y1=0.0, y2=1.0; xylock=true, axes=true, xlabel=missing, ylabel=missing, axescolor=:black, grid=false, overlay=false, cursor=false, kwargs...)
  x2 <= x1 && error("Bad x range")
  y2 <= y1 && error("Bad y range")
  viz = ifigure()
  datasrc = datasource(;
    generate! = (b, c) -> apply!(f, b, c),
    xyrect = ℛ{Float64}(x1, y1, x2-x1, y2-y1),
    xylock = xylock,
    kwargs...
  )
  ijrect = axes && !overlay ? inset(viz.ijrect, 100, 50, 100, 50) : viz.ijrect
  c = addcanvas!(HeatmapCanvas, viz, datasrc; rect=ijrect, kwargs...)
  axes && addaxes!(c; xlabel=xlabel, ylabel=ylabel, grid=grid, color=axescolor, inset = overlay ? 100 : 0, border = overlay ? 0 : 150)
  cursor && addcursor!(c, color=axescolor)
  c
end

### helpers

function autorange(y, margin=0.0)
  y1 = minimum(y)
  y2 = maximum(y)
  yr = y2 - y1
  if yr == 0
    y1 -= 0.5
    y2 += 0.5
  elseif margin != 0
    y1 -= margin*yr
    y2 += margin*yr
  end
  y1, y2
end

inset(r::ℛ, li, ri, bi, ti) = ℛ(left(r) + li, bottom(r) + bi, width(r) - li - ri, height(r) - bi - ti)
inset(r::Node{ℛ{T}}, li, ri, bi, ti) where T = lift(r -> ℛ(left(r) + li, bottom(r) + bi, width(r) - li - ri, height(r) - bi - ti), r)

function datasource(; kwargs...)
  dargs = Dict()
  for k ∈ fieldnames(DataSource)
    k ∈ keys(kwargs) && (dargs[k] = kwargs[k])
  end
  DataSource(; dargs...)
end


text:
"""
    generate_population(icn, pop_size
Generate a pôpulation of weigths (individuals) for the genetic algorithm weigthing `icn`.
"""
function generate_population(icn, pop_size)
    population = Vector{BitVector}()
    foreach(_ -> push!(population, falses(nbits(icn))), 1:pop_size)
    return population
end

"""
    _optimize!(icn, X, X_sols; metric = hamming, pop_size = 200)
Optimize and set the weigths of an ICN with a given set of configuration `X` and solutions `X_sols`.
"""
function _optimize!(
    icn,
    solutions,
    non_sltns,
    dom_size,
    param,
    metric,
    pop_size,
    iterations;
    samples=nothing,
    memoize=false,
)
    inplace = zeros(dom_size, max_icn_length())
    _non_sltns = isnothing(samples) ? non_sltns : rand(non_sltns, samples)

    function fitness(w)
        compo = compose(icn, w)
        f = composition(compo)
        S = Iterators.flatten((solutions, _non_sltns))
        return sum(x -> abs(f(x; X=inplace, param, dom_size) - metric(x, solutions)), S) +
               regularization(icn) +
               weigths_bias(w)
    end
    _fitness = memoize ? (@memoize Dict memoize_fitness(w) = fitness(w)) : fitness

    _icn_ga = GA(;
        populationSize=pop_size,
        crossoverRate=0.8,
        epsilon=0.05,
        selection=tournament(2),
        crossover=singlepoint,
        mutation=flip,
        mutationRate=1.0,
    )

    pop = generate_population(icn, pop_size)
    r = Evolutionary.optimize(_fitness, pop, _icn_ga, Evolutionary.Options(; iterations))
    return weights!(icn, Evolutionary.minimizer(r))
end

"""
    optimize!(icn, X, X_sols, global_iter, local_iter; metric=hamming, popSize=100)
Optimize and set the weigths of an ICN with a given set of configuration `X` and solutions `X_sols`. The best weigths among `global_iter` will be set.
"""
function optimize!(
    icn,
    solutions,
    non_sltns,
    global_iter,
    iter,
    dom_size,
    param,
    metric,
    pop_size;
    sampler=nothing,
    memoize=false,
)
    results = Dictionary{BitVector,Int}()
    aux_results = Vector{BitVector}(undef, global_iter)
    nt = Base.Threads.nthreads()

    @info """Starting optimization of weights$(nt > 1 ? " (multithreaded)" : "")"""
    samples = isnothing(sampler) ? nothing : sampler(length(solutions) + length(non_sltns))
    @qthreads for i in 1:global_iter
        @info "Iteration $i"
        aux_icn = deepcopy(icn)
        _optimize!(
            aux_icn,
            solutions,
            non_sltns,
            dom_size,
            param,
            eval(metric),
            pop_size,
            iter;
            samples,
            memoize,
        )
        aux_results[i] = weigths(aux_icn)
    end
    foreach(bv -> incsert!(results, bv), aux_results)
    best = rand(findall(x -> x == maximum(results), results))
    weights!(icn, best)
    return best, results
end


text:
module TestConvert
    using Base.Test
    using CategoricalArrays

    pool = CategoricalPool([1, 2, 3])
    @test convert(CategoricalPool{Int, CategoricalArrays.DefaultRefType}, pool) === pool
    @test convert(CategoricalPool{Int}, pool) === pool
    @test convert(CategoricalPool, pool) === pool
    convert(CategoricalPool{Float64, UInt8}, pool)
    convert(CategoricalPool{Float64}, pool)
    convert(CategoricalPool, pool)

    v1 = CategoricalValue(1, pool)
    v2 = CategoricalValue(2, pool)
    v3 = CategoricalValue(3, pool)

    @test convert(Int32, v1) === Int32(1)
    @test convert(Int32, v2) === Int32(2)
    @test convert(Int32, v3) === Int32(3)

    @test convert(UInt8, v1) === 0x01
    @test convert(UInt8, v2) === 0x02
    @test convert(UInt8, v3) === 0x03

    @test convert(CategoricalValue, v1) === v1
    @test convert(CategoricalValue{Int}, v1) === v1
    @test convert(CategoricalValue{Int, CategoricalArrays.DefaultRefType}, v1) === v1
    @test convert(Any, v1) === v1

    @test get(v1) === 1
    @test get(v2) === 2
    @test get(v3) === 3

    @test promote(1, v1) === (1, 1)
    @test promote(1.0, v1) === (1.0, 1.0)
    @test promote(0x1, v1) === (1, 1)

    # Test that ordered property is preserved
    pool = CategoricalPool([1, 2, 3], true)
    @test convert(CategoricalPool{Float64, UInt8}, pool).ordered === true
end


text:
import FFTW

@doc raw"""
Determine the minimal grid size for the cubic basis set to be able to
represent product of orbitals (with the default `supersampling=2`).

Optionally optimize the grid afterwards for the FFT procedure by
ensuring factorization into small primes.

The function will determine the smallest cube containing the wave vectors
 ``|G|^2/2 \leq E_\text{cut} ⋅ \text{supersampling}^2``.
For an exact representation of the density resulting from wave functions
represented in the spherical basis sets, `supersampling` should be at least `2`.
"""
function determine_grid_size(lattice::AbstractMatrix{T}, Ecut; supersampling=2, tol=1e-8,
                             ensure_smallprimes=true) where T
    cutoff_Gsq = 2 * supersampling^2 * Ecut
    Gmax = [norm(lattice[:, i]) / 2T(π) * sqrt(cutoff_Gsq) for i in 1:3]

    # Round up, unless exactly zero (in which case keep it zero in
    # order to just have one G vector for 1D or 2D systems)
    for i = 1:3
        if Gmax[i] != 0
            Gmax[i] = ceil.(Int, Gmax[i] .- tol)
        end
    end

    # Optimize FFT grid size: Make sure the number factorises in small primes only
    if ensure_smallprimes
        Vec3([nextprod([2, 3, 5], 2gs + 1) for gs in Gmax])
    else
        Vec3([2gs + 1 for gs in Gmax])
    end
end
function determine_grid_size(model::Model, Ecut; kwargs...)
    determine_grid_size(model.lattice, Ecut; kwargs...)
end

"""
Plan a FFT of type `T` and size `fft_size`, spending some time on finding an
optimal algorithm. Both an inplace and an out-of-place FFT plan are returned.
"""
function build_fft_plans(T::Type{Float32}, fft_size)
    # TODO For Float32 there are issues with aligned FFTW plans, so we
    #      fall back to unaligned FFTW plans (which are generally discouraged).
    tmp = Array{ComplexF32}(undef, fft_size...)
    ipFFT = FFTW.plan_fft!(tmp, flags=FFTW.MEASURE | FFTW.UNALIGNED)
    opFFT = FFTW.plan_fft(tmp, flags=FFTW.MEASURE | FFTW.UNALIGNED)
    ipFFT, opFFT
end

function build_fft_plans(T::Type{Float64}, fft_size)
    tmp = Array{ComplexF64}(undef, fft_size...)
    ipFFT = FFTW.plan_fft!(tmp, flags=FFTW.MEASURE)
    opFFT = FFTW.plan_fft(tmp, flags=FFTW.MEASURE)
    ipFFT, opFFT
end

# TODO Some grid sizes are broken in the generic FFT implementation
# in FourierTransforms, for more details see fft_generic.jl
# This function is needed to provide a noop fallback for grid adjustment for
# for floating-point types natively supported by FFTW
next_working_fft_size(::Type{Float32}, size) = size
next_working_fft_size(::Type{Float64}, size) = size


text:
module AssociativeArray

import Base: indices, pairs, getindex, setindex!, similar, start, next, done, length, size,
             pairs, IndexStyle

using Base: OneTo, @propagate_inbounds

export Assoc, Dic, AbstractArr, AbstractVec, AbstractMat, Arr, Vec, Mat
export indextype, empty, issetequal
export IndexDirect, IndexToken, IndexStyle # for my own convenience

eltype(::Base.KeyIterator{<:Associative{K, V}}) where {K, V} = K

include("issetequal.jl")
include("Associative.jl")
include("Dict.jl")
include("AbstractArray.jl")
include("Array.jl")

end # module


text:
function getJointsAndObject3DsWithoutParents!(evaluatedParameters,
                                              object3DWithoutParents::Vector{Object3D},
                                              jointObjects::Vector{Object3D},
                                              path::String)::Nothing
    for (key,value) in evaluatedParameters   # zip(keys(evaluatedParameters), evaluatedParameters)
        if typeof(value) == Object3D
            if value.parent === value
                push!(object3DWithoutParents, value)
            elseif typeof(value.feature) == Modia3D.Composition.SceneOptions
                error(value.path, ": Object3D has a parent and defines SceneOptions!")
            end

        elseif typeof(value) <: Modia3D.AbstractJoint
            push!(jointObjects, value.obj2)

        elseif typeof(value) <: OrderedDict
            getJointsAndObject3DsWithoutParents!(value, object3DWithoutParents, jointObjects, path)
        end
    end
    return nothing
end

multiBodyName(instantiatedModel, mbsPath) = mbsPath == "" ? instantiatedModel.modelName :
                                                            instantiatedModel.modelName * "." * mbsPath

"""
    (world, jointObjects) = checkMultibodySystemAndGetWorldAndJoints(instantiatedModel, id)

Recursively traverse model and perform the following actions:

- Search for an OrderedDict that has `key = :_id, value = id`
  (identifies the root of the multibody system).

- Search from the root of the multibody system and perform the following actions:
  - Check that from all Object3Ds exactly one of them has no parent.
    Return this parent as `world`.
  - Check that only `world` has potentially a `feature` entry that
    is a SceneOption.
  - Return a vector of joint objects as `joints`.
"""
function checkMultibodySystemAndGetWorldAndJoints(instantiatedModel::ModiaLang.SimulationModel{FloatType,ParType,EvaluatedParType,TimeType}, id::Int) where {FloatType,ParType,EvaluatedParType,TimeType}
    # Find root mbs of multibody system
    (mbsRoot, mbsPath) = ModiaLang.getIdParameter(instantiatedModel.evaluatedParameters, ParType, id)
    if isnothing(mbsRoot)
        error(instantiatedModel.modelName, ": did not find _id = ", id, " in the evaluated parameters!")
    end

    object3DWithoutParents = Object3D[]
    jointObjects = Object3D[]
    getJointsAndObject3DsWithoutParents!(mbsRoot, object3DWithoutParents, jointObjects, mbsPath)

    if length(object3DWithoutParents) == 0
        error("\n", multiBodyName(instantiatedModel, mbsPath), ": There is either no Object3D or all of them have a parent!\n",
              "(Note, there must be exactly one Object3D that has no parent.)")
    elseif length(object3DWithoutParents) > 1
        object3DNames = "   " * object3DWithoutParents[1].path
        for i = 2:length(object3DWithoutParents)
            object3DNames *= "\n   " * object3DWithoutParents[i].path
        end
        error("\n", instantiatedModel.modelName, ": The following ", length(object3DWithoutParents), " Object3Ds have no parent\n",
              "(note, there must be exactly one Object3D that has no parent):\n", object3DNames, "\n")
    end
    return (object3DWithoutParents[1], jointObjects)
end

function initAnalysis2!(world)
    # Construct Scene(..) object
    Modia3D.Composition.EmptyObject3DFeature
    if typeof(world.feature) <: Modia3D.SceneOptions
        scene = Modia3D.Scene(world.feature)
    else
        scene = Modia3D.Scene()
    end
    scene.analysis = Modia3D.DynamicAnalysis

    Modia3D.Composition.chooseAndBuildUpTree(world, scene)

    if !scene.options.enableContactDetection
        scene.collide = false
    end

    return scene
end

struct MultibodyData{FloatType}
    nqdd::Int                       # Length of qdd vector
    world::Object3D                 # Pointer to world object
    scene::Scene                    # Pointer to scene
    jointObjects::Vector{Object3D}  # References to Object3Ds that have a joint
    jointStartIndex::Vector{Int}    # Start index of joint in qdd
    jointNdof::Vector{Int}          # Number-of-degrees-of-freedom of joint
    zStartIndex::Int                # eventHandler.z[zStartIndex] is first index of crossing function
                                    # (or zero, if nableContactDetection=false)
    nz::Int                         # Number of used zero crossing functions
    residuals::Vector{FloatType}    # Residuals - length(residuals) = nqdd
    cache_h::Vector{FloatType}      # Cached vector: = h(q,qd,gravity,contact-forces)
end

function multibodyResiduals(id::Int, _leq_mode, simulationModel::ModiaLang.SimulationModel{FloatType}, time, qdd, args...)::Vector{FloatType} where {FloatType}
     TimerOutputs.@timeit simulationModel.timer "Modia3D" begin
        separateObjects = simulationModel.separateObjects  # is emptied for every new simulate! call
        parameters      = simulationModel.evaluatedParameters
        if haskey(separateObjects, id)
            first = false
            mbs::MultibodyData{FloatType} = separateObjects[id]
            @assert(length(qdd) == mbs.nqdd)
            world           = mbs.world
            scene           = mbs.scene
            jointObjects    = mbs.jointObjects
            jointStartIndex = mbs.jointStartIndex
            jointNdof       = mbs.jointNdof
            residuals       = mbs.residuals
            cache_h         = mbs.cache_h
            zStartIndex     = mbs.zStartIndex
            nz              = mbs.nz
        else
            first = true

            # Search in parameters and retrieve the name of the object with the required id
            # Instantiate the Modia3D system
            #mbsPar = getIdParameter(parameters, id)
            #mbsObj = instantiateDependentObjects(simulationModel.modelModule, mbsPar)
            (world, jointObjects) = checkMultibodySystemAndGetWorldAndJoints(simulationModel, id)

            # Construct startIndex vector and number of degrees of freedom per joint
            nJoints         = length(jointObjects)
            jointStartIndex = fill(0,nJoints)
            jointNdof       = fill(0,nJoints)
            j = 1
            for (i, jointObject) in enumerate(jointObjects)
                jointStartIndex[i] = j
                jointNdof[i]       = jointObject.joint.ndof
                j += jointNdof[i]
            end
            nqdd = j-1
            if length(qdd) != nqdd
                error("Bug in Modia3D/src/Composition/dynamics.jl or in the generated code:\n",
                    "   length(qdd) = ", length(qdd), " is not identical to nqdd = ", nqdd)
            end

            # Construct MultibodyData
            scene = initAnalysis2!(world)
            residuals = zeros(FloatType,nqdd)
            cache_h   = zeros(FloatType,nqdd)
            if scene.options.enableContactDetection
                nz = 2
                zStartIndex = ModiaLang.addZeroCrossings(simulationModel, nz)
            else
                nz = 0
                zStartIndex = 0
            end
            mbs = MultibodyData{FloatType}(length(qdd), world, scene, jointObjects, jointStartIndex,
                                        jointNdof, zStartIndex, nz, residuals, cache_h)
            separateObjects[id] = mbs

            if simulationModel.options.log && scene.visualize
                println("          Modia3D: Number of visual shapes = ", length(scene.allVisuElements))
            end

            # Print
            if false
                printScene(scene)
            end
        end

        # Copy generalized position, velocity, acceleration, force values in to the joints
        @assert(length(args) == length(jointObjects))
        setJointVariables!(scene, jointObjects, args, qdd, jointStartIndex, jointNdof)::Nothing

        # Compute residuals
        leq_mode = isnothing(_leq_mode) ? -1 : _leq_mode.mode

        # Method with improved efficiency
        multibodyResiduals3!(simulationModel, scene, world, time, simulationModel.storeResult, first,
                            ModiaLang.isTerminal(simulationModel) && leq_mode == -1, leq_mode)

        # Copy the joint residuals in to the residuals vector
        if leq_mode == 0
            getJointResiduals_for_leq_mode_0!(scene, jointObjects, residuals, jointStartIndex, jointNdof, cache_h)
        elseif leq_mode > 0
            getJointResiduals_for_leq_mode_pos!(scene, jointObjects, residuals, jointStartIndex, jointNdof, cache_h)
        else
            residuals .= convert(FloatType, 0)
        end
    end
    return residuals
end

#= Efficiency improvement, by computing only terms that are needed in the respective phase
   of ModiaBase/src/EquationAndStateInfo.jl - Base.iterate(iterator::LinearEquationsIterator, ...)

 nf: Number of degrees of freedom
  q: Generalized joint coordinates of the tree joints
  u: Generalized joint forces of the tree joints
res: Residuals
e_i: i-th unit vector of length(q)

qd  = der(q)
qdd = der(qd)
res = M(q)*qdd + h(q,qd,gravity,contact forces) - u(q,qd,qdd)

multibodyResidual(..., qdd, ...) is called inside the LinearEquationsIterator to solve a linear equation system:

    res := A*x - b

that is called in the following way:

    leq_mode = 0        : Return "res := -b"        or "res :=            h(q,qd,gravity,contact forces) - u(q,qd,qdd=0)"
    leq_mode = i (1..nf): Return "res := A*e_i - b" or "res := M(q)*e_i + h(q,qd,gravity,contact forces) - u(q,qd,qdd=e_i)"
    leq_mode = -1       : Return "res := A*x - b"   or "res := 0.0"   # since residual is not used, it need not be computed.

For Modia3D:
    cache_h = h(q,qd,gravity,contact forces)

    leq_mode =  0: Compute only terms that depend on q, qd, gravity and contact forces.
                   All position and velocity terms, such as Rrel and Rabs are computed
                   and stored in Object3D.
                   res := h(q,qd,gravity) - u(q,qd,qdd)   # qdd = 0
                   cache_h := res + u
                   return res

    leq_mode =  i: res := M(q)*qdd - u(q,qd,qdd)         # qdd = e_i
                   return res := res + cache_h

    leq_mode = -1: res := M(q)*qdd - u(q,qd,qdd)         # qdd = A\b
                   return res := res + cache_h
=#

function multibodyResiduals3!(sim, scene, world, time, storeResults, isFirstInitial, isTerminal, leq_mode)
    tree            = scene.treeForComputation
    visualize       = scene.visualize   # && sim.model.visualiz
    exportAnimation = scene.exportAnimation

    # Handle initialization and termination of model
    # if Modia.isFirstInitialOfAllSegments(sim)
    if isFirstInitial
        if visualize
            TimerOutputs.@timeit sim.timer "Modia3D_0 initializeVisualization" Modia3D.Composition.initializeVisualization(Modia3D.renderer[1], scene.allVisuElements)
        end
    end

    if isTerminal  #if Modia.isTerminalOfAllSegments(sim)
        TimerOutputs.@timeit sim.timer "Modia3D_4 isTerminal" begin
            if exportAnimation
                Modia3D.exportAnimation(scene)
            end
            if visualize
                closeVisualization(Modia3D.renderer[1])
            end
        end
        #if scene.collide
        #    closeContactDetection!(m.assembly)
        #end
        # Do not return, since otherwise the linear system of equations cannot be solved -> error
        #return nothing
    end

    # Initialize force/torque of world-frame
    world.f = Modia3D.ZeroVector3D
    world.t = Modia3D.ZeroVector3D

    compute_vis = storeResults && !isTerminal && (visualize || exportAnimation)

    # Computation depending on leq_mode (the mode of the LinearEquationsIterator)
    if leq_mode == 0
        TimerOutputs.@timeit sim.timer "Modia3D_1" begin
            # Compute only terms that depend on q, qd, gravity and contact forces.
            # All position and velocity terms, such as Rrel and Rabs are computed
            # and stored in Object3D.
            #   res := h(q,qd,gravity) - u(q,qd,qdd)   # qdd = 0
            #   cache_h := res + u
            #   return res

            # Compute kinematics
            TimerOutputs.@timeit sim.timer "Modia3D_1 computeKinematics!" computeKinematics!(scene, tree, time)

            # Compute mass/inertia forces in a forward recursion and initializes forces/torques
            for obj in tree
                if compute_vis && length( obj.visualizationFrame ) == 1
                    obj.visualizationFrame[1].r_abs = obj.r_abs
                    obj.visualizationFrame[1].R_abs = obj.R_abs
                end

                if obj.hasMass
                    # Compute inertia forces / torques
                    w     = obj.w
                    grav  = gravityAcceleration(scene.options.gravityField, obj.r_abs)
                    obj.f = -obj.m*( obj.R_abs*(obj.a0 - grav) + cross(obj.z, obj.r_CM) + cross(w, cross(w, obj.r_CM)))
                    obj.t = -(obj.I_CM*obj.z + cross(w, obj.I_CM*w)) + cross(obj.r_CM, obj.f)
                else
                    obj.f = Modia3D.ZeroVector3D
                    obj.t = Modia3D.ZeroVector3D
                end
            end # end forward recursion

            # Compute contact forces/torques
            if scene.collide
                computeContactForcesAndTorques(sim, scene, world, time, nothing)
            end

            # Backward recursion (Compute forces/torques and residues)
            TimerOutputs.@timeit sim.timer "Modia3D_1 computeForcesAndResiduals" computeForcesTorquesAndResiduals!(scene, tree, time)
        end
        return nothing

    elseif leq_mode > 0
        TimerOutputs.@timeit sim.timer "Modia3D_2" begin
            # Compute   res := M(q)*e_i - u(q,qd,qdd=e_i)
            #           res := res + cache_h

            # Compute kinematics
            TimerOutputs.@timeit sim.timer "Modia3D_2 computeKinematics!" computeKinematics_for_leq_mode_pos!(scene, tree, time)

            for obj in tree
                if obj.hasMass
                    # Compute inertia forces / torques
                    obj.f = -obj.m*( obj.R_abs*obj.a0 + cross(obj.z, obj.r_CM) )
                    obj.t = -obj.I_CM*obj.z + cross(obj.r_CM, obj.f)
                else
                    obj.f = Modia3D.ZeroVector3D
                    obj.t = Modia3D.ZeroVector3D
                end
            end # end forward recursion

            # Compute forces/torques and residues in a backward recursion
            TimerOutputs.@timeit sim.timer "Modia3D_2 computeForcesAndResiduals" computeForcesTorquesAndResiduals!(scene, tree,time)
        end
        return nothing

    elseif leq_mode == -1
        TimerOutputs.@timeit sim.timer "Modia3D_3" begin
            if compute_vis  # Visualize at a communication point
                # Compute positions of frames that are only used for visualization
                TimerOutputs.@timeit sim.timer "Modia3D_3 visualize!" begin
                    if scene.options.useOptimizedStructure
                        for obj in scene.updateVisuElements
                            parent = obj.parent
                            obj.r_abs = obj.r_rel ≡ Modia3D.ZeroVector3D ? parent.r_abs : parent.r_abs + parent.R_abs'*obj.r_rel
                            obj.R_abs = obj.R_rel ≡ Modia3D.NullRotation ? parent.R_abs : obj.R_rel*parent.R_abs
                            # is executed only if an internal Object3D called
                            if length( obj.visualizationFrame ) == 1
                                obj.visualizationFrame[1].r_abs = obj.r_abs
                                obj.visualizationFrame[1].R_abs = obj.R_abs
                            end
                            for mesh in obj.fileMeshConvexShapes
                                mesh.r_abs = obj.r_abs
                                mesh.R_abs = obj.R_abs
                            end
                        end
                    end
                    if visualize
                        Modia3D.visualize!(Modia3D.renderer[1], time)
                    end
                end
                if exportAnimation
                    TimerOutputs.@timeit sim.timer "Modia3D_3 exportAnimation" begin
                        objectData = []
                        for obj in scene.allVisuElements
                            pos = obj.r_abs
                            ori = Modia3D.from_R(obj.R_abs)
                            dat = animationData(pos, ori)
                            push!(objectData, dat)
                        end
                        aniStep = animationStep(time, objectData)
                        push!(scene.animation, aniStep)
                    end
                end
            end
        end

    else
        error("Bug in dynamics.jl: leq_mode = $leq_mode (must be >= -1)")
    end

    # Copy variables to residues
    #Modia3D.copy_variables_to_residue!(var, _x, _derx, _r)
    #end # from open
    return nothing
end