text:
module Klio

using HTTP
using Sockets
using TimeZones
using Dates

include("json_handler.jl")
include("mattermost_types.jl")

include("cmd_calc.jl")
include("cmd_choose.jl")
include("cmd_julia.jl")
include("cmd_time.jl")
include("cmd_expl.jl")

using .Mattermost
using .JSON

mutable struct Settings
    server_host::Sockets.IPAddr # HTTP server's IP Address to listen on
    server_port::UInt16 # HTTP server's TCP Port to listen on
    server_verbose::Bool # HTTP server's verbosity
    expl_sqlite_file::String # SQLite path and file
    expl_time_zone::TimeZone # time zone used for DateTime formatting in !expl responses
    expl_datetime_format::DateFormat # format for DateTime in !expl

    Settings(;
        server_host = Sockets.localhost,
        server_port = 8000,
        server_verbose = true,
        expl_sqlite_file = "/tmp/expl.sqlite",
        expl_time_zone = tz"Europe/Berlin",
        expl_datetime_format = dateformat"dd.mm.YYYY HH:MM") =
            new(server_host, server_port, server_verbose, expl_sqlite_file, expl_time_zone, expl_datetime_format)
end

settings = Settings()

function run()
    klioRouter = HTTP.Router()

    wrap(h) = JSONHandler(OutgoingWebhookRequest, h)

    HTTP.@register(klioRouter, "POST", "/calc", wrap(Calc.calc))
    HTTP.@register(klioRouter, "POST", "/choose", wrap(Choose.choose))
    HTTP.@register(klioRouter, "POST", "/julia", wrap(Julia.julia))
    HTTP.@register(klioRouter, "POST", "/time", wrap(Time.time))

    HTTP.@register(klioRouter, "POST", "/add", wrap(Expl.add))
    HTTP.@register(klioRouter, "POST", "/expl", wrap(Expl.expl))
    HTTP.@register(klioRouter, "POST", "/del", wrap(Expl.del))
    HTTP.@register(klioRouter, "POST", "/find", wrap(Expl.find))
    HTTP.@register(klioRouter, "POST", "/topexpl", wrap(Expl.topexpl))

    HTTP.@register(klioRouter, "GET", "/expl", Expl.www_expl)
    HTTP.@register(klioRouter, "GET", "/find", Expl.www_find)

    HTTP.serve(klioRouter, settings.server_host, settings.server_port, verbose = settings.server_verbose)
end

end


text:
"""
    AbstractPredictiveAsymmetryTest

An abstract type for predictive asymmetry causality tests.

Concrete subtypes must implement the following fields:

- **`predictive_test`**. A predictive test to use for computing the
    asymmetries.
"""
abstract type AbstractPredictiveAsymmetryTest <: CausalityTest end

"""
    get_ηs

Get the prediction lags for a predictive asymmetry test.
"""
get_ηs(x::CT) where {CT <: AbstractPredictiveAsymmetryTest} = x.predictive_test.ηs

export AbstractPredictiveAsymmetryTest, get_ηs

ext: jl


text:
using ClimateMachine.Mesh.Elements
using ClimateMachine.Mesh.Metrics
using Test

const VGEO2D = (x1 = 1, x2 = 2, J = 3, ξ1x1 = 4, ξ2x1 = 5, ξ1x2 = 6, ξ2x2 = 7)
const SGEO2D = (sJ = 1, n1 = 2, n2 = 3)

const VGEO3D = (
    x1 = 1,
    x2 = 2,
    x3 = 3,
    J = 4,
    ξ1x1 = 5,
    ξ2x1 = 6,
    ξ3x1 = 7,
    ξ1x2 = 8,
    ξ2x2 = 9,
    ξ3x2 = 10,
    ξ1x3 = 11,
    ξ2x3 = 12,
    ξ3x3 = 13,
)
const SGEO3D = (sJ = 1, n1 = 2, n2 = 3, n3 = 4)

@testset "1-D Metric terms" begin
    for T in (Float32, Float64)
        #{{{
        let
            N = 4

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 1
            nfaces = 4
            e2c = Array{T, 3}(undef, 1, 2, 2)
            e2c[:, :, 1] = [-1 0]
            e2c[:, :, 2] = [0 10]
            nelem = size(e2c, 3)

            (x1,) = Metrics.creategrid1d(e2c, r)
            @test x1[:, 1] ≈ (r .- 1) / 2
            @test x1[:, 2] ≈ 5 * (r .+ 1)

            metric = Metrics.computemetric(x1, D)
            @test metric.J[:, 1] ≈ ones(T, Nq) / 2
            @test metric.J[:, 2] ≈ 5 * ones(T, Nq)
            @test metric.ξ1x1[:, 1] ≈ 2 * ones(T, Nq)
            @test metric.ξ1x1[:, 2] ≈ ones(T, Nq) / 5
            @test metric.n1[1, 1, :] ≈ -ones(T, nelem)
            @test metric.n1[1, 2, :] ≈ ones(T, nelem)
            @test metric.sJ[1, 1, :] ≈ ones(T, nelem)
            @test metric.sJ[1, 2, :] ≈ ones(T, nelem)
        end
        #}}}
    end
end

@testset "2-D Metric terms" begin
    for T in (Float32, Float64)
        # linear and rotation test
        #{{{
        let
            N = 2

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 2
            nfaces = 4
            e2c = Array{T, 3}(undef, 2, 4, 4)
            e2c[:, :, 1] = [0 2 0 2; 0 0 2 2]
            e2c[:, :, 2] = [2 2 0 0; 0 2 0 2]
            e2c[:, :, 3] = [2 0 2 0; 2 2 0 0]
            e2c[:, :, 4] = [0 0 2 2; 2 0 2 0]
            nelem = size(e2c, 3)

            x_exact = Array{Int, 3}(undef, 3, 3, 4)
            x_exact[:, :, 1] = [0 0 0; 1 1 1; 2 2 2]
            x_exact[:, :, 2] = rotr90(x_exact[:, :, 1])
            x_exact[:, :, 3] = rotr90(x_exact[:, :, 2])
            x_exact[:, :, 4] = rotr90(x_exact[:, :, 3])

            y_exact = Array{Int, 3}(undef, 3, 3, 4)
            y_exact[:, :, 1] = [0 1 2; 0 1 2; 0 1 2]
            y_exact[:, :, 2] = rotr90(y_exact[:, :, 1])
            y_exact[:, :, 3] = rotr90(y_exact[:, :, 2])
            y_exact[:, :, 4] = rotr90(y_exact[:, :, 3])

            J_exact = ones(Int, 3, 3, 4)

            ξ1x1_exact = zeros(Int, 3, 3, 4)
            ξ1x1_exact[:, :, 1] .= 1
            ξ1x1_exact[:, :, 3] .= -1

            ξ1x2_exact = zeros(Int, 3, 3, 4)
            ξ1x2_exact[:, :, 2] .= 1
            ξ1x2_exact[:, :, 4] .= -1

            ξ2x1_exact = zeros(Int, 3, 3, 4)
            ξ2x1_exact[:, :, 2] .= -1
            ξ2x1_exact[:, :, 4] .= 1

            ξ2x2_exact = zeros(Int, 3, 3, 4)
            ξ2x2_exact[:, :, 1] .= 1
            ξ2x2_exact[:, :, 3] .= -1

            sJ_exact = ones(Int, Nq, nfaces, nelem)

            nx_exact = zeros(Int, Nq, nfaces, nelem)
            nx_exact[:, 1, 1] .= -1
            nx_exact[:, 2, 1] .= 1
            nx_exact[:, 3, 2] .= 1
            nx_exact[:, 4, 2] .= -1
            nx_exact[:, 1, 3] .= 1
            nx_exact[:, 2, 3] .= -1
            nx_exact[:, 3, 4] .= -1
            nx_exact[:, 4, 4] .= 1

            ny_exact = zeros(Int, Nq, nfaces, nelem)
            ny_exact[:, 3, 1] .= -1
            ny_exact[:, 4, 1] .= 1
            ny_exact[:, 1, 2] .= -1
            ny_exact[:, 2, 2] .= 1
            ny_exact[:, 3, 3] .= 1
            ny_exact[:, 4, 3] .= -1
            ny_exact[:, 1, 4] .= 1
            ny_exact[:, 2, 4] .= -1

            vgeo = Array{T, 4}(undef, Nq, Nq, length(VGEO2D), nelem)
            sgeo = Array{T, 4}(undef, Nq, nfaces, length(SGEO2D), nelem)
            Metrics.creategrid!(ntuple(j -> (@view vgeo[:, :, j, :]), d)..., e2c, r)
            Metrics.computemetric!(
                ntuple(j -> (@view vgeo[:, :, j, :]), length(VGEO2D))...,
                ntuple(j -> (@view sgeo[:, :, j, :]), length(SGEO2D))...,
                D,
            )

            @test (@view vgeo[:, :, VGEO2D.x1, :]) ≈ x_exact
            @test (@view vgeo[:, :, VGEO2D.x2, :]) ≈ y_exact
            @test (@view vgeo[:, :, VGEO2D.J, :]) ≈ J_exact
            @test (@view vgeo[:, :, VGEO2D.ξ1x1, :]) ≈ ξ1x1_exact
            @test (@view vgeo[:, :, VGEO2D.ξ1x2, :]) ≈ ξ1x2_exact
            @test (@view vgeo[:, :, VGEO2D.ξ2x1, :]) ≈ ξ2x1_exact
            @test (@view vgeo[:, :, VGEO2D.ξ2x2, :]) ≈ ξ2x2_exact
            @test (@view sgeo[:, :, SGEO2D.sJ, :]) ≈ sJ_exact
            @test (@view sgeo[:, :, SGEO2D.n1, :]) ≈ nx_exact
            @test (@view sgeo[:, :, SGEO2D.n2, :]) ≈ ny_exact

            nothing
        end
        #}}}

        # Polynomial 2-D test
        #{{{
        let
            N = 4

            f(r, s) = (
                9 .* r - (1 .+ r) .* s .^ 2 +
                (r .- 1) .^ 2 .* (1 .- s .^ 2 .+ s .^ 3),
                10 .* s .+ r .^ 4 .* (1 .- s) .+ r .^ 2 .* s .* (1 .+ s),
            )
            fx1ξ1(r, s) =
                7 .+ s .^ 2 .- 2 .* s .^ 3 .+ 2 .* r .* (1 .- s .^ 2 .+ s .^ 3)
            fx1ξ2(r, s) =
                -2 .* (1 .+ r) .* s .+ (-1 .+ r) .^ 2 .* s .* (-2 .+ 3 .* s)
            fx2ξ1(r, s) = -4 .* r .^ 3 .* (-1 .+ s) .+ 2 .* r .* s .* (1 .+ s)
            fx2ξ2(r, s) = 10 .- r .^ 4 .+ r .^ 2 .* (1 .+ 2 .* s)

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 2
            nfaces = 4
            e2c = Array{T, 3}(undef, 2, 4, 1)
            e2c[:, :, 1] = [-1 1 -1 1; -1 -1 1 1]
            nelem = size(e2c, 3)

            vgeo = Array{T, 4}(undef, Nq, Nq, length(VGEO2D), nelem)
            sgeo = Array{T, 4}(undef, Nq, nfaces, length(SGEO2D), nelem)

            Metrics.creategrid!(ntuple(j -> (@view vgeo[:, :, j, :]), d)..., e2c, r)
            x1 = @view vgeo[:, :, VGEO2D.x1, :]
            x2 = @view vgeo[:, :, VGEO2D.x2, :]

            (x1ξ1, x1ξ2, x2ξ1, x2ξ2) =
                (fx1ξ1(x1, x2), fx1ξ2(x1, x2), fx2ξ1(x1, x2), fx2ξ2(x1, x2))
            J = x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1
            foreach(j -> (x1[j], x2[j]) = f(x1[j], x2[j]), 1:length(x1))

            Metrics.computemetric!(
                ntuple(j -> (@view vgeo[:, :, j, :]), length(VGEO2D))...,
                ntuple(j -> (@view sgeo[:, :, j, :]), length(SGEO2D))...,
                D,
            )
            @test J ≈ (@view vgeo[:, :, VGEO2D.J, :])
            @test (@view vgeo[:, :, VGEO2D.ξ1x1, :]) ≈ x2ξ2 ./ J
            @test (@view vgeo[:, :, VGEO2D.ξ2x1, :]) ≈ -x2ξ1 ./ J
            @test (@view vgeo[:, :, VGEO2D.ξ1x2, :]) ≈ -x1ξ2 ./ J
            @test (@view vgeo[:, :, VGEO2D.ξ2x2, :]) ≈ x1ξ1 ./ J

            # check the normals?
            n1 = @view sgeo[:, :, SGEO2D.n1, :]
            n2 = @view sgeo[:, :, SGEO2D.n2, :]
            sJ = @view sgeo[:, :, SGEO2D.sJ, :]
            @test hypot.(n1, n2) ≈ ones(T, size(n1))
            @test sJ[:, 1, :] .* n1[:, 1, :] ≈ -x2ξ2[1, :, :]
            @test sJ[:, 1, :] .* n2[:, 1, :] ≈ x1ξ2[1, :, :]
            @test sJ[:, 2, :] .* n1[:, 2, :] ≈ x2ξ2[Nq, :, :]
            @test sJ[:, 2, :] .* n2[:, 2, :] ≈ -x1ξ2[Nq, :, :]
            @test sJ[:, 3, :] .* n1[:, 3, :] ≈ x2ξ1[:, 1, :]
            @test sJ[:, 3, :] .* n2[:, 3, :] ≈ -x1ξ1[:, 1, :]
            @test sJ[:, 4, :] .* n1[:, 4, :] ≈ -x2ξ1[:, Nq, :]
            @test sJ[:, 4, :] .* n2[:, 4, :] ≈ x1ξ1[:, Nq, :]
        end
        #}}}

        #{{{
        let
            N = 4

            f(r, s) = (
                9 .* r - (1 .+ r) .* s .^ 2 +
                (r .- 1) .^ 2 .* (1 .- s .^ 2 .+ s .^ 3),
                10 .* s .+ r .^ 4 .* (1 .- s) .+ r .^ 2 .* s .* (1 .+ s),
            )
            fx1ξ1(r, s) =
                7 .+ s .^ 2 .- 2 .* s .^ 3 .+ 2 .* r .* (1 .- s .^ 2 .+ s .^ 3)
            fx1ξ2(r, s) =
                -2 .* (1 .+ r) .* s .+ (-1 .+ r) .^ 2 .* s .* (-2 .+ 3 .* s)
            fx2ξ1(r, s) = -4 .* r .^ 3 .* (-1 .+ s) .+ 2 .* r .* s .* (1 .+ s)
            fx2ξ2(r, s) = 10 .- r .^ 4 .+ r .^ 2 .* (1 .+ 2 .* s)

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 2
            nfaces = 4
            e2c = Array{T, 3}(undef, 2, 4, 1)
            e2c[:, :, 1] = [-1 1 -1 1; -1 -1 1 1]
            nelem = size(e2c, 3)

            (x1, x2) = Metrics.creategrid2d(e2c, r)

            (x1ξ1, x1ξ2, x2ξ1, x2ξ2) =
                (fx1ξ1(x1, x2), fx1ξ2(x1, x2), fx2ξ1(x1, x2), fx2ξ2(x1, x2))
            J = x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1
            foreach(j -> (x1[j], x2[j]) = f(x1[j], x2[j]), 1:length(x1))

            metric = Metrics.computemetric(x1, x2, D)
            @test J ≈ metric.J
            @test metric.ξ1x1 ≈ x2ξ2 ./ J
            @test metric.ξ2x1 ≈ -x2ξ1 ./ J
            @test metric.ξ1x2 ≈ -x1ξ2 ./ J
            @test metric.ξ2x2 ≈ x1ξ1 ./ J

            # check the normals?
            n1 = metric.n1
            n2 = metric.n2
            sJ = metric.sJ
            @test hypot.(n1, n2) ≈ ones(T, size(n1))
            @test sJ[:, 1, :] .* n1[:, 1, :] ≈ -x2ξ2[1, :, :]
            @test sJ[:, 1, :] .* n2[:, 1, :] ≈ x1ξ2[1, :, :]
            @test sJ[:, 2, :] .* n1[:, 2, :] ≈ x2ξ2[Nq, :, :]
            @test sJ[:, 2, :] .* n2[:, 2, :] ≈ -x1ξ2[Nq, :, :]
            @test sJ[:, 3, :] .* n1[:, 3, :] ≈ x2ξ1[:, 1, :]
            @test sJ[:, 3, :] .* n2[:, 3, :] ≈ -x1ξ1[:, 1, :]
            @test sJ[:, 4, :] .* n1[:, 4, :] ≈ -x2ξ1[:, Nq, :]
            @test sJ[:, 4, :] .* n2[:, 4, :] ≈ x1ξ1[:, Nq, :]
        end
        #}}}
    end

    # Constant preserving test
    #{{{
    let
        N = 4
        T = Float64

        f(r, s) = (
            9 .* r - (1 .+ r) .* s .^ 2 +
            (r .- 1) .^ 2 .* (1 .- s .^ 2 .+ s .^ 3),
            10 .* s .+ r .^ 4 .* (1 .- s) .+ r .^ 2 .* s .* (1 .+ s),
        )
        fx1ξ1(r, s) =
            7 .+ s .^ 2 .- 2 .* s .^ 3 .+ 2 .* r .* (1 .- s .^ 2 .+ s .^ 3)
        fx1ξ2(r, s) =
            -2 .* (1 .+ r) .* s .+ (-1 .+ r) .^ 2 .* s .* (-2 .+ 3 .* s)
        fx2ξ1(r, s) = -4 .* r .^ 3 .* (-1 .+ s) .+ 2 .* r .* s .* (1 .+ s)
        fx2ξ2(r, s) = 10 .- r .^ 4 .+ r .^ 2 .* (1 .+ 2 .* s)

        r, w = Elements.lglpoints(T, N)
        D = Elements.spectralderivative(r)
        Nq = N + 1

        d = 2
        nfaces = 4
        e2c = Array{T, 3}(undef, 2, 4, 1)
        e2c[:, :, 1] = [-1 1 -1 1; -1 -1 1 1]
        nelem = size(e2c, 3)

        vgeo = Array{T, 4}(undef, Nq, Nq, length(VGEO2D), nelem)
        sgeo = Array{T, 4}(undef, Nq, nfaces, length(SGEO2D), nelem)

        Metrics.creategrid!(ntuple(j -> (@view vgeo[:, :, j, :]), d)..., e2c, r)
        x1 = @view vgeo[:, :, VGEO2D.x1, :]
        x2 = @view vgeo[:, :, VGEO2D.x2, :]

        foreach(j -> (x1[j], x2[j]) = f(x1[j], x2[j]), 1:length(x1))

        Metrics.computemetric!(
            ntuple(j -> (@view vgeo[:, :, j, :]), length(VGEO2D))...,
            ntuple(j -> (@view sgeo[:, :, j, :]), length(SGEO2D))...,
            D,
        )

        (Cx1, Cx2) = (zeros(T, Nq, Nq), zeros(T, Nq, Nq))

        J = @view vgeo[:, :, VGEO2D.J, :]
        ξ1x1 = @view vgeo[:, :, VGEO2D.ξ1x1, :]
        ξ1x2 = @view vgeo[:, :, VGEO2D.ξ1x2, :]
        ξ2x1 = @view vgeo[:, :, VGEO2D.ξ2x1, :]
        ξ2x2 = @view vgeo[:, :, VGEO2D.ξ2x2, :]

        e = 1
        for n in 1:Nq
            Cx1[:, n] += D * (J[:, n, e] .* ξ1x1[:, n, e])
            Cx1[n, :] += D * (J[n, :, e] .* ξ2x1[n, :, e])

            Cx2[:, n] += D * (J[:, n, e] .* ξ1x2[:, n, e])
            Cx2[n, :] += D * (J[n, :, e] .* ξ2x2[n, :, e])
        end
        @test maximum(abs.(Cx1)) ≤ 1000 * eps(T)
        @test maximum(abs.(Cx2)) ≤ 1000 * eps(T)
    end
    #}}}
end

@testset "3-D Metric terms" begin
    # linear test
    #{{{
    for T in (Float32, Float64)
        let
            N = 2

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 3
            nfaces = 6
            e2c = Array{T, 3}(undef, d, 8, 2)
            e2c[:, :, 1] = [
                0 2 0 2 0 2 0 2
                0 0 2 2 0 0 2 2
                0 0 0 0 2 2 2 2
            ]
            e2c[:, :, 2] = [
                2 2 0 0 2 2 0 0
                0 2 0 2 0 2 0 2
                0 0 0 0 2 2 2 2
            ]

            nelem = size(e2c, 3)

            x_exact = Array{Int, 4}(undef, 3, 3, 3, nelem)
            x_exact[1, :, :, 1] .= 0
            x_exact[2, :, :, 1] .= 1
            x_exact[3, :, :, 1] .= 2
            x_exact[:, 1, :, 2] .= 2
            x_exact[:, 2, :, 2] .= 1
            x_exact[:, 3, :, 2] .= 0

            ξ1x1_exact = zeros(Int, 3, 3, 3, nelem)
            ξ1x1_exact[:, :, :, 1] .= 1

            ξ1x2_exact = zeros(Int, 3, 3, 3, nelem)
            ξ1x2_exact[:, :, :, 2] .= 1

            ξ2x1_exact = zeros(Int, 3, 3, 3, nelem)
            ξ2x1_exact[:, :, :, 2] .= -1

            ξ2x2_exact = zeros(Int, 3, 3, 3, nelem)
            ξ2x2_exact[:, :, :, 1] .= 1

            ξ3x3_exact = ones(Int, 3, 3, 3, nelem)

            y_exact = Array{Int, 4}(undef, 3, 3, 3, nelem)
            y_exact[:, 1, :, 1] .= 0
            y_exact[:, 2, :, 1] .= 1
            y_exact[:, 3, :, 1] .= 2
            y_exact[1, :, :, 2] .= 0
            y_exact[2, :, :, 2] .= 1
            y_exact[3, :, :, 2] .= 2

            z_exact = Array{Int, 4}(undef, 3, 3, 3, nelem)
            z_exact[:, :, 1, 1:2] .= 0
            z_exact[:, :, 2, 1:2] .= 1
            z_exact[:, :, 3, 1:2] .= 2

            J_exact = ones(Int, 3, 3, 3, nelem)

            sJ_exact = ones(Int, Nq, Nq, nfaces, nelem)

            nx_exact = zeros(Int, Nq, Nq, nfaces, nelem)
            nx_exact[:, :, 1, 1] .= -1
            nx_exact[:, :, 2, 1] .= 1
            nx_exact[:, :, 3, 2] .= 1
            nx_exact[:, :, 4, 2] .= -1

            ny_exact = zeros(Int, Nq, Nq, nfaces, nelem)
            ny_exact[:, :, 3, 1] .= -1
            ny_exact[:, :, 4, 1] .= 1
            ny_exact[:, :, 1, 2] .= -1
            ny_exact[:, :, 2, 2] .= 1

            nz_exact = zeros(Int, Nq, Nq, nfaces, nelem)
            nz_exact[:, :, 5, 1:2] .= -1
            nz_exact[:, :, 6, 1:2] .= 1

            vgeo = Array{T, 5}(undef, Nq, Nq, Nq, length(VGEO3D), nelem)
            sgeo = Array{T, 4}(undef, Nq^2, nfaces, length(SGEO3D), nelem)
            Metrics.creategrid!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), d)...,
                e2c,
                r,
            )
            Metrics.computemetric!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), length(VGEO3D))...,
                ntuple(j -> (@view sgeo[:, :, j, :]), length(SGEO3D))...,
                D,
            )
            sgeo = reshape(sgeo, Nq, Nq, nfaces, length(SGEO3D), nelem)

            @test (@view vgeo[:, :, :, VGEO3D.x1, :]) ≈ x_exact
            @test (@view vgeo[:, :, :, VGEO3D.x2, :]) ≈ y_exact
            @test (@view vgeo[:, :, :, VGEO3D.x3, :]) ≈ z_exact
            @test (@view vgeo[:, :, :, VGEO3D.J, :]) ≈ J_exact
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x1, :]) ≈ ξ1x1_exact
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x2, :]) ≈ ξ1x2_exact
            @test maximum(abs.(@view vgeo[
                :,
                :,
                :,
                VGEO3D.ξ1x3,
                :,
            ])) ≤ 10 * eps(T)
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x1, :]) ≈ ξ2x1_exact
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x2, :]) ≈ ξ2x2_exact
            @test maximum(abs.(@view vgeo[
                :,
                :,
                :,
                VGEO3D.ξ2x3,
                :,
            ])) ≤ 10 * eps(T)
            @test maximum(abs.(@view vgeo[
                :,
                :,
                :,
                VGEO3D.ξ3x1,
                :,
            ])) ≤ 10 * eps(T)
            @test maximum(abs.(@view vgeo[
                :,
                :,
                :,
                VGEO3D.ξ3x2,
                :,
            ])) ≤ 10 * eps(T)
            @test (@view vgeo[:, :, :, VGEO3D.ξ3x3, :]) ≈ ξ3x3_exact
            @test (@view sgeo[:, :, :, SGEO3D.sJ, :]) ≈ sJ_exact
            @test (@view sgeo[:, :, :, SGEO3D.n1, :]) ≈ nx_exact
            @test (@view sgeo[:, :, :, SGEO3D.n2, :]) ≈ ny_exact
            @test (@view sgeo[:, :, :, SGEO3D.n3, :]) ≈ nz_exact
        end
    end
    #}}}

    # linear test with rotation
    #{{{
    for T in (Float32, Float64)
        θ1 = 2 * T(π) * T(0.9)
        θ2 = 2 * T(π) * T(-0.56)
        θ3 = 2 * T(π) * T(0.33)
        #=
        θ1 = 2 * T(π) * rand(T)
        θ2 = 2 * T(π) * rand(T)
        θ3 = 2 * T(π) * rand(T)
        =#
        #=
        θ1 = T(π) / 6
        θ2 = T(π) / 12
        θ3 = 4 * T(π) / 5
        =#
        Q = [cos(θ1) -sin(θ1) 0; sin(θ1) cos(θ1) 0; 0 0 1]
        Q *= [cos(θ2) 0 -sin(θ2); 0 1 0; sin(θ2) 0 cos(θ2)]
        Q *= [1 0 0; 0 cos(θ3) -sin(θ3); 0 sin(θ3) cos(θ3)]
        let
            N = 2

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 3
            nfaces = 6
            e2c = Array{T, 3}(undef, d, 8, 1)
            e2c[:, :, 1] = [
                0 2 0 2 0 2 0 2
                0 0 2 2 0 0 2 2
                0 0 0 0 2 2 2 2
            ]
            @views (x1, x2, x3) = (e2c[1, :, 1], e2c[2, :, 1], e2c[3, :, 1])
            for i in 1:length(x1)
                (x1[i], x2[i], x3[i]) = Q * [x1[i]; x2[i]; x3[i]]
            end

            nelem = size(e2c, 3)

            x1e = Array{T, 4}(undef, 3, 3, 3, nelem)
            x1e[1, :, :, 1] .= 0
            x1e[2, :, :, 1] .= 1
            x1e[3, :, :, 1] .= 2

            x2e = Array{T, 4}(undef, 3, 3, 3, nelem)
            x2e[:, 1, :, 1] .= 0
            x2e[:, 2, :, 1] .= 1
            x2e[:, 3, :, 1] .= 2

            x3e = Array{T, 4}(undef, 3, 3, 3, nelem)
            x3e[:, :, 1, 1] .= 0
            x3e[:, :, 2, 1] .= 1
            x3e[:, :, 3, 1] .= 2

            for i in 1:length(x1e)
                (x1e[i], x2e[i], x3e[i]) = Q * [x1e[i]; x2e[i]; x3e[i]]
            end

            Je = ones(Int, 3, 3, 3, nelem)

            # By construction
            # Q = [x1ξ1 x1ξ2 x2ξ2; x2ξ1 x2ξ2 x2ξ3; x3ξ1 x3ξ2 x3ξ3] = [ξ1x1 ξ2x1 ξ3x1; ξ1x2 ξ2x2 ξ3x2; ξ1x3 ξ2x3 ξ3x3]
            ξ1x1e = fill(Q[1, 1], 3, 3, 3, nelem)
            ξ1x2e = fill(Q[2, 1], 3, 3, 3, nelem)
            ξ1x3e = fill(Q[3, 1], 3, 3, 3, nelem)
            ξ2x1e = fill(Q[1, 2], 3, 3, 3, nelem)
            ξ2x2e = fill(Q[2, 2], 3, 3, 3, nelem)
            ξ2x3e = fill(Q[3, 2], 3, 3, 3, nelem)
            ξ3x1e = fill(Q[1, 3], 3, 3, 3, nelem)
            ξ3x2e = fill(Q[2, 3], 3, 3, 3, nelem)
            ξ3x3e = fill(Q[3, 3], 3, 3, 3, nelem)

            sJe = ones(Int, Nq, Nq, nfaces, nelem)

            n1e = zeros(T, Nq, Nq, nfaces, nelem)
            n2e = zeros(T, Nq, Nq, nfaces, nelem)
            n3e = zeros(T, Nq, Nq, nfaces, nelem)

            fill!(@view(n1e[:, :, 1, :]), -Q[1, 1])
            fill!(@view(n1e[:, :, 2, :]), Q[1, 1])
            fill!(@view(n1e[:, :, 3, :]), -Q[1, 2])
            fill!(@view(n1e[:, :, 4, :]), Q[1, 2])
            fill!(@view(n1e[:, :, 5, :]), -Q[1, 3])
            fill!(@view(n1e[:, :, 6, :]), Q[1, 3])
            fill!(@view(n2e[:, :, 1, :]), -Q[2, 1])
            fill!(@view(n2e[:, :, 2, :]), Q[2, 1])
            fill!(@view(n2e[:, :, 3, :]), -Q[2, 2])
            fill!(@view(n2e[:, :, 4, :]), Q[2, 2])
            fill!(@view(n2e[:, :, 5, :]), -Q[2, 3])
            fill!(@view(n2e[:, :, 6, :]), Q[2, 3])
            fill!(@view(n3e[:, :, 1, :]), -Q[3, 1])
            fill!(@view(n3e[:, :, 2, :]), Q[3, 1])
            fill!(@view(n3e[:, :, 3, :]), -Q[3, 2])
            fill!(@view(n3e[:, :, 4, :]), Q[3, 2])
            fill!(@view(n3e[:, :, 5, :]), -Q[3, 3])
            fill!(@view(n3e[:, :, 6, :]), Q[3, 3])

            vgeo = Array{T, 5}(undef, Nq, Nq, Nq, length(VGEO3D), nelem)
            sgeo = Array{T, 4}(undef, Nq^2, nfaces, length(SGEO3D), nelem)
            Metrics.creategrid!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), d)...,
                e2c,
                r,
            )
            Metrics.computemetric!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), length(VGEO3D))...,
                ntuple(j -> (@view sgeo[:, :, j, :]), length(SGEO3D))...,
                D,
            )
            sgeo = reshape(sgeo, Nq, Nq, nfaces, length(SGEO3D), nelem)

            @test (@view vgeo[:, :, :, VGEO3D.x1, :]) ≈ x1e
            @test (@view vgeo[:, :, :, VGEO3D.x2, :]) ≈ x2e
            @test (@view vgeo[:, :, :, VGEO3D.x3, :]) ≈ x3e
            @test (@view vgeo[:, :, :, VGEO3D.J, :]) ≈ Je
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x1, :]) ≈ ξ1x1e
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x2, :]) ≈ ξ1x2e
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x3, :]) ≈ ξ1x3e
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x1, :]) ≈ ξ2x1e
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x2, :]) ≈ ξ2x2e
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x3, :]) ≈ ξ2x3e
            @test (@view vgeo[:, :, :, VGEO3D.ξ3x1, :]) ≈ ξ3x1e
            @test (@view vgeo[:, :, :, VGEO3D.ξ3x2, :]) ≈ ξ3x2e
            @test (@view vgeo[:, :, :, VGEO3D.ξ3x3, :]) ≈ ξ3x3e
            @test (@view sgeo[:, :, :, SGEO3D.sJ, :]) ≈ sJe
            @test (@view sgeo[:, :, :, SGEO3D.n1, :]) ≈ n1e
            @test (@view sgeo[:, :, :, SGEO3D.n2, :]) ≈ n2e
            @test (@view sgeo[:, :, :, SGEO3D.n3, :]) ≈ n3e
        end
    end
    #}}}

    # Polynomial 3-D test
    #{{{
    for T in (Float32, Float64)
        f(r, s, t) = @.( (
            s + r * t - (r^2 * s^2 * t^2) / 4,
            t - ((r * s * t) / 2 + 1 / 2)^3 + 1,
            r + (r / 2 + 1 / 2)^6 * (s / 2 + 1 / 2)^6 * (t / 2 + 1 / 2)^6,
        ))

        fx1ξ1(r, s, t) = @.(t - (r * s^2 * t^2) / 2)
        fx1ξ2(r, s, t) = @.(1 - (r^2 * s * t^2) / 2)
        fx1ξ3(r, s, t) = @.(r - (r^2 * s^2 * t) / 2)
        fx2ξ1(r, s, t) = @.(-(3 * s * t * ((r * s * t) / 2 + 1 / 2)^2) / 2)
        fx2ξ2(r, s, t) = @.(-(3 * r * t * ((r * s * t) / 2 + 1 / 2)^2) / 2)
        fx2ξ3(r, s, t) = @.(1 - (3 * r * s * ((r * s * t) / 2 + 1 / 2)^2) / 2)
        fx3ξ1(r, s, t) = @.(
            3 * (r / 2 + 1 / 2)^5 * (s / 2 + 1 / 2)^6 * (t / 2 + 1 / 2)^6 + 1
        )
        fx3ξ2(r, s, t) =
            @.(3 * (r / 2 + 1 / 2)^6 * (s / 2 + 1 / 2)^5 * (t / 2 + 1 / 2)^6)
        fx3ξ3(r, s, t) =
            @.(3 * (r / 2 + 1 / 2)^6 * (s / 2 + 1 / 2)^6 * (t / 2 + 1 / 2)^5)

        let
            N = 9

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 3
            nfaces = 6
            e2c = Array{T, 3}(undef, d, 8, 1)
            e2c[:, :, 1] = [
                -1 1 -1 1 -1 1 -1 1
                -1 -1 1 1 -1 -1 1 1
                -1 -1 -1 -1 1 1 1 1
            ]

            nelem = size(e2c, 3)

            vgeo = Array{T, 5}(undef, Nq, Nq, Nq, length(VGEO3D), nelem)
            sgeo = Array{T, 4}(undef, Nq^2, nfaces, length(SGEO3D), nelem)
            Metrics.creategrid!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), d)...,
                e2c,
                r,
            )
            x1 = @view vgeo[:, :, :, VGEO3D.x1, :]
            x2 = @view vgeo[:, :, :, VGEO3D.x2, :]
            x3 = @view vgeo[:, :, :, VGEO3D.x3, :]

            (x1ξ1, x1ξ2, x1ξ3, x2ξ1, x2ξ2, x2ξ3, x3ξ1, x3ξ2, x3ξ3) = (
                fx1ξ1(x1, x2, x3),
                fx1ξ2(x1, x2, x3),
                fx1ξ3(x1, x2, x3),
                fx2ξ1(x1, x2, x3),
                fx2ξ2(x1, x2, x3),
                fx2ξ3(x1, x2, x3),
                fx3ξ1(x1, x2, x3),
                fx3ξ2(x1, x2, x3),
                fx3ξ3(x1, x2, x3),
            )
            J = (
                x1ξ1 .* (x2ξ2 .* x3ξ3 - x2ξ3 .* x3ξ2) +
                x2ξ1 .* (x3ξ2 .* x1ξ3 - x3ξ3 .* x1ξ2) +
                x3ξ1 .* (x1ξ2 .* x2ξ3 - x1ξ3 .* x2ξ2)
            )

            ξ1x1 = (x2ξ2 .* x3ξ3 - x2ξ3 .* x3ξ2) ./ J
            ξ1x2 = (x3ξ2 .* x1ξ3 - x3ξ3 .* x1ξ2) ./ J
            ξ1x3 = (x1ξ2 .* x2ξ3 - x1ξ3 .* x2ξ2) ./ J
            ξ2x1 = (x2ξ3 .* x3ξ1 - x2ξ1 .* x3ξ3) ./ J
            ξ2x2 = (x3ξ3 .* x1ξ1 - x3ξ1 .* x1ξ3) ./ J
            ξ2x3 = (x1ξ3 .* x2ξ1 - x1ξ1 .* x2ξ3) ./ J
            ξ3x1 = (x2ξ1 .* x3ξ2 - x2ξ2 .* x3ξ1) ./ J
            ξ3x2 = (x3ξ1 .* x1ξ2 - x3ξ2 .* x1ξ1) ./ J
            ξ3x3 = (x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1) ./ J

            foreach(
                j -> (x1[j], x2[j], x3[j]) = f(x1[j], x2[j], x3[j]),
                1:length(x1),
            )

            Metrics.computemetric!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), length(VGEO3D))...,
                ntuple(j -> (@view sgeo[:, :, j, :]), length(SGEO3D))...,
                D,
            )
            sgeo = reshape(sgeo, Nq, Nq, nfaces, length(SGEO3D), nelem)

            @test (@view vgeo[:, :, :, VGEO3D.J, :]) ≈ J
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x1, :]) ≈ ξ1x1
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x2, :]) ≈ ξ1x2
            @test (@view vgeo[:, :, :, VGEO3D.ξ1x3, :]) ≈ ξ1x3
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x1, :]) ≈ ξ2x1
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x2, :]) ≈ ξ2x2
            @test (@view vgeo[:, :, :, VGEO3D.ξ2x3, :]) ≈ ξ2x3
            @test (@view vgeo[:, :, :, VGEO3D.ξ3x1, :]) ≈ ξ3x1
            @test (@view vgeo[:, :, :, VGEO3D.ξ3x2, :]) ≈ ξ3x2
            @test (@view vgeo[:, :, :, VGEO3D.ξ3x3, :]) ≈ ξ3x3
            n1 = @view sgeo[:, :, :, SGEO3D.n1, :]
            n2 = @view sgeo[:, :, :, SGEO3D.n2, :]
            n3 = @view sgeo[:, :, :, SGEO3D.n3, :]
            sJ = @view sgeo[:, :, :, SGEO3D.sJ, :]
            @test hypot.(n1, n2, n3) ≈ ones(T, size(n1))
            @test (
                [
                    sJ[:, :, 1, :] .* n1[:, :, 1, :],
                    sJ[:, :, 1, :] .* n2[:, :, 1, :],
                    sJ[:, :, 1, :] .* n3[:, :, 1, :],
                ] ≈ [
                    -J[1, :, :, :] .* ξ1x1[1, :, :, :],
                    -J[1, :, :, :] .* ξ1x2[1, :, :, :],
                    -J[1, :, :, :] .* ξ1x3[1, :, :, :],
                ]
            )
            @test (
                [
                    sJ[:, :, 2, :] .* n1[:, :, 2, :],
                    sJ[:, :, 2, :] .* n2[:, :, 2, :],
                    sJ[:, :, 2, :] .* n3[:, :, 2, :],
                ] ≈ [
                    J[Nq, :, :, :] .* ξ1x1[Nq, :, :, :],
                    J[Nq, :, :, :] .* ξ1x2[Nq, :, :, :],
                    J[Nq, :, :, :] .* ξ1x3[Nq, :, :, :],
                ]
            )
            @test sJ[:, :, 3, :] .* n1[:, :, 3, :] ≈
                  -J[:, 1, :, :] .* ξ2x1[:, 1, :, :]
            @test sJ[:, :, 3, :] .* n2[:, :, 3, :] ≈
                  -J[:, 1, :, :] .* ξ2x2[:, 1, :, :]
            @test sJ[:, :, 3, :] .* n3[:, :, 3, :] ≈
                  -J[:, 1, :, :] .* ξ2x3[:, 1, :, :]
            @test sJ[:, :, 4, :] .* n1[:, :, 4, :] ≈
                  J[:, Nq, :, :] .* ξ2x1[:, Nq, :, :]
            @test sJ[:, :, 4, :] .* n2[:, :, 4, :] ≈
                  J[:, Nq, :, :] .* ξ2x2[:, Nq, :, :]
            @test sJ[:, :, 4, :] .* n3[:, :, 4, :] ≈
                  J[:, Nq, :, :] .* ξ2x3[:, Nq, :, :]
            @test sJ[:, :, 5, :] .* n1[:, :, 5, :] ≈
                  -J[:, :, 1, :] .* ξ3x1[:, :, 1, :]
            @test sJ[:, :, 5, :] .* n2[:, :, 5, :] ≈
                  -J[:, :, 1, :] .* ξ3x2[:, :, 1, :]
            @test sJ[:, :, 5, :] .* n3[:, :, 5, :] ≈
                  -J[:, :, 1, :] .* ξ3x3[:, :, 1, :]
            @test sJ[:, :, 6, :] .* n1[:, :, 6, :] ≈
                  J[:, :, Nq, :] .* ξ3x1[:, :, Nq, :]
            @test sJ[:, :, 6, :] .* n2[:, :, 6, :] ≈
                  J[:, :, Nq, :] .* ξ3x2[:, :, Nq, :]
            @test sJ[:, :, 6, :] .* n3[:, :, 6, :] ≈
                  J[:, :, Nq, :] .* ξ3x3[:, :, Nq, :]
        end
    end
    #}}}

    #{{{
    for T in (Float32, Float64)
        f(r, s, t) = @.( (
            s + r * t - (r^2 * s^2 * t^2) / 4,
            t - ((r * s * t) / 2 + 1 / 2)^3 + 1,
            r + (r / 2 + 1 / 2)^6 * (s / 2 + 1 / 2)^6 * (t / 2 + 1 / 2)^6,
        ))

        fx1ξ1(r, s, t) = @.(t - (r * s^2 * t^2) / 2)
        fx1ξ2(r, s, t) = @.(1 - (r^2 * s * t^2) / 2)
        fx1ξ3(r, s, t) = @.(r - (r^2 * s^2 * t) / 2)
        fx2ξ1(r, s, t) = @.(-(3 * s * t * ((r * s * t) / 2 + 1 / 2)^2) / 2)
        fx2ξ2(r, s, t) = @.(-(3 * r * t * ((r * s * t) / 2 + 1 / 2)^2) / 2)
        fx2ξ3(r, s, t) = @.(1 - (3 * r * s * ((r * s * t) / 2 + 1 / 2)^2) / 2)
        fx3ξ1(r, s, t) = @.(
            3 * (r / 2 + 1 / 2)^5 * (s / 2 + 1 / 2)^6 * (t / 2 + 1 / 2)^6 + 1
        )
        fx3ξ2(r, s, t) =
            @.(3 * (r / 2 + 1 / 2)^6 * (s / 2 + 1 / 2)^5 * (t / 2 + 1 / 2)^6)
        fx3ξ3(r, s, t) =
            @.(3 * (r / 2 + 1 / 2)^6 * (s / 2 + 1 / 2)^6 * (t / 2 + 1 / 2)^5)

        let
            N = 9

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 3
            nfaces = 6
            e2c = Array{T, 3}(undef, d, 8, 1)
            e2c[:, :, 1] = [
                -1 1 -1 1 -1 1 -1 1
                -1 -1 1 1 -1 -1 1 1
                -1 -1 -1 -1 1 1 1 1
            ]

            nelem = size(e2c, 3)

            (x1, x2, x3) = Metrics.creategrid3d(e2c, r)

            (x1ξ1, x1ξ2, x1ξ3, x2ξ1, x2ξ2, x2ξ3, x3ξ1, x3ξ2, x3ξ3) = (
                fx1ξ1(x1, x2, x3),
                fx1ξ2(x1, x2, x3),
                fx1ξ3(x1, x2, x3),
                fx2ξ1(x1, x2, x3),
                fx2ξ2(x1, x2, x3),
                fx2ξ3(x1, x2, x3),
                fx3ξ1(x1, x2, x3),
                fx3ξ2(x1, x2, x3),
                fx3ξ3(x1, x2, x3),
            )
            J = (
                x1ξ1 .* (x2ξ2 .* x3ξ3 - x2ξ3 .* x3ξ2) +
                x2ξ1 .* (x3ξ2 .* x1ξ3 - x3ξ3 .* x1ξ2) +
                x3ξ1 .* (x1ξ2 .* x2ξ3 - x1ξ3 .* x2ξ2)
            )

            ξ1x1 = (x2ξ2 .* x3ξ3 - x2ξ3 .* x3ξ2) ./ J
            ξ1x2 = (x3ξ2 .* x1ξ3 - x3ξ3 .* x1ξ2) ./ J
            ξ1x3 = (x1ξ2 .* x2ξ3 - x1ξ3 .* x2ξ2) ./ J
            ξ2x1 = (x2ξ3 .* x3ξ1 - x2ξ1 .* x3ξ3) ./ J
            ξ2x2 = (x3ξ3 .* x1ξ1 - x3ξ1 .* x1ξ3) ./ J
            ξ2x3 = (x1ξ3 .* x2ξ1 - x1ξ1 .* x2ξ3) ./ J
            ξ3x1 = (x2ξ1 .* x3ξ2 - x2ξ2 .* x3ξ1) ./ J
            ξ3x2 = (x3ξ1 .* x1ξ2 - x3ξ2 .* x1ξ1) ./ J
            ξ3x3 = (x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1) ./ J

            foreach(
                j -> (x1[j], x2[j], x3[j]) = f(x1[j], x2[j], x3[j]),
                1:length(x1),
            )

            metric = Metrics.computemetric(x1, x2, x3, D)

            @test metric.J ≈ J
            @test metric.ξ1x1 ≈ ξ1x1
            @test metric.ξ1x2 ≈ ξ1x2
            @test metric.ξ1x3 ≈ ξ1x3
            @test metric.ξ2x1 ≈ ξ2x1
            @test metric.ξ2x2 ≈ ξ2x2
            @test metric.ξ2x3 ≈ ξ2x3
            @test metric.ξ3x1 ≈ ξ3x1
            @test metric.ξ3x2 ≈ ξ3x2
            @test metric.ξ3x3 ≈ ξ3x3
            ind = LinearIndices((1:Nq, 1:Nq))
            n1 = metric.n1
            n2 = metric.n2
            n3 = metric.n3
            sJ = metric.sJ
            @test hypot.(n1, n2, n3) ≈ ones(T, size(n1))
            @test (
                [
                    sJ[ind, 1, :] .* n1[ind, 1, :],
                    sJ[ind, 1, :] .* n2[ind, 1, :],
                    sJ[ind, 1, :] .* n3[ind, 1, :],
                ] ≈ [
                    -J[1, :, :, :] .* ξ1x1[1, :, :, :],
                    -J[1, :, :, :] .* ξ1x2[1, :, :, :],
                    -J[1, :, :, :] .* ξ1x3[1, :, :, :],
                ]
            )
            @test (
                [
                    sJ[ind, 2, :] .* n1[ind, 2, :],
                    sJ[ind, 2, :] .* n2[ind, 2, :],
                    sJ[ind, 2, :] .* n3[ind, 2, :],
                ] ≈ [
                    J[Nq, :, :, :] .* ξ1x1[Nq, :, :, :],
                    J[Nq, :, :, :] .* ξ1x2[Nq, :, :, :],
                    J[Nq, :, :, :] .* ξ1x3[Nq, :, :, :],
                ]
            )
            @test sJ[ind, 3, :] .* n1[ind, 3, :] ≈
                  -J[:, 1, :, :] .* ξ2x1[:, 1, :, :]
            @test sJ[ind, 3, :] .* n2[ind, 3, :] ≈
                  -J[:, 1, :, :] .* ξ2x2[:, 1, :, :]
            @test sJ[ind, 3, :] .* n3[ind, 3, :] ≈
                  -J[:, 1, :, :] .* ξ2x3[:, 1, :, :]
            @test sJ[ind, 4, :] .* n1[ind, 4, :] ≈
                  J[:, Nq, :, :] .* ξ2x1[:, Nq, :, :]
            @test sJ[ind, 4, :] .* n2[ind, 4, :] ≈
                  J[:, Nq, :, :] .* ξ2x2[:, Nq, :, :]
            @test sJ[ind, 4, :] .* n3[ind, 4, :] ≈
                  J[:, Nq, :, :] .* ξ2x3[:, Nq, :, :]
            @test sJ[ind, 5, :] .* n1[ind, 5, :] ≈
                  -J[:, :, 1, :] .* ξ3x1[:, :, 1, :]
            @test sJ[ind, 5, :] .* n2[ind, 5, :] ≈
                  -J[:, :, 1, :] .* ξ3x2[:, :, 1, :]
            @test sJ[ind, 5, :] .* n3[ind, 5, :] ≈
                  -J[:, :, 1, :] .* ξ3x3[:, :, 1, :]
            @test sJ[ind, 6, :] .* n1[ind, 6, :] ≈
                  J[:, :, Nq, :] .* ξ3x1[:, :, Nq, :]
            @test sJ[ind, 6, :] .* n2[ind, 6, :] ≈
                  J[:, :, Nq, :] .* ξ3x2[:, :, Nq, :]
            @test sJ[ind, 6, :] .* n3[ind, 6, :] ≈
                  J[:, :, Nq, :] .* ξ3x3[:, :, Nq, :]
        end
    end
    #}}}

    # Constant preserving test
    #{{{
    for T in (Float32, Float64)
        f(r, s, t) = @.( (
            s + r * t - (r^2 * s^2 * t^2) / 4,
            t - ((r * s * t) / 2 + 1 / 2)^3 + 1,
            r + (r / 2 + 1 / 2)^6 * (s / 2 + 1 / 2)^6 * (t / 2 + 1 / 2)^6,
        ))
        let
            N = 5

            r, w = Elements.lglpoints(T, N)
            D = Elements.spectralderivative(r)
            Nq = N + 1

            d = 3
            nfaces = 6
            e2c = Array{T, 3}(undef, d, 8, 1)
            e2c[:, :, 1] = [
                -1 1 -1 1 -1 1 -1 1
                -1 -1 1 1 -1 -1 1 1
                -1 -1 -1 -1 1 1 1 1
            ]

            nelem = size(e2c, 3)

            vgeo = Array{T, 5}(undef, Nq, Nq, Nq, length(VGEO3D), nelem)
            sgeo = Array{T, 4}(undef, Nq^2, nfaces, length(SGEO3D), nelem)
            Metrics.creategrid!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), d)...,
                e2c,
                r,
            )
            x1 = @view vgeo[:, :, :, VGEO3D.x1, :]
            x2 = @view vgeo[:, :, :, VGEO3D.x2, :]
            x3 = @view vgeo[:, :, :, VGEO3D.x3, :]

            foreach(
                j -> (x1[j], x2[j], x3[j]) = f(x1[j], x2[j], x3[j]),
                1:length(x1),
            )

            Metrics.computemetric!(
                ntuple(j -> (@view vgeo[:, :, :, j, :]), length(VGEO3D))...,
                ntuple(j -> (@view sgeo[:, :, j, :]), length(SGEO3D))...,
                D,
            )
            sgeo = reshape(sgeo, Nq, Nq, nfaces, length(SGEO3D), nelem)

            (Cx1, Cx2, Cx3) = (
                zeros(T, Nq, Nq, Nq),
                zeros(T, Nq, Nq, Nq),
                zeros(T, Nq, Nq, Nq),
            )

            J = @view vgeo[:, :, :, VGEO3D.J, :]
            ξ1x1 = @view vgeo[:, :, :, VGEO3D.ξ1x1, :]
            ξ1x2 = @view vgeo[:, :, :, VGEO3D.ξ1x2, :]
            ξ1x3 = @view vgeo[:, :, :, VGEO3D.ξ1x3, :]
            ξ2x1 = @view vgeo[:, :, :, VGEO3D.ξ2x1, :]
            ξ2x2 = @view vgeo[:, :, :, VGEO3D.ξ2x2, :]
            ξ2x3 = @view vgeo[:, :, :, VGEO3D.ξ2x3, :]
            ξ3x1 = @view vgeo[:, :, :, VGEO3D.ξ3x1, :]
            ξ3x2 = @view vgeo[:, :, :, VGEO3D.ξ3x2, :]
            ξ3x3 = @view vgeo[:, :, :, VGEO3D.ξ3x3, :]

            e = 1
            for m in 1:Nq
                for n in 1:Nq
                    Cx1[:, n, m] += D * (J[:, n, m, e] .* ξ1x1[:, n, m, e])
                    Cx1[n, :, m] += D * (J[n, :, m, e] .* ξ2x1[n, :, m, e])
                    Cx1[n, m, :] += D * (J[n, m, :, e] .* ξ3x1[n, m, :, e])

                    Cx2[:, n, m] += D * (J[:, n, m, e] .* ξ1x2[:, n, m, e])
                    Cx2[n, :, m] += D * (J[n, :, m, e] .* ξ2x2[n, :, m, e])
                    Cx2[n, m, :] += D * (J[n, m, :, e] .* ξ3x2[n, m, :, e])

                    Cx3[:, n, m] += D * (J[:, n, m, e] .* ξ1x3[:, n, m, e])
                    Cx3[n, :, m] += D * (J[n, :, m, e] .* ξ2x3[n, :, m, e])
                    Cx3[n, m, :] += D * (J[n, m, :, e] .* ξ3x3[n, m, :, e])
                end
            end
            @test maximum(abs.(Cx1)) ≤ 1000 * eps(T)
            @test maximum(abs.(Cx2)) ≤ 1000 * eps(T)
            @test maximum(abs.(Cx3)) ≤ 1000 * eps(T)
        end
    end
    #}}}
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using BackwardsLinalg
using Test, Random

@testset "qr Q complex" begin
    T = ComplexF64
    Random.seed!(3)
    for (M, N) in [(2, 6), (4, 4), (6, 2)]
        A = randn(T, M, N)
        op = randn(M, M)
        op += op'
        op2 = randn(N, N)
        op2 += op2'

        function tfunc(x)
            Q, R = qr(x)
            v = Q[:,1]
            v2 = R[2,:]
            (v'*op*v + v2'*op2*v2)[] |> real
        end
        @test gradient_check(tfunc, A)
    end
    a = [1+1im 2+1im; 3-1im 4+2im]
    @test copyltu!(a) ≈ [1 3+1im; 3-1im 4]
end

@testset "lq" begin
    Random.seed!(3)
    T = ComplexF64
    for (M, N) in [(2, 6), (4, 4), (6, 2)]
        @show M, N
        A = randn(T, M, N)

        op = randn(M, M)
        op += op'
        op2 = randn(N, N)
        op2 += op2'

        function tfunc(x)
            L, Q = lq(x)
            v = L[:,1]
            v2 = Q[2,:]
            (v'*op*v + v2'*op2*v2)[] |> real
        end
        @test gradient_check(tfunc, A)
    end
end

@testset "qr Q complex pivot" begin
    T = Float64
    Random.seed!(3)
    for (M, N) in [(20, 60), (40, 40), (60, 20)]
        A = randn(T, M, N)
        op = randn(M, M)
        op += op'
        op2 = randn(N, N)
        op2 += op2'

        function tfunc(x)
            Q, R, P = qr(x, Val(true))
            v = Q[:,1]
            v2 = R[2,:]
            @show size(v), size(v2), size(op), size(op2)
            (v'*op*v + v2'*op2*v2)[] |> real
        end
        @show tfunc(A)
        @test gradient_check(tfunc, A)
    end
    a = [1+1im 2+1im; 3-1im 4+2im]
    @test copyltu!(a) ≈ [1 3+1im; 3-1im 4]
end


text:
"""
This yields:
0: 3
1: 0
2: 0
3: 0
4: 4
5: 4
6: 3
7: 2
8: 3
9: 3
10: 4
11: 5
12: 0
13: 0
14: 0
15: 0
16: 6
17: 7
18: 7

At this point we need to understand how to go from these values to
bit strings. This is done with Huffman coding.
https://en.wikipedia.org/wiki/Huffman_coding

I won't go into too much detail but the important idea is that we can git a
bunch of bit strings whereby each string can be uniquely identified.

For example the following, when read left to right are unique:
010
001
000

The following however would not be ok:
01
01011

Since if we saw 01 we wouldn't know if it was its own code or if it was a
subset of 01011

The algorithm is as follows:

- Start with 0 at the minimum length.
- Increment by 1 for every code at that length
- Once done go back to the next length and start at 2*last_value + 1
"""


text:
function randomstress(npts::Int = 100, ndgs::Int = 200)

    V = rand(2, npts)
    model = Lar.Model(rand(2, npts))

    cv = unique([rand(1 : npts, 2) for i = 1 : 10*npts])[1 : ndgs*2]
    cv = [[c[1] c[2]] for c in cv if c[1] != c[2]][1 : ndgs]
    I = []
    J = []
    K = ones(Int8, 2*ndgs)
    for i = 1 : ndgs
      push!(J, cv[i]...)
      push!(I, [i; i]...)
    end
    Lar.addModelCells!(model, 1, SparseArrays.sparse(I, J, K, ndgs, npts))

    GL.VIEW([ GL.GLLines(model.G, LAR.cop2lar(model.T[1])) ])

    model1, edge_map = LarA.planar_arrangement_1(model, spzeros(Int8, 0), true)

    W, EW = LAR.fraglines(1.5,1.5,1.5)((model1.G, LAR.cop2lar(model1.T[1])))
    GL.VIEW([ GL.GLLines(W, EW, GL.COLORS[1]) ])

    bicon_comps = LarA.biconnected_components(model1.T[1])

	#V, EVs = LAR.biconnectedComponent((model1.G, LAR.cop2lar(model1.T[1]))) # 2-connected components
	#VV = [[k] for k = 1 : size(V,2)]
	#meshes = [ GL.numbering(.01)((V, (VV,EVs[k])), GL.COLORS[k],1) for k=1:length(EVs) ]
	#GL.VIEW( cat(meshes) );

    edge_map = LarA.remove_mapping_dangling_edges(edge_map, bicon_comps)
    model2 = LarA.planar_arrangement_2(model1, bicon_comps)

    EVs = Lar.FV2EVs(model2)
	FVs = Lar.triangulate2D(model2)
	GL.VIEW(GL.GLExplode(model2.G, FVs, 1.2, 1.2, 1.2, 99, 1));
	GL.VIEW(GL.GLExplode(model2.G, FVs,   1,   1,   1, 99, 1));
	GL.VIEW(GL.GLExplode(model2.G, EVs,   1,   1,   1, 99, 1));

end


text:
export twins

"""
`twins(G,u,v)` determines if `u` and `v` are twin vertices of `G`.
That is, if `G[u]-v == G[v]-u`. This is an equivalence relation.

`twins(G)` returns a partition of the graph's vertex set into twin
equivalence classes.
"""
function twins(G::SimpleGraph,u,v)::Bool
    @assert has(G,u) "vertex $u not in the graph"
    @assert has(G,v) "vertex $v not in the graph"

    if u==v
        return true
    end

    Nu = Set(G[u])
    Nv = Set(G[v])

    SimpleGraphs.delete!(Nu,v)
    SimpleGraphs.delete!(Nv,u)

    return Nu==Nv
end

function twins(G)::Partition
    if cache_check(G,:twins)
      return cache_recall(G,:twins)
    end
    VV = vlist(G)
    n  = NV(G)

    P = Partition(VV)
    for i=1:n-1
        u = VV[i]
        for j=i+1:n
            v = VV[j]
            if twins(G,u,v)
                merge_parts!(P,u,v)
            end
        end
    end
    cache_save(G,:twins,P)
    return P
end


text:
using ProceduralNoise
using Plots

scale = (2, 2) .^ 2
resolution = (2, 2) .^ 9
A = perlin(resolution, scale)
heatmap(A)


text:
using Test
using DFTK
using Random

#
# This test suite test arguments. For example:
#     Pkg.test("DFTK"; test_args = ["fast"])
# only runs the "fast" tests (i.e. not the expensive ones)
#     Pkg.test("DFTK"; test_args = ["example"])
# runs only the tests tagged as "example" and
#     Pkg.test("DFTK"; test_args = ["example", "all"])
# runs all tests plus the "example" tests.
#

# By default run expensive tests, but not if in CI environment
# If user supplies the "fast" tag
const FAST_TESTS = ifelse("CI" in keys(ENV), parse(Bool, get(ENV, "CI", "false")),
                          "fast" in ARGS)

# Tags supplied by the user ... filter out "fast" (already dealt with)
TAGS = filter(e -> !(e in ["fast"]), ARGS)
isempty(TAGS) && (TAGS = ["all"])

if FAST_TESTS
    println("   Running fast tests (TAGS = $(join(TAGS, ", "))).")
else
    println("   Running tests (TAGS = $(join(TAGS, ", "))).")
end

# Initialize seed
Random.seed!(0)

# Wrap in an outer testset to get a full report if one test fails
@testset "DFTK.jl" begin
    # Synthetic tests at the beginning, so it fails faster if
    # something has gone badly wrong
    if "all" in TAGS || "functionality" in TAGS
        include("hydrogen_all_electron.jl")
        include("silicon_redHF.jl")
        include("silicon_lda.jl")
        include("silicon_pbe.jl")
        include("scf_compare.jl")
    end

    if "all" in TAGS
        include("determine_grid_size.jl")
        include("fourier_transforms.jl")
        include("PlaneWaveBasis.jl")
        include("interpolation.jl")
        include("load_psp.jl")
        include("PspHgh.jl")
        include("elements.jl")
        include("bzmesh.jl")
        include("bzmesh_symmetry.jl")
        include("external_pymatgen.jl")
        include("external_ase.jl")
    end

    if "all" in TAGS
        include("hamiltonian_consistency.jl")
    end

    if "all" in TAGS
        include("lobpcg.jl")
        include("diag_compare.jl")
        include("xc_fallback.jl")
        include("interval_arithmetic.jl")
    end

    if "all" in TAGS
        include("ewald.jl")
        include("energy_nuclear.jl")
        include("occupation.jl")
        include("energies_guess_density.jl")
        include("compute_density.jl")
        include("forces.jl")
    end

    if "all" in TAGS
        include("variational.jl")
        include("compute_bands.jl")
        include("chi0.jl")
        include("kernel.jl")
    end

    ("example" in TAGS) && include("runexamples.jl")
end


text:
module MangalTestDataset
    using Mangal
    using Test
    include("setup.jl")

    # Datasets with no arguments
    @test typeof(datasets()) <: Vector{MangalDataset}

    # Datasets with count
    count_2 = datasets(Pair("count", 2))
    @test typeof(count_2) <: Vector{MangalDataset}
    @test length(count_2) == 2

    # Dataset by id
    @test typeof(dataset(1)) <: MangalDataset

    # Dataset by name
    @test typeof(dataset("roberson_1929")) <: MangalDataset

    @test isnothing(dataset("This_DOESNOT_exist"))

end


text:
type Payer <: SwapType end
type Receiver <: SwapType end

type Buyer <: CDSProtectionSide end
type Seller <: CDSProtectionSide end

type SwapResults <: Results
  legNPV::Vector{Float64}
  legBPS::Vector{Float64}
  npvDateDiscount::Float64
  startDiscounts::Vector{Float64}
  endDiscounts::Vector{Float64}
  fairRate::Float64
  value::Float64

  function SwapResults{I <: Integer}(n::I)
    legNPV = zeros(n)
    legBPS = zeros(n)
    startDiscounts = zeros(n)
    endDiscounts = zeros(n)

    new(legNPV, legBPS, 0.0, startDiscounts, endDiscounts, -1.0, 0.0)
  end
end

type VanillaSwapArgs
  fixedResetDates::Vector{Date}
  fixedPayDates::Vector{Date}
  floatingResetDates::Vector{Date}
  floatingPayDates::Vector{Date}
  floatingAccrualTimes::Vector{Float64}
  floatingSpreads::Vector{Float64}
  fixedCoupons::Vector{Float64}
  floatingCoupons::Vector{Float64}
end

function VanillaSwapArgs{L <: Leg}(legs::Vector{L})
  fixedCoups = legs[1].coupons
  floatingCoups = legs[2].coupons
  fixedCoupons = [amount(coup) for coup in fixedCoups]
  # floatingCoupons = [amount(coup) for coup in floatingCoups]
  floatingCoupons = zeros(length(floatingCoups))
  floatingAccrualTimes = [accrual_period(coup) for coup in floatingCoups]
  floatingSpreads = [coup.spread for coup in floatingCoups]
  return VanillaSwapArgs(get_reset_dates(fixedCoups), get_pay_dates(fixedCoups), get_reset_dates(floatingCoups), get_pay_dates(floatingCoups), floatingAccrualTimes, floatingSpreads, fixedCoupons, floatingCoupons)
end

function reset!(sr::SwapResults)
  n = length(sr.legNPV)
  sr.legNPV = zeros(n)
  sr.legBPS = zeros(n)
  sr.npvDateDiscount = 0.0
  sr.startDiscounts = zeros(n)
  sr.endDiscounts = zeros(n)
  sr.value = 0.0
  sr.fairRate = 0.0

  return sr
end

type NonstandardSwapArgs
  vSwapArgs::VanillaSwapArgs
  fixedIsRedemptionFlow::Vector{Bool}
  floatingIsRedemptionFlow::Vector{Bool}
  floatingGearings::Vector{Float64}
end

function NonstandardSwapArgs{L <: Leg}(legs::Vector{L})
  vSwapArgs = VanillaSwapArgs(legs)
  fixedIsRedemptionFlow = Bool[!check_coupon(cf) for cf in legs[1].coupons]
  floatingIsRedemptionFlow = Bool[!check_coupon(cf) for cf in legs[2].coupons]
  floatingGearings = get_gearings(legs[2].coupons)

  return NonstandardSwapArgs(vSwapArgs, fixedIsRedemptionFlow, floatingIsRedemptionFlow, floatingGearings)
end

type CDSResults
  upfrontNPV::Float64
  couponLegNPV::Float64
  defaultLegNPV::Float64
  fairSpread::Float64
  fairUpfront::Float64
  couponLegBPS::Float64
  upfrontBPS::Float64
  value::Float64
end

CDSResults() = CDSResults(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

function reset!(cr::CDSResults)
  cr.upfrontNPV = 0.0
  cr.couponLegNPV = 0.0
  cr.defaultLegNPV = 0.0
  cr.fairSpread = 0.0
  cr.fairUpfront = 0.0
  cr.couponLegBPS = 0.0
  cr.upfrontBPS = 0.0
  cr.value = 0.0

  return cr
end

type VanillaSwap{ST <: SwapType, DC_fix <: DayCount, DC_float <: DayCount, B <: BusinessDayConvention, L <: Leg, P <: PricingEngine} <: Swap
  lazyMixin::LazyMixin
  swapT::ST
  nominal::Float64
  fixedSchedule::Schedule
  fixedRate::Float64
  fixedDayCount::DC_fix
  iborIndex::IborIndex
  spread::Float64
  floatSchedule::Schedule
  floatDayCount::DC_float
  paymentConvention::B
  legs::Vector{L}
  payer::Vector{Float64}
  pricingEngine::P
  results::SwapResults
  args::VanillaSwapArgs
end

# Constructors
function VanillaSwap{ST <: SwapType, DC_fix <: DayCount, DC_float <: DayCount, B <: BusinessDayConvention, P <: PricingEngine}(swapT::ST, nominal::Float64, fixedSchedule::Schedule, fixedRate::Float64,
                    fixedDayCount::DC_fix, iborIndex::IborIndex, spread::Float64, floatSchedule::Schedule, floatDayCount::DC_float, pricingEngine::P = NullPricingEngine(), paymentConvention::B = floatSchedule.convention)
  # build swap cashflows
  legs = Vector{Leg}(2)
  # first leg is fixed
  legs[1] = FixedRateLeg(fixedSchedule, nominal, fixedRate, fixedSchedule.cal, paymentConvention, fixedDayCount; add_redemption=false)
  # second leg is floating
  legs[2] = IborLeg(floatSchedule, nominal, iborIndex, floatDayCount, paymentConvention; add_redemption=false)

  payer = _build_payer(swapT)

  results = SwapResults(2)

  return VanillaSwap{ST, DC_fix, DC_float, B, Leg, P}(LazyMixin(), swapT, nominal, fixedSchedule, fixedRate, fixedDayCount, iborIndex, spread, floatSchedule, floatDayCount,
                    paymentConvention, legs, payer, pricingEngine, results, VanillaSwapArgs(legs))
end

# accessor methods #
get_fixed_reset_dates(swap::VanillaSwap) = swap.args.fixedResetDates
get_fixed_pay_dates(swap::VanillaSwap) = swap.args.fixedPayDates
get_floating_reset_dates(swap::VanillaSwap) = swap.args.floatingResetDates
get_floating_pay_dates(swap::VanillaSwap) = swap.args.floatingPayDates
get_floating_accrual_times(swap::VanillaSwap) = swap.args.floatingAccrualTimes
get_floating_spreads(swap::VanillaSwap) = swap.args.floatingSpreads
get_fixed_coupons(swap::VanillaSwap) = swap.args.fixedCoupons
get_floating_coupons(swap::VanillaSwap) = swap.args.floatingCoupons

type NonstandardSwap{ST <: SwapType, DC_fix <: DayCount, DC_float <: DayCount, B <: BusinessDayConvention, L <: Leg, P <: PricingEngine} <: Swap
  lazyMixin::LazyMixin
  swapT::ST
  fixedNominal::Vector{Float64}
  floatingNominal::Vector{Float64}
  fixedSchedule::Schedule
  fixedRate::Vector{Float64}
  fixedDayCount::DC_fix
  iborIndex::IborIndex
  spread::Float64
  gearing::Float64
  floatSchedule::Schedule
  floatDayCount::DC_float
  paymentConvention::B
  intermediateCapitalExchange::Bool
  finalCapitalExchange::Bool
  legs::Vector{L}
  payer::Vector{Float64}
  pricingEngine::P
  results::SwapResults
  args::NonstandardSwapArgs
end

# Constructor #
function NonstandardSwap(vs::VanillaSwap)
  # build swap cashflows
  legs = Vector{Leg}(2)
  # first leg is fixed
  legs[1] = FixedRateLeg(vs.fixedSchedule, vs.nominal, vs.fixedRate, vs.fixedSchedule.cal, vs.paymentConvention, vs.fixedDayCount; add_redemption=false)
  # second leg is floating
  legs[2] = IborLeg(vs.floatSchedule, vs.nominal, vs.iborIndex, vs.floatDayCount, vs.paymentConvention; add_redemption=false)

  payer = _build_payer(vs.swapT)

  results = SwapResults(2)

  fixedSize = length(vs.legs[1].coupons)
  floatSize = length(vs.legs[2].coupons)

  return NonstandardSwap(LazyMixin(), vs.swapT, fill(vs.nominal, fixedSize), fill(vs.nominal, floatSize), vs.fixedSchedule, fill(vs.fixedRate, fixedSize),
                        vs.fixedDayCount, vs.iborIndex, vs.spread, 1.0, vs.floatSchedule, vs.floatDayCount, vs.paymentConvention, false, false, legs,
                        payer, vs.pricingEngine, results, NonstandardSwapArgs(legs))
end

# accessor methods #
get_fixed_reset_dates(swap::NonstandardSwap) = swap.args.vSwapArgs.fixedResetDates
get_fixed_pay_dates(swap::NonstandardSwap) = swap.args.vSwapArgs.fixedPayDates
get_floating_reset_dates(swap::NonstandardSwap) = swap.args.vSwapArgs.floatingResetDates
get_floating_pay_dates(swap::NonstandardSwap) = swap.args.vSwapArgs.floatingPayDates
get_floating_accrual_times(swap::NonstandardSwap) = swap.args.vSwapArgs.floatingAccrualTimes
get_floating_spreads(swap::NonstandardSwap) = swap.args.vSwapArgs.floatingSpreads
get_fixed_coupons(swap::NonstandardSwap) = swap.args.vSwapArgs.fixedCoupons
get_floating_coupons(swap::NonstandardSwap) = swap.args.vSwapArgs.floatingCoupons

# CDS #
type CreditDefaultSwap{S <: CDSProtectionSide, B <: BusinessDayConvention, DC <: DayCount, P <: PricingEngine} <: Swap
  lazyMixin::LazyMixin
  side::S
  notional::Float64
  spread::Float64
  schedule::Schedule
  convention::B
  dc::DC
  leg::FixedRateLeg
  upfrontPayment::SimpleCashFlow
  settlesAccrual::Bool
  paysAtDefaultTime::Bool
  protectionStart::Date
  pricingEngine::P
  claim::FaceValueClaim
  results::CDSResults

  CreditDefaultSwap(lazyMixin::LazyMixin,
                    side::S,
                    notional::Float64,
                    spread::Float64,
                    schedule::Schedule,
                    convention::B,
                    dc::DC,
                    leg::FixedRateLeg,
                    upfrontPayment::SimpleCashFlow,
                    settlesAccrual::Bool,
                    paysAtDefaultTime::Bool,
                    protectionStart::Date,
                    pricingEngine::P) = new(lazyMixin, side, notional, spread, schedule, convention, dc, leg, upfrontPayment, settlesAccrual, paysAtDefaultTime, protectionStart, pricingEngine,
                                            FaceValueClaim(), CDSResults())
end

function CreditDefaultSwap{S <: CDSProtectionSide, B <: BusinessDayConvention, DC <: DayCount, P <: PricingEngine}(side::S, notional::Float64, spread::Float64, schedule::Schedule,
                          convention::B, dc::DC, settlesAccrual::Bool, paysAtDefaultTime::Bool, protectionStart::Date, pricingEngine::P)
  # build leg
  leg = FixedRateLeg(schedule, notional, spread, schedule.cal, convention, dc; add_redemption = false)

  # build upfront payment
  upfrontPayment = SimpleCashFlow(0.0, schedule.dates[1])
  return CreditDefaultSwap{S, B, DC, P}(LazyMixin(), side, notional, spread, schedule, convention, dc, leg, upfrontPayment, settlesAccrual, paysAtDefaultTime, protectionStart, pricingEngine)
end

# Swap Helper methods
function _build_payer(swapT::Payer)
  x = ones(2)
  x[1] = -1.0
  return x
end

function _build_payer(swapT::Receiver)
  x = ones(2)
  x[2] = -1.0
  return x
end

# Swap methods #
function maturity_date(swap::Swap)
  d = maturity_date(swap.legs[1])
  for i = 2:length(swap.legs)
    d = max(d, maturity_date(swap.legs[i]))
  end

  return d
end

# Calculation method #
function perform_calculations!(swap::VanillaSwap)
  reset!(swap.results) # reset - TODO this will be expanded
  swap.args.floatingCoupons = [amount(coup) for coup in swap.legs[2].coupons]
  _calculate!(swap.pricingEngine, swap)

  return swap
end

function perform_calculations!(cds::CreditDefaultSwap)
  reset!(cds.results)
  _calculate!(cds.pricingEngine, cds)

  return cds
end

floating_leg_NPV(swap::VanillaSwap) = swap.results.legNPV[2]
floating_leg_BPS(swap::VanillaSwap) = swap.results.legBPS[2]

fixed_leg_BPS(swap::VanillaSwap) = swap.results.legBPS[1]

function fair_rate(swap::VanillaSwap)
  calculate!(swap)

  return swap.results.fairRate
end

# some helper methods #
# function clone(swap::VanillaSwap, pe::PricingEngine = swap.pricingEngine)
#   lazyMixin, res, args = pe == swap.pricingEngine ? (swap.lazyMixin, swap.results, swap.args) : (LazyMixin(), SwapResults(2), VanillaSwapArgs(swap.legs))
#   # args = VanillaSwapArgs(swap.legs)
#   # res = SwapResults(2)
#
#   return VanillaSwap(lazyMixin, swap.swapT, swap.nominal, swap.fixedSchedule, swap.fixedRate, swap.fixedDayCount, swap.iborIndex, swap.spread,
#                     swap.floatSchedule, swap.floatDayCount, swap.paymentConvention, swap.legs, swap.payer, pe, res, args)
# end

function clone(swap::VanillaSwap, pe::PricingEngine = swap.pricingEngine, ts::TermStructure = swap.iborIndex.ts)
  # is_new = pe != swap.pricingEngine || ts != swap.iborIndex.ts

  lazyMixin, res, args = pe == swap.pricingEngine ? (swap.lazyMixin, swap.results, swap.args) : (LazyMixin(), SwapResults(2), VanillaSwapArgs(swap.legs))

  # we need a new ibor and to rebuild floating rate coupons
  if ts != swap.iborIndex.ts
    newIbor = clone(swap.iborIndex, ts)
    newLegs = Vector{Leg}(2)
    newLegs[1] = swap.legs[1]
    newLegs[2] = IborLeg(swap.floatSchedule, swap.nominal, newIbor, swap.floatDayCount, swap.paymentConvention; add_redemption=false)
  else
    newIbor = swap.iborIndex
    newLegs = swap.legs
  end

  return VanillaSwap(lazyMixin, swap.swapT, swap.nominal, swap.fixedSchedule, swap.fixedRate, swap.fixedDayCount, newIbor, swap.spread,
                    swap.floatSchedule, swap.floatDayCount, swap.paymentConvention, newLegs, swap.payer, pe, res, args)
end

get_pricing_engine_type{ST, DC_fix, DC_float, B, L, P}(::VanillaSwap{ST, DC_fix, DC_float, B, L, P}) = P

function update_ts_idx!(swap::VanillaSwap, ts::TermStructure)
  typeof(ts) == typeof(swap.iborIndex.ts) || error("Term Structure mismatch for swap between ts and index ts")
  newIborIdx = clone(swap.iborIndex, ts)
  swap.iborIndex = newIborIdx

  # update legs
  for coup in swap.legs[2].coupons
    coup.iborIndex = newIborIdx
  end

  swap.lazyMixin.calculated = false

  return swap
end

function update_ts_pe!(swap::VanillaSwap, ts::TermStructure)
  typeof(ts) == typeof(swap.pricingEngine.yts) || error("Term Structure mismatch for swap between ts and pric engine ts")

  swap.pricingEngine.yts = ts

  swap.lazyMixin.calculated = false

  return swap
end

function update_all_ts!(swap::VanillaSwap, ts::TermStructure)
  # this will update the ts of the pricing engine and ibor index
  update_ts_idx!(swap, ts)
  update_ts_pe!(swap, ts)

  return swap
end


text:
using Genie
using Genie.Router, Genie.Renderer, Genie.Renderer.Html, Genie.Renderer.Json, Genie.Requests
using MySQL
using JSON
using DataFrames
using JSONTables

create_table_sql = """
    CREATE TABLE IF NOT EXISTS USERS(
        user_id VARCHAR(50) NOT NULL PRIMARY KEY,
        user_status ENUM('active', 'inactive') DEFAULT 'active'
    );
    """
conn = DBInterface.connect(MySQL.Connection, "localhost", "<user_name>", "<password>", db = "test_user")

cur = DBInterface.execute(conn, create_table_sql)

function createUser(user_id, user_status)

DBInterface.execute(conn, "INSERT INTO USERS (user_id, user_status) VALUES ('$user_id', '$user_status')")
end

function getUsers()
    cur = DBInterface.execute(conn, "SELECT user_id, user_status FROM users")
    df = DataFrame(cur)
    stringified_json = arraytable(df)
    return stringified_json

end

route("/createuser", method = POST) do

    createUser(postpayload(:user_id), postpayload(:user_status, "active"))

    return "POST OK"

end

route("/getusers") do

    obj = getUsers()

    return obj

end

up(8001, async = false)


text:
"""
    householderqr(R1, R2) computes the QR decompsition of
    (m+l)-by-l block matrix [R1;R2]

```math
          l               l
    m  (  R1  )  =  m  (  Q1  ) * R
    l  (  R2  )     l  (  Q2  )
```

- Q1 is m-by-l orthogonal matrix,
- Q2 is l-by-l orthogonal matrix,
- R is l-by-l upper triangular matrix.

# arguments:
R1: m-by-l matrix,
R2: m-by-l matrix.

# returns:
Q1: m-by-l matrix,
Q2: l-by-l matrix,
R: l-by-l matrix.
"""
function householderqr(R1, R2)
    r, c = size(R1)
    F = qr!([R1;R2])
    Q = Matrix(F.Q)
    Q1 = Q[1:r,:]
    Q2 = Q[r+1:end,:]
    return Q1, Q2, F.R
end


text:
# lbm.jl
# Lattice Boltzmann Method (D2Q9) written in Julia
# Vincent San Miguel

# Rectangular domain dimensions, x₁ in East and x₂ in North
x₁ = 300
x₂ = 100

# Viscosity
ν = 0.02
# Relaxation parameter
ω = 1 / (3 * ν + 0.5)
# Inlet velocity
u₀ = 0.1

# For a 2D gas, the continuous distribution of the thermal velocities, v,
# is given by the Boltzmann distribution
#
# Lattice is D2Q9. We specify 9 velocity vectors e₀ through e₈
# corresponding to zero, the 4 cardinal directions, and 4 diagonal directions.
# We can attach probabilities to these vectors. Since the fluid is at rest,
# the probability distribution is symmetrical, and velocities with equal
# magnitude - regardless of direction - will have equal weightings.
#
# Lattice total velocities, e = Δx / Δt
# e₀ = [0, 0], e₁ = [1, 0] (E), e₂ = [0, 1] (N), e₃ = [-1, 0] (W),
# e₄ = [0, -1] (S), e₅ = [1, 1] (NE), e₆ = [-1, 1] (NW), e₇ = [-1, -1] (SW),
# e₈ = [1, -1] (SE)
#
# We can then choose probability weightings, w, to approximate the
# continuous Boltzmann distribution such that these average velocity
# values are preserved

# Lattice probability weightings
w₀ = 4.0 / 9.0
w₁ = w₂ = w₃ = w₄ = 1.0 / 9.0
w₅ = w₆ = w₇ = w₈ = 1.0 / 36.0

# We can obtain number densities for a fluid in motion by rewriting the
# continuous Boltzmann distribution, replacing the thermal velocity with
# the scaled total velocity minus the macroscopic velocity, v = e * c - u,
# and taking the Taylor series

# Initial number densities - solution domain initialized with steady
# rightward flow
n₀ = w₀ * ones(x₂, x₁) - 1.5 * u₀^2
n₁ = w₁ * ones(x₂, x₁) + 3.0 * u₀ + 4.5 * u₀^2 - 1.5 * u₀^2
n₂ = w₂ * ones(x₂, x₁) - 1.5 * u₀^2
n₃ = w₃ * ones(x₂, x₁) - 3.0 * u₀ + 4.5 * u₀^2 - 1.5 * u₀^2
n₄ = w₄ * ones(x₂, x₁) - 1.5 * u₀^2
n₅ = w₅ * ones(x₂, x₁) + 3.0 * u₀ + 4.5 * u₀^2 - 1.5 * u₀^2
n₆ = w₆ * ones(x₂, x₁) - 3.0 * u₀ + 4.5 * u₀^2 - 1.5 * u₀^2
n₇ = w₇ * ones(x₂, x₁) - 3.0 * u₀ + 4.5 * u₀^2 - 1.5 * u₀^2
n₈ = w₈ * ones(x₂, x₁) + 3.0 * u₀ + 4.5 * u₀^2 - 1.5 * u₀^2

# Macroscopic density
ρ = n₀ + n₁ + n₂ + n₃ + n₄ + n₅ + n₆ + n₇ + n₈

# Macroscopic x-velocity and y-velocity
ux = ((n₁ + n₅ + n₈) - (n₃ + n₆ + n₇)) / ρ
uy = ((n₂ + n₅ + n₆) - (n₄ + n₇ + n₈)) / ρ

# Initialize barriers
barrier = falses(x₂, x₁)

barrier[convert(Int, 0.5 * x₂ - 8) : convert(Int, 0.5 * x₂ + 8),
        convert(Int, 0.5 * x₂)] = true

# Sites around barrier
b₁ = circshift(barrier, [0, 1])
b₂ = circshift(barrier, [1, 0])
b₃ = circshift(barrier, [0, -1])
b₄ = circshift(barrier, [-1, 0])
b₅ = circshift(barrier, [1, 1])
b₆ = circshift(barrier, [1, -1])
b₇ = circshift(barrier, [-1, -1])
b₈ = circshift(barrier, [-1, 1])

function collide!()
    # Collides all particles

    rho = n₀ + n₁ + n₂ + n₃ + n₄ + n₅ + n₆ + n₇ + n₈

    ux = ((n₁ + n₅ + n₈) - (n₃ + n₆ + n₇)) / ρ
    uy = ((n₂ + n₅ + n₆) - (n₄ + n₇ + n₈)) / ρ

    ux² = ux^2
    uy² = uy^2
    u² = ux² * uy²
    uxuy = ux * uy

    n₀ = (1 - ω) * n₀ + ω * w₀ * ρ * (1 - 1.5 * u²)
    n₁ = (1 - ω) * n₁ + ω * w₁ * ρ * (1 - 1.5 * u² + 3 * ux + 4.5 * ux^2)
    n₂ = (1 - ω) * n₂ + ω * w₂ * ρ * (1 - 1.5 * u² + 3 * uy + 4.5 * uy^2)
    n₃ = (1 - ω) * n₃ + ω * w₃ * ρ * (1 - 1.5 * u² - 3 * ux + 4.5 * ux^2)
    n₄ = (1 - ω) * n₄ + ω * w₄ * ρ * (1 - 1.5 * u² - 3 * uy + 4.5 * uy^2)
    n₅ = (1 - ω) * n₅ + ω * w₅ * ρ * (1 - 1.5 * u² + 3 * (ux + uy) +
                                      4.5 * (u² + 2 * uxuy))
    n₆ = (1 - ω) * n₆ + ω * w₆ * ρ * (1 - 1.5 * u² + 3 * (uy - ux) +
                                      4.5 * (u² - 2 * uxuy))
    n₇ = (1 - ω) * n₇ + ω * w₇ * ρ * (1 - 1.5 * u² + 3 * (-ux - uy) +
                                      4.5 * (u² + 2 * uxuy))
    n₈ = (1 - ω) * n₈ + ω * w₈ * ρ * (1 - 1.5 * u² + 3 * (ux - uy) +
                                      4.5 * (u² - 2 * uxuy))

    # Force steady inlet flow
    fac = 3 * u₀ + 4.5 * u₀^2 - 1.5 * u₀^2

    n₁[:, 1] = w₁ * (1 + fac)
    n₃[:, 1] = w₃ * (1 - fac)
    n₅[:, 1] = w₅ * (1 + fac)
    n₆[:, 1] = w₆ * (1 - fac)
    n₇[:, 1] = w₇ * (1 - fac)
    n₈[:, 1] = w₈ * (1 + fac)

end

function stream!()
    # Move particles by one step in the direction of motion

    circshift!(n₁, [0, 1])
    circshift!(n₂, [1, 0])
    circshift!(n₃, [0, -1])
    circshift!(n₄, [-1, 0])
    circshift!(n₅, [1, 1])
    circshift!(n₆, [1, -1])
    circshift!(n₇, [-1, -1])
    circshift!(n₈, [-1, 1])

    # Handle barrier collisions
    n₁[b₁] = n₃[barrier]
    n₂[b₂] = n₄[barrier]
    n₃[b₃] = n₁[barrier]
    n₄[b₄] = n₂[barrier]
    n₅[b₅] = n₇[barrier]
    n₆[b₆] = n₈[barrier]
    n₇[b₇] = n₅[barrier]
    n₈[b₈] = n₆[barrier]

end

function curl(ux, uy)
    # Compute curl of macroscopic velocity

    return circshift(uy, [0, -1]) - circshift(uy, [0, 1]) -
           circshift(ux, [-1, 0]) + circshift(ux, [1, 0])

end

# Plot
Pkg.add("PyPlot")
using PyPlot

# Set Plot
fig = PyPlot.figure(figsize=(8,3))
fluidImg = PyPlot.imshow(curl(ux, uy),
                         origin="lower",
                         cmap=PyPlot.get_cmap("jet"),
                         interpolation="none")

# Set RGBA image
rgbaImg =  zeros(4, x₂, x₁)
# Color barrier
rgbaImg[barrier, 3] = 255
# Barrier image
barrierImg = PyPlot.imshow(rgbaImg, origin='lower', interpolation='none')

function nextFrame(frame)
    # Plot next frame

    for step = 1 : 20
        stream()
        collide()
    end

    fluidImg.set_array(curl(ux, uy))

    return fluidImg, barrierImg
end

animate = PyPlot.animation.FuncAnimation(fig, nextFrame, interval=1, blit=true)
PyPlot.show()

max_stars_repo_name: vincentjs/boltzmann

max_issues_repo_name: vincentjs/boltzmann

max_forks_repo_name: vincentjs/boltzmann


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using FEMBasis
using Test

basis = NSeg()
@test size(basis) == (1, 2)
@test length(basis) == 2
N = zeros(1, 2)
eval_basis!(basis, N, (0.0, ))
@test isapprox(N, [0.5 0.5])

basis = NSurf()
@test size(basis) == (2, 4)
@test length(basis) == 4
N = zeros(1, 4)
eval_basis!(basis, N, (0.0, 0.0))
@test isapprox(N, [0.25 0.25 0.25 0.25])

basis = NSolid()
@test size(basis) == (3, 8)
@test length(basis) == 8
N = zeros(1, 8)
eval_basis!(basis, N, (0.0, 0.0, 0.0))
@test isapprox(N, [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125])


text:
module TermMarkdown

using Markdown
import OrderedCollections: OrderedDict
import MyterialColors:
    yellow_light, indigo_light, blue, light_blue, cyan_light, cyan_lighter
import UnicodeFun: to_latex

import Term: reshape_text, highlight_syntax, do_by_line, fillin, escape_brackets
import ..Tables: Table
import ..Style: apply_style
import ..Layout: pad, hLine, vLine
import ..Consoles: console_width
import ..Renderables: RenderableText
import ..Renderables
import ..Tprint: tprint, tprintln
import ..Tprint
import ..Panels: Panel
import ..Segments: Segment

export parse_md

# ---------------------------------------------------------------------------- #
#                                   PARSE MD                                   #
# ---------------------------------------------------------------------------- #
function parse_md(text::Markdown.MD; kwargs...)::String
    elements = parse_md.(text.content; kwargs...)
    return join(elements, "\n\n")
end

function parse_md(header::Markdown.Header{l}; width = console_width(), kwargs...) where {l}
    header_styles = Dict(
        1 => "bold $indigo_light",
        2 => "bold $blue underline",
        3 => "bold $blue",
        4 => "bold $light_blue",
        5 => "bold $cyan_light",
        6 => "bold $cyan_lighter",
    )

    header_justify =
        Dict(1 => :center, 2 => :center, 3 => :center, 4 => :left, 5 => :left, 6 => :left)

    style = header_styles[l]
    header_text = chomp(join(map(ln -> "{$style}$ln{/$style}\n", header.text)))
    if l > 1
        header_text = reshape_text(header_text, width - 2)
        return pad(header_text, width, header_justify[l])
    else
        return string(
            Panel(
                header_text;
                box = :HEAVY,
                style = "dim",
                width = width,
                justify = :center,
                padding = (2, 2, 0, 0),
            ),
        )
    end
end

function parse_md(paragraph::Markdown.Paragraph; width = console_width(), kwargs...)::String
    out = join(parse_md.(paragraph.content; inline = true))
    return reshape_text(out, width)
end

parse_md(italic::Markdown.Italic; kwargs...)::String =
    join(map(ln -> "{italic}$(ln){/italic}", italic.text))
parse_md(bold::Markdown.Bold; kwargs...)::String =
    join(map(ln -> "{bold}$(ln){/bold}", bold.text))
parse_md(lb::Markdown.LineBreak; kwargs...)::String = "\n"

function parse_md(
    ltx::Markdown.LaTeX;
    inline = false,
    width = console_width(),
    kwargs...,
)::String
    formula = ""
    try
        formula = escape_brackets(to_latex(escape_string(string(ltx.formula))))
    catch
        formula = string(ltx.formula)
    end

    if inline
        return "{$yellow_light italic}" * formula * "{/$yellow_light italic}"
    else
        txt = "\n     {$yellow_light italic}" * formula * "{/$yellow_light italic}\n"
        return reshape_text(txt, width) * "\e[0m"
    end
end

function parse_md(
    code::Markdown.Code;
    width = console_width(),
    inline = false,
    kwargs...,
)::String
    syntax = highlight_syntax(code.code)
    if inline
        return "{$yellow_light italic}`$(code.language){/$yellow_light italic}" *
               syntax *
               "{$yellow_light italic}`{/$yellow_light italic}"
    else
        txt = fillin(reshape_text(syntax, width - 20))
        txt = do_by_line(ln -> apply_style(ln, "on_#262626"), txt)
        panel = Panel(
            txt;
            style = "white dim on_#262626",
            box = :SQUARE,
            subtitle = length(code.language) > 0 ? code.language : nothing,
            width = width - 4,
            background = "on_#262626",
            subtitle_justify = :right,
        )

        return string("    " * panel)
    end
end

function parse_md(
    qt::Markdown.BlockQuote;
    width = min(88, console_width() - 30),
    kwargs...,
)::String
    content = parse_md.(qt.content; inline = true)
    content = length(content) > 1 ? join(content) : content[1]
    # content = reshape_text(content, width)

    content = "{bold white}“{/bold white}" * content * "{bold white}”{/bold white}\e[0m"

    content = RenderableText(content; width = max(30, width))
    line =
        content.measure.h > 1 ?
        (
            "{#5a74f2}>{/#5a74f2}" /
            vLine(content.measure.h - 1; style = "#6b77b0 dim", box = :HEAVY)
        ) : "{#5a74f2}>{/#5a74f2}"
    string("  " * line * " " * content)
end

parse_md(hl::Markdown.HorizontalRule; width = console_width(), kwargs...)::String =
    string(hLine(width; style = "dim", box = :HEAVY))

function parse_md(
    list::Markdown.List;
    width = console_width(),
    inline = false,
    space = "",
)::String
    list_elements = []
    for item in list.items
        push!(list_elements, join(parse_md.(item)))
    end

    rendered = ""
    for (i, item) in enumerate(list.items)
        bullet = if Markdown.isordered(list)
            "$space{bold}" * "  $(i + list.ordered - 1). " * "{/bold}"
        else
            "$space{white bold}" * "  • " * "{/white bold}"
        end

        item_content = map(
            elem ->
                elem isa Markdown.List ? "\n" * parse_md(elem; space = "   ") :
                parse_md(elem; inline = true),
            item,
        )
        item_content = length(item_content) > 1 ? join(item_content) : item_content[1]
        rendered *= bullet * item_content * "\n"
    end
    return reshape_text(rendered, width)
end

function parse_md(img::Markdown.Image; kwargs...)::String
    "{dim} 🌄 image: {/dim}" * img.alt * " {dim}| at: " * img.url * "{/dim}"
end

function parse_md(note::Markdown.Footnote; width = console_width(), inline = false)
    if isnothing(note.text)
        return id = "{#9aacdb}[$(note.id)]{/#9aacdb}"
    else
        id = (inline ? "\n" : "") * "{#9aacdb}[$(note.id)]{/#9aacdb}"
        content = parse_md.(note.text)
        return if length(content) == 1
            string(RenderableText("$id: " * content[1]; width = width))
        else
            string(RenderableText("$id:\n" * join(content); width = width))
        end
    end
end

function parse_md(content::Vector; kwargs...)::String
    content = parse_md.(content; kwargs...)
    return if length(content) == 1
        content[1]
    else
        join(content)
    end
end

function parse_md(tb::Markdown.Table; width = console_width())::String
    just = Dict(:l => :left, :r => :right, :c => :center)
    header = parse_md.(tb.rows[1])
    table_content = OrderedDict(
        header[i] => [parse_md(r[i]; inline = true) for r in tb.rows[2:end]] for
        i in 1:length(header)
    )

    return string(
        pad(
            Table(
                table_content;
                columns_justify = [just[j] for j in tb.align],
                box = :ROUNDED,
                header_style = "bold yellow",
                style = "dim",
            );
            width = width,
        ),
    )
end

parse_md(link::Markdown.Link; kwargs...)::String =
    "{white bold}$(parse_md(link.text; inline=true)){/white bold} {dim}($(link.url)){/dim}"

function parse_md(ad::Markdown.Admonition; width = console_width(), kwargs...)::String
    title_styles = Dict(
        "note" => "blue",
        "info" => "blue",
        "warning" => yellow_light,
        "danger" => "red",
        "tip" => "green",
    )
    has_title = length(ad.title) > 0
    content = parse_md(ad.content)
    content = reshape_text(content, width - 10)
    style = get(title_styles, ad.category, "white")
    return string(
        "    " * Panel(
            content;
            title = has_title ? parse_md(ad.title) : "",
            title_style = has_title ? style * " default" : "",
            style = style * " dim",
            width = width - 4,
        ),
    )
end

parse_md(x; kwargs...)::String = string(x)

# ---------------------------------------------------------------------------- #
#                              RENDERABLE & TPRINT                             #
# ---------------------------------------------------------------------------- #

Renderables.RenderableText(md::Markdown.MD; width = console_width() - 2, kwargs...) =
    RenderableText(parse_md(md; width = width); width = width, kwargs...)

Tprint.tprint(md::Markdown.MD; kwargs...) = tprint(parse_md(md); kwargs...)
Tprint.tprint(io::IO, md::Markdown.MD; kwargs...) = tprint(io, parse_md(md); kwargs...)

Tprint.tprintln(md::Markdown.MD; kwargs...) = tprintln(parse_md(md); kwargs...)
Tprint.tprintln(io::IO, md::Markdown.MD; kwargs...) = tprintln(io, parse_md(md); kwargs...)
end


text:
@discord_object struct VoiceState
    guild_id::Snowflake
    channel_id::Snowflake
    user_id::Snowflake
    member::GuildMember
    session_id::String
    deaf::Bool
    mute::Bool
    self_deaf::Bool
    self_mute::Bool
    self_stream::Bool
    suppress::Bool
end


text:
using Test; using Laplacians; using MAT; using LinearAlgebra; using BenchmarkTools
include("../src/CombinatorialMultigrid.jl");

using Main.CombinatorialMultigrid

## load example matrix

file = matopen("../example/X.mat"); X = read(file, "X"); close(file)
LX = lap(X);
b1 = rand(Float64, size(X, 1));
b1 = b1 .- sum(b1)/length(b1);

@info "###  Running CMG ###"
t = @elapsed (pfunc, h) = cmg_preconditioner_lap(LX);
@info "Time Required to build CMG Solver: $(t) seconds"
t = @elapsed x = pfunc(b1);
@info "Time Required to find x: $(t) seconds"

@info "###  Running PCG with CMG    ### "
## solve with pcg and cmg preconditioner
f1 = pcgSolver(LX,pfunc);
t = @elapsed x = f1(b1, maxits=40, tol=1e-6,verbose=true);
@info "Time Required to solve system: $(t) seconds"

## solve with approxchol_lap from laplacians
@info "###  Running approxchol_lap    ###"
t = @elapsed solver = approxchol_lap(X; tol=1e-6, maxits=1000, maxtime=Inf, verbose=true, pcgIts=Int[], params=ApproxCholParams());
@info "Time Required to build Lap Solver: $(t) seconds"
t = @elapsed x = solver(b1);
@info "Time Required to find x: $(t) seconds"


text:
import LinearAlgebra: I, Hermitian, tr, pinv
import ForneyLab: unsafeCov, unsafeMean, unsafePrecision, VariateType

export ruleVariationalLARXOutNPPPPP,
       ruleVariationalLARXIn1PNPPPP,
       ruleVariationalLARXIn2PPNPPP,
       ruleVariationalLARXIn3PPPNPP,
	   ruleVariationalLARXIn4PPPPNP,
	   ruleVariationalLARXIn5PPPPPN

# Autoregression order and bookkeeping matrices
order = Nothing
S = Nothing
s = Nothing

function defineOrder(dim::Int64)
    global order, s, S

	# Set autoregression order
    order = dim

	# Set bookkeeping matrices
    s = shift_vec(order)
    S = shift_mat(order)
end

function ruleVariationalLARXOutNPPPPP(marg_y :: Nothing,
                                      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_θ :: ProbabilityDistribution{Multivariate},
                                      marg_η :: ProbabilityDistribution{Univariate},
                                      marg_u :: ProbabilityDistribution{Univariate},
                                      marg_γ :: ProbabilityDistribution{Univariate})

    # Expectations of incoming marginal beliefs
    mz = unsafeMean(marg_z)
    mθ = unsafeMean(marg_θ)
    mη = unsafeMean(marg_η)
    mu = unsafeMean(marg_u)
    mγ = unsafeMean(marg_γ)

    # Check order
	if order == Nothing; defineOrder(length(mz)); end

		# Construct precision matrix
    Wmγ = AR_precisionmat(mγ, order)

    # Parameters of outgoing message
	Φ = Wmγ
    ϕ = (S + s*mθ')*mz + s*mη*mu

	# Set outgoing message
	return Message(Multivariate, GaussianMeanPrecision, m=ϕ, w=Φ)
end

function ruleVariationalLARXIn1PNPPPP(marg_y :: ProbabilityDistribution{Multivariate},
                                      marg_z :: Nothing,
                                      marg_θ :: ProbabilityDistribution{Multivariate},
                                      marg_η :: ProbabilityDistribution{Univariate},
                                      marg_u :: ProbabilityDistribution{Univariate},
                                      marg_γ :: ProbabilityDistribution{Univariate})

    # Expectations of marginal beliefs
    my = unsafeMean(marg_y)
    mθ,Vθ = unsafeMeanCov(marg_θ)
    mη = unsafeMean(marg_η)
    mu = unsafeMean(marg_u)
    mγ = unsafeMean(marg_γ)

	# Check order
	if order == Nothing; defineOrder(length(mz)); end

    # Construct precision matrix
	Wmγ = AR_precisionmat(mγ, order)

	# Parameters of outgoing message
    Φ = S'*Wmγ*S + mγ*(mθ*mθ'+ Vθ)
    ϕ = (S + s*mθ')'*Wmγ*(my - s*mη*mu)

    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalLARXIn2PPNPPP(marg_y :: ProbabilityDistribution{Multivariate},
                                      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_θ :: Nothing,
							  	      marg_η :: ProbabilityDistribution{Univariate},
								      marg_u :: ProbabilityDistribution{Univariate},
                                      marg_γ :: ProbabilityDistribution{Univariate})

	# Expectations of marginal beliefs
	my = unsafeMean(marg_y)
	mz,Vz = unsafeMeanCov(marg_z)
	mη = unsafeMean(marg_η)
	mu = unsafeMean(marg_u)
	mγ = unsafeMean(marg_γ)

	# Check order
	if order == Nothing; defineOrder(length(mz)); end

    # Construct precision matrix
	Wmγ = AR_precisionmat(mγ, order)

    # Parameters of outgoing message
    Φ = mγ*(mz*mz' + Vz)
    ϕ = mz*s'*Wmγ*(my - s*mη*mu)

    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalLARXIn3PPPNPP(marg_y :: ProbabilityDistribution{Multivariate},
                                      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_θ :: ProbabilityDistribution{Multivariate},
								      marg_η :: Nothing,
								      marg_u :: ProbabilityDistribution{Univariate},
                                      marg_γ :: ProbabilityDistribution{Univariate})

 	# Expectations of marginal beliefs
	mθ = unsafeMean(marg_θ)
    my = unsafeMean(marg_y)
    mz = unsafeMean(marg_z)
	mu,vu = unsafeMeanCov(marg_u)
	mγ = unsafeMean(marg_γ)

	# Check order
	if order == Nothing; defineOrder(length(mz)); end

    # Construct precision matrix
    Wmγ = AR_precisionmat(mγ, order)

	# Parameters of outgoing message
	Φ = mγ*(mu^2 + vu)
    ϕ = mu*s'*Wmγ*(my - (S + s*mθ')*mz)

	return Message(Univariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalLARXIn4PPPPNP(marg_y :: ProbabilityDistribution{Multivariate},
                                      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_θ :: ProbabilityDistribution{Multivariate},
							  	      marg_η :: ProbabilityDistribution{Univariate},
								      marg_u :: Nothing,
                                      marg_γ :: ProbabilityDistribution{Univariate})

 	# Expectations of marginal beliefs
	mθ = unsafeMean(marg_θ)
    my = unsafeMean(marg_y)
    mz = unsafeMean(marg_z)
	mη,vη = unsafeMeanCov(marg_η)
	mγ = unsafeMean(marg_γ)

	# Check order
	if order == Nothing; defineOrder(length(mz)); end

    # Construct precision matrix
    Wmγ = AR_precisionmat(mγ, order)

	# Parameters of outgoing message
	Φ = mγ*(mη^2 + vη)
    ϕ = mη*s'*Wmγ*(my - (S + s*mθ')*mz)

	return Message(Univariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalLARXIn5PPPPPN(marg_y :: ProbabilityDistribution{Multivariate},
                                      marg_z :: ProbabilityDistribution{Multivariate},
                                      marg_θ :: ProbabilityDistribution{Multivariate},
							  	      marg_η :: ProbabilityDistribution{Univariate},
								      marg_u :: ProbabilityDistribution{Univariate},
                                      marg_γ :: Nothing)

    # Expectations of marginal beliefs
	mθ,Vθ = unsafeMeanCov(marg_θ)
    my,Vy = unsafeMeanCov(marg_y)
    mz,Vz = unsafeMeanCov(marg_z)
	mη,vη = unsafeMeanCov(marg_η)
	mu,vu = unsafeMeanCov(marg_u)

	# Check order
	if order == Nothing; defineOrder(length(mz)); end

	# Convenience variables
	EA = (S + s*mθ')
	EB = s*mη

	# Intermediate terms
	term1 = (my*my' + Vy)[1,1]
	term2 = -(EA*mz*my')[1,1]
	term3 = -(EB*mu*my')[1,1]
	term4 = -(my*mz'*EA)[1,1]
	term5 = ((S'*S + mθ*mθ' + Vθ)*(mz*mz' + Vz))[1,1]
	term6 = (EB*mu*mz'*EA')[1,1]
	term7 = (my*mu*EB')[1,1]
	term8 = (EA*mz*mu*EB')[1,1]
	term9 = (mu^2*s*(mη*mη' + vη)*s')[1,1]

	# print("term 1 = ")
	# println(term1)
	# print("term 2 = ")
	# println(term2)
	# print("term 3 = ")
	# println(term3)
	# print("term 4 = ")
	# println(term4)
	# print("term 5 = ")
	# println(term5)
	# print("term 6 = ")
	# println(term6)
	# print("term 7 = ")
	# println(term7)
	# print("term 8 = ")
	# println(term8)
	# print("term 9 = ")
	# println(term9)

	# Parameters of outgoing message
	a = 3/2
    B = 1/2*(term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8 + term9)

	return Message(Gamma, a=a, b=B)
end

max_stars_repo_name: wmkouw/LARX

max_issues_repo_name: wmkouw/LARX

max_forks_repo_name: wmkouw/LARX


text:
module RectifierSimulation

using Modia
@usingModiaPlot

setLogMerge(false)

include("$(Modia.modelsPath)/Electric.jl")

# Sinusoidal voltage source
SineVoltage = OnePort | Model( V = 1.0u"V", f = 1.0u"Hz", equations = :[ v = V*sin(2*3.14*f*time) ] )

# Ideal diode
IdealDiode = OnePort | Model(
        Ron   = 1e-5u"Ω",
        Goff  = 1e-5u"1/Ω",
        s = Var(start = 0.0),
        equations = :[
            closed = positive(s)   # closed = s > 0
            v = s*u"V"*(if closed; Ron*u"1/Ω" else 1         end)
            i = s*u"A"*(if closed; 1          else Goff*u"Ω" end)
        ]
    )

Rectifier1 = Model(
    R1 = Resistor | Map(R=1.0u"Ω"),
    R2 = Resistor | Map(R=100.0u"Ω"),
    C = Capacitor | Map(C=0.1u"F", v=Var(init=0.0u"V")),
    D = IdealDiode,
    V = SineVoltage | Map(V=5.0u"V", f=1.5u"Hz"),
    ground = Ground,
    connect = :[
      (V.p , R1.p)
      (R1.n, D.p)
      (D.n , R2.p, C.p)
      (ground.p, R2.n, C.n, V.n)
    ]
)

rectifier1 = @instantiateModel(Rectifier1)
@time simulate!(rectifier1, Tsit5(), stopTime = 3, requiredFinalStates=[4.5665505086780565])
plot(rectifier1, [("V.v", "D.v", "C.v"), "D.i"], figure=1)

Rectifier2 = Model(
    V  = SineVoltage | Map(V=220.0u"V", f=50.0u"Hz"),
    R1 = Resistor    | Map(R=20.0u"Ω"),
    R2 = Resistor    | Map(R=500.0u"Ω"),
    C  = Capacitor   | Map(C=1e-4u"F", v=Var(init=0.0u"V")),
    D1 = IdealDiode,
    D2 = IdealDiode,
    D3 = IdealDiode,
    D4 = IdealDiode,
    ground = Ground,
    connect = :[
      (V.p , R1.p)
      (R1.n, D1.p, D3.n)
      (V.n , D2.p, D4.n, ground.p)
      (D1.n, D2.n, R2.p, C.p)
      (D3.p, D4.p, R2.n, C.n)
    ]
)

rectifier2 = @instantiateModel(Rectifier2, unitless=true)
@time simulate!(rectifier2, Tsit5(), stopTime = 0.1, log=true, requiredFinalStates=[183.9899542497182])
plot(rectifier2, [("V.v", "C.v"), "V.i"], figure=2)

@time simulate!(rectifier2, stopTime = 0.1, log=true)  #, requiredFinalStates=[183.9899542497182])
plot(rectifier2, [("V.v", "C.v"), "V.i"], figure=3)

end