text:
# static calculations used by the full model (static and dynamic)
function L_tilde(g, z_hat, Ω, E, S, parameters)
    @unpack N, θ, κ, χ, ζ = parameters
    return Ω * ((N-1) * z_hat^(-θ)*κ + ζ*(S + E / χ)) # (33) (nests the stationary case with E = δ)
end

function L_tilde_x(z_hat, Ω, parameters)
    @unpack N, θ, κ = parameters
    return Ω * (N - 1) * z_hat^(-θ) * κ # (34)
end

L_tilde_E(Ω, E, parameters) = Ω * parameters.ζ * E / parameters.χ # (35) (nests the stationary case with E = δ)
L_tilde_a(Ω, S, parameters) = Ω * parameters.ζ * S # (36)

function S(g, parameters)
    @unpack θ, μ, υ = parameters
    return θ * (g - μ - θ * υ^2 /2) # (32)
end

w(z_bar, parameters) = (parameters.σ - 1)/(parameters.σ) * z_bar # (C.13)

function x(w, parameters)
    @unpack Theta, ζ, η = parameters
    return ζ * (1 - η + η * Theta / w) # (C.14)
end

function λ_ii(z_hat, parameters)
    @unpack N, σ, θ, d = parameters
    return  1/(1 + (N-1)*z_hat^(σ-1-θ)*d^(1-σ)) # (C.79)
end

c(L_tilde, Ω, z_bar) = (1 - L_tilde)*z_bar # (C.15)
function entry_residual(v_1, Ξ₁, parameters)
    @assert parameters.η ≈ 0.
    Ξ₁*v_1 - parameters.ζ*(1- parameters.χ)/parameters.χ # (56)
end
π_min(L_tilde, z_bar, parameters) = (1 - L_tilde) / ((parameters.σ-1)*z_bar^(parameters.σ-1)) # (38)

function π_rat(z_hat, parameters)
    @unpack θ, σ, d, κ, N = parameters
    return θ/(1 + θ - σ) + (N - 1) * d^(1-σ) * (σ - 1)*z_hat^(σ - 1 - θ)/(1 + θ - σ) # (C.17)
end

function z_bar(z_hat, Ω, parameters)
    @unpack θ, σ, N, d = parameters
    return (Ω*θ/(1 + θ - σ) * (1 + (N - 1)*d^(1 - σ)*z_hat^(σ - 1 - θ)))^(1/(σ-1)) # (C.11)
end

# subset of the above used in the dynamic model
function static_equilibrium(Ξ₁, v_1, g, z_hat, E, Ω, z, parameters)
    @unpack θ, σ, N, d, χ, ζ, κ = parameters
    S_t = S(g, parameters)
    L_tilde_t = L_tilde(g, z_hat, Ω, E, S_t, parameters)
    z_bar_t = z_bar(z_hat, Ω, parameters)
    w_t = w(z_bar_t, parameters)
    π_min_t = (1 - L_tilde_t) / ((σ-1)*z_bar_t^(σ-1))  # (C.20)

    i_vectorized = z .>= log(z_hat) # Vectorized indicator function
    π_t = π_min_t * (1.0.+(N-1)*d^(1-σ)*i_vectorized) - (N-1)*κ*exp.(-(σ-1).*z).*i_vectorized  # (39)
    entry_residual_t = entry_residual(v_1, Ξ₁, parameters)

    return (S = S_t, L_tilde = L_tilde_t, z_bar = z_bar_t, π_min = π_min_t, π = π_t, entry_residual = entry_residual_t,
            w = w_t)
end

consumption_equivalent(U, U_old, parameters) = isapprox(parameters.γ, 1.0) ?  exp(parameters.ρ*(U - U_old)) - 1 : (U/U_old)^(1.0/(1.0 - parameters.γ)) -1# H.23 main paper
ACR(new_λ, old_λ, params) = (-1/params.θ)*log(new_λ/old_λ)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestSchur

using Test, LinearAlgebra, Random
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

Random.seed!(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2

@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int)
    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    asym = a' + a                 # symmetric indefinite
    apd  = a' * a                 # symmetric positive-definite
    for (a, asym, apd) in ((a, asym, apd),
                           (view(a, 1:n, 1:n),
                            view(asym, 1:n, 1:n),
                            view(apd, 1:n, 1:n)))
        ε = εa = eps(abs(float(one(eltya))))

        d,v = eigen(a)
        f   = schur(a)
        @test f.vectors*f.Schur*f.vectors' ≈ a
        @test sort(real(f.values)) ≈ sort(real(d))
        @test sort(imag(f.values)) ≈ sort(imag(d))
        @test istriu(f.Schur) || eltype(a)<:Real
        @test convert(Array, f) ≈ a
        @test_throws ErrorException f.A

        sch, vecs, vals = schur(UpperTriangular(triu(a)))
        @test vecs*sch*vecs' ≈ triu(a)
        sch, vecs, vals = schur(UnitUpperTriangular(triu(a)))
        @test vecs*sch*vecs' ≈ UnitUpperTriangular(triu(a))
        sch, vecs, vals = schur(LowerTriangular(tril(a)))
        @test vecs*sch*vecs' ≈ tril(a)
        sch, vecs, vals = schur(UnitLowerTriangular(tril(a)))
        @test vecs*sch*vecs' ≈ UnitLowerTriangular(tril(a))
        sch, vecs, vals = schur(Hermitian(asym))
        @test vecs*sch*vecs' ≈ asym
        sch, vecs, vals = schur(Symmetric(a + transpose(a)))
        @test vecs*sch*vecs' ≈ a + transpose(a)
        sch, vecs, vals = schur(Tridiagonal(a + transpose(a)))
        @test vecs*sch*vecs' ≈ Tridiagonal(a + transpose(a))
        sch, vecs, vals = schur(Bidiagonal(a, :U))
        @test vecs*sch*vecs' ≈ Bidiagonal(a, :U)
        sch, vecs, vals = schur(Bidiagonal(a, :L))
        @test vecs*sch*vecs' ≈ Bidiagonal(a, :L)

        tstring = sprint((t, s) -> show(t, "text/plain", s), f.T)
        zstring = sprint((t, s) -> show(t, "text/plain", s), f.Z)
        vstring = sprint((t, s) -> show(t, "text/plain", s), f.values)
        fstring = sprint((t, s) -> show(t, "text/plain", s), f)
        @test fstring == "$(summary(f))\nT factor:\n$tstring\nZ factor:\n$(zstring)\neigenvalues:\n$vstring"
        @testset "Reorder Schur" begin
            # use asym for real schur to enforce tridiag structure
            # avoiding partly selection of conj. eigenvalues
            ordschura = eltya <: Complex ? a : asym
            S = schur(ordschura)
            select = bitrand(n)
            O = ordschur(S, select)
            sum(select) != 0 && @test S.values[findall(select)] ≈ O.values[1:sum(select)]
            @test O.vectors*O.Schur*O.vectors' ≈ ordschura
            @test_throws ErrorException f.A
            Snew = LinearAlgebra.Schur(S.T, S.Z, S.values)
            SchurNew = ordschur!(copy(Snew), select)
            @test O.vectors ≈ SchurNew.vectors
            @test O.Schur ≈ SchurNew.Schur
        end

        if isa(a, Array)
            a1_sf = a[1:n1, 1:n1]
            a2_sf = a[n1+1:n2, n1+1:n2]
        else
            a1_sf = view(a, 1:n1, 1:n1)
            a2_sf = view(a, n1+1:n2, n1+1:n2)
        end
        @testset "Generalized Schur" begin
            f = schur(a1_sf, a2_sf)
            @test f.Q*f.S*f.Z' ≈ a1_sf
            @test f.Q*f.T*f.Z' ≈ a2_sf
            @test istriu(f.S) || eltype(a)<:Real
            @test istriu(f.T) || eltype(a)<:Real
            @test_throws ErrorException f.A

            sstring = sprint((t, s) -> show(t, "text/plain", s), f.S)
            tstring = sprint((t, s) -> show(t, "text/plain", s), f.T)
            qstring = sprint((t, s) -> show(t, "text/plain", s), f.Q)
            zstring = sprint((t, s) -> show(t, "text/plain", s), f.Z)
            αstring = sprint((t, s) -> show(t, "text/plain", s), f.α)
            βstring = sprint((t, s) -> show(t, "text/plain", s), f.β)
            fstring = sprint((t, s) -> show(t, "text/plain", s), f)
            @test fstring == "$(summary(f))\nS factor:\n$sstring\nT factor:\n$(tstring)\nQ factor:\n$(qstring)\nZ factor:\n$(zstring)\nα:\n$αstring\nβ:\n$βstring"
        end
        @testset "Reorder Generalized Schur" begin
            NS = schur(a1_sf, a2_sf)
            # Currently just testing with selecting gen eig values < 1
            select = abs2.(NS.values) .< 1
            m = sum(select)
            S = ordschur(NS, select)
            # Make sure that the new factorization still factors matrix
            @test S.Q*S.S*S.Z' ≈ a1_sf
            @test S.Q*S.T*S.Z' ≈ a2_sf
            # Make sure that we have sorted it correctly
            @test NS.values[findall(select)] ≈ S.values[1:m]

            Snew = LinearAlgebra.GeneralizedSchur(NS.S, NS.T, NS.alpha, NS.beta, NS.Q, NS.Z)
            SchurNew = ordschur!(copy(Snew), select)
            @test S.Q ≈ SchurNew.Q
            @test S.S ≈ SchurNew.S
            @test S.T ≈ SchurNew.T
            @test S.Z ≈ SchurNew.Z
            @test S.alpha ≈ SchurNew.alpha
            @test S.beta  ≈ SchurNew.beta
            sS,sT,sQ,sZ = schur(a1_sf,a2_sf)
            @test NS.Q ≈ sQ
            @test NS.T ≈ sT
            @test NS.S ≈ sS
            @test NS.Z ≈ sZ
        end
    end
    @testset "0x0 matrix" for A in (zeros(eltya, 0, 0), view(rand(eltya, 2, 2), 1:0, 1:0))
        T, Z, λ = LinearAlgebra.schur(A)
        @test T == A
        @test Z == A
        @test λ == zeros(0)
    end

    if eltya <: Real
        @testset "quasitriangular to triangular" begin
            S = schur(a)
            SC = Schur{Complex}(S)
            @test eltype(SC) == complex(eltype(S))
            @test istriu(SC.T)
            @test SC.Z*SC.Z' ≈ I
            @test SC.Z*SC.T*SC.Z' ≈ a
            @test sort(SC.values,by=LinearAlgebra.eigsortby) ≈ sort(S.values,by=LinearAlgebra.eigsortby)
            @test Schur{Complex}(SC) === SC === Schur{eltype(SC)}(SC)
            @test Schur{eltype(S)}(S) === S
            if eltype(S) === Float32
                S64 = Schur{Float64}(S)
                @test eltype(S64) == Float64
                @test S64.Z == S.Z
                @test S64.T == S.T
                @test S64.values == S.values
            end
        end
    end

    @testset "0x0 $eltya matrices" begin
        A = zeros(eltya, 0, 0)
        B = zeros(eltya, 0, 0)
        S = LinearAlgebra.schur(A, B)
        @test S.S == A
        @test S.T == A
        @test S.Q == A
        @test S.Z == A
        @test S.alpha == zeros(0)
        @test S.beta == zeros(0)
    end
end

@testset "Generalized Schur convergence" begin
    # Check for convergence issues, #40279
    problematic_pencils = [
        (   ComplexF64[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 3.7796350217469814 -3.3125635598133054 0.0 0.0 0.0 0.0 0.0 0.0 6.418270043493963 -6.625127119626611 0.0 0.0 0.0 0.0 0.0 -1.0; -3.312563559813306 3.779635021746982 0.0 0.0 0.0 0.0 0.0 0.0 -6.625127119626612 6.418270043493964 -1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 3.7796350217469814 0.0 0.0 -3.3125635598133054 0.0 0.0 0.0 -1.0 6.418270043493963 0.0 0.0 -6.625127119626611 0.0 0.0; 0.0 0.0 0.0 3.779635021746982 -3.312563559813306 0.0 0.0 0.0 0.0 0.0 0.0 6.418270043493964 -6.625127119626612 0.0 -1.0 0.0; 0.0 0.0 0.0 -3.3125635598133054 3.7796350217469814 0.0 0.0 0.0 0.0 0.0 0.0 -6.625127119626611 6.418270043493963 -1.0 0.0 0.0; 0.0 0.0 -3.312563559813306 0.0 0.0 3.779635021746982 0.0 0.0 0.0 0.0 -6.625127119626612 0.0 -1.0 6.418270043493964 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 3.7796350217469814 -3.3125635598133054 0.0 0.0 0.0 -1.0 0.0 0.0 6.418270043493963 -6.625127119626611; 0.0 0.0 0.0 0.0 0.0 0.0 -3.312563559813306 3.779635021746982 -1.0 0.0 0.0 0.0 0.0 0.0 -6.625127119626612 6.418270043493964],
            ComplexF64[1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.7796350217469814 3.312563559813306 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.3125635598133054 -3.779635021746982 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.7796350217469814 0.0 0.0 3.312563559813306 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.779635021746982 3.3125635598133054 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.312563559813306 -3.7796350217469814 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.3125635598133054 0.0 0.0 -3.779635021746982 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.7796350217469814 3.312563559813306; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.3125635598133054 -3.779635021746982]
        ),
        (   ComplexF64[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.62 -1.0 0.0 0.0 0.0 0.0 -1.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 -1.0 -2.62 0.0 0.0 0.0 0.0 0.0; 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0 0.0; 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62 0.0; 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 -2.62],
            ComplexF64[1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
        ),
        (   ComplexF64[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 0.33748484079831426 -0.10323794456968927 0.0 0.0 0.0 0.0 0.0 0.0 -2.5940303184033713 -0.20647588913937853 0.0 0.0 0.0 0.0 0.0 -1.0; -0.10323794456968927 0.3374848407983142 0.0 0.0 0.0 0.0 0.0 0.0 -0.20647588913937853 -2.5940303184033713 -1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.33748484079831426 0.0 0.0 -0.10323794456968927 0.0 0.0 0.0 -1.0 -2.5940303184033713 0.0 0.0 -0.20647588913937853 0.0 0.0; 0.0 0.0 0.0 0.3374848407983142 -0.10323794456968927 0.0 0.0 0.0 0.0 0.0 0.0 -2.5940303184033713 -0.20647588913937853 0.0 -1.0 0.0; 0.0 0.0 0.0 -0.10323794456968927 0.33748484079831426 0.0 0.0 0.0 0.0 0.0 0.0 -0.20647588913937853 -2.5940303184033713 -1.0 0.0 0.0; 0.0 0.0 -0.10323794456968927 0.0 0.0 0.3374848407983142 0.0 0.0 0.0 0.0 -0.20647588913937853 0.0 -1.0 -2.5940303184033713 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.33748484079831426 -0.10323794456968927 0.0 0.0 0.0 -1.0 0.0 0.0 -2.5940303184033713 -0.20647588913937853; 0.0 0.0 0.0 0.0 0.0 0.0 -0.10323794456968927 0.3374848407983142 -1.0 0.0 0.0 0.0 0.0 0.0 -0.20647588913937853 -2.5940303184033713],
            ComplexF64[1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.33748484079831426 0.10323794456968927 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.10323794456968927 -0.3374848407983142 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.33748484079831426 0.0 0.0 0.10323794456968927 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.3374848407983142 0.10323794456968927 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.10323794456968927 -0.33748484079831426 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.10323794456968927 0.0 0.0 -0.3374848407983142 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.33748484079831426 0.10323794456968927; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.10323794456968927 -0.3374848407983142]
        ),
        (   ComplexF64[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 1.7391668762048442 -1.309613611600033 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.150333752409688 -2.619227223200066 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0; -1.3096136116000332 1.739166876204844 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.6192272232000664 2.150333752409688 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.739166876204844 0.0 0.0 -1.3096136116000332 0.0 0.0 0.0 0.0 0.0 -1.0 2.150333752409688 0.0 0.0 -2.6192272232000664 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.739166876204844 0.0 0.0 0.0 0.0 -1.3096136116000332 0.0 -1.0 0.0 0.0 2.150333752409688 0.0 0.0 0.0 0.0 -2.6192272232000664 0.0; 0.0 0.0 0.0 0.0 1.7391668762048442 0.0 0.0 0.0 0.0 -1.309613611600033 0.0 0.0 0.0 0.0 2.150333752409688 -1.0 0.0 0.0 0.0 -2.619227223200066; 0.0 0.0 -1.309613611600033 0.0 0.0 1.7391668762048442 0.0 0.0 0.0 0.0 0.0 0.0 -2.619227223200066 0.0 -1.0 2.150333752409688 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.739166876204844 -1.3096136116000332 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.150333752409688 -2.6192272232000664 0.0 -1.0; 0.0 0.0 0.0 0.0 0.0 0.0 -1.309613611600033 1.7391668762048442 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.619227223200066 2.150333752409688 -1.0 0.0; 0.0 0.0 0.0 -1.309613611600033 0.0 0.0 0.0 0.0 1.7391668762048442 0.0 0.0 0.0 0.0 -2.619227223200066 0.0 0.0 0.0 -1.0 2.150333752409688 0.0; 0.0 0.0 0.0 0.0 -1.3096136116000332 0.0 0.0 0.0 0.0 1.739166876204844 0.0 0.0 0.0 0.0 -2.6192272232000664 0.0 -1.0 0.0 0.0 2.150333752409688],
            ComplexF64[1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.7391668762048442 1.3096136116000332 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.309613611600033 -1.739166876204844 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.739166876204844 0.0 0.0 1.309613611600033 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.739166876204844 0.0 0.0 0.0 0.0 1.309613611600033 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.7391668762048442 0.0 0.0 0.0 0.0 1.3096136116000332; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.3096136116000332 0.0 0.0 -1.7391668762048442 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.739166876204844 1.309613611600033 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.3096136116000332 -1.7391668762048442 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.3096136116000332 0.0 0.0 0.0 0.0 -1.7391668762048442 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.309613611600033 0.0 0.0 0.0 0.0 -1.739166876204844]
        ),
        (   ComplexF64[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230788 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007; 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230788 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230788 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769246 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230784 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230788 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230788 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230788 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769246 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230784 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230788 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.90076923076925 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0000000000000007 -12.019230769230788; -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 -6.009615384615393 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384622 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769244 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615393 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384622 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769244 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.490384615384624 -1.0000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -12.019230769230784 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.900769230769248],
            ComplexF64[1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615393 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615393 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384622 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615392 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384622 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.009615384615394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.490384615384624]
        )]

    for (A, B) in problematic_pencils
        f = schur(A, B)
        @test f.Q*f.S*f.Z' ≈ A
        @test f.Q*f.T*f.Z' ≈ B
    end
end

@testset "adjoint and transpose for schur (#40941)" begin
    A = rand(3, 3)
    B = schur(A', A)
    C = B.left*B.S*B.right'
    D = schur(transpose(A), A)
    E = D.left*D.S*D.right'
    @test A' ≈ C ≈ E
end

end # module TestSchur

ext: jl

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
using Unitful
import Unitful: ms, s, kHz, Hz
export dB, ms, s, kHz, Hz, frames, uconvert, ustrip, inseconds, inHz, inframes

const TimeDim = Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)}
const FreqDim = Unitful.Dimensions{(Unitful.Dimension{:Time}(-1//1),)}
const Time{N} = Quantity{N,TimeDim}
const Freq{N} = Quantity{N,FreqDim}

@dimension 𝐅r "Fr" Frame
@refunit frames "frames" Frames 𝐅r false
const FrameDim = Unitful.Dimensions{(Unitful.Dimension{:Frame}(1//1),)}
const FrameQuant{N} = Quantity{N,FrameDim}

"""
    inframes(quantity,rate)

Translate the given quantity (usually a time) to a sample index, given
a particualr samplerate.

# Example

> inframes(1s,44100Hz)
44100

"""
inframes{N <: Integer}(time::FrameQuant{N},rate) = ustrip(time)
inframes{N}(time::FrameQuant{N},rate) =
  error("Cannot use non-integer frames.")
function inframes(time,rate)
  r = inHz(rate)
  floor(Int,ustrip(inseconds(time,r)*r))
end
function inframes{N,M}(time::Time{N},rate::Freq{M})
  floor(Int,ustrip(uconvert(s,time)*uconvert(Hz,rate)))
end

"""
    inHz(quantity)

Translate a particular quantity (usually a frequency) to a value in Hz.

# Example

> inHz(1.0kHz)
1000.0 Hz

"""
inHz(x::Quantity) = uconvert(Hz,x)
inHz(typ::Type{N},x::Q) where {N <: Number,Q <: Quantity} =
  floor(N,ustrip(inHz(x)))*Hz
inHz(typ::Type{N},x::N) where {N <: Number} = inHz(x)

"""
   inseconds(quantity)

Translate a particular quantity (usually a time) to a value in seconds.

# Example
> inseconds(50.0ms)
0.05 s

"""
inseconds(x::FrameQuant,R) = (ustrip(x) / R)*s
inseconds(x::Time) = uconvert(s,x)
inseconds(x::Quantity,R) = uconvert(s,x)
inseconds(x::Number,R) = inseconds(x)
inseconds(x::FrameQuant) =
  error("Expected second argument, specifying sample rate.")


text:
include("utils.jl")

@testset "Stiefel" begin
    @testset "Real" begin
        M = Stiefel(3,2)
        @testset "Basics" begin
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            @test representation_size(M) == (3,2)
            @test manifold_dimension(M) == 3
            @test_throws DomainError is_manifold_point(M, [1., 0., 0., 0.],true)
            @test_throws DomainError is_manifold_point(M, 1im*[1.0 0.0; 0.0 1.0; 0.0 0.0],true)
            @test !is_tangent_vector(M, x, [0., 0., 1., 0.])
            @test_throws DomainError is_tangent_vector(M, x, 1 * im * zero_tangent_vector(M,x), true)
        end

        types = [
            Matrix{Float64},
            MMatrix{3, 2, Float64},
            Matrix{Float32},
        ]
        @testset "Type $T" for T in types
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            y = exp(M,x, [0.0 0.0; 0.0 0.0; 1.0 1.0])
            z = exp(M,x,[ 0.0 0.0; 0.0 0.0; -1.0 1.0])
            pts = convert.(T, [x,y,z])
            v = inverse_retract(M,x,y,PolarInverseRetraction())
            @test !is_manifold_point(M,2*x)
            @test_throws DomainError !is_manifold_point(M,2*x,true)
            @test !is_tangent_vector(M,2*x,v)
            @test_throws DomainError !is_tangent_vector(M,2*x,v,true)
            @test !is_tangent_vector(M,x,y)
            @test_throws DomainError is_tangent_vector(M,x,y,true)
            test_manifold(
                M,
                pts,
                test_exp_log = false,
                default_inverse_retraction_method = PolarInverseRetraction(),
                test_injectivity_radius = false,
                test_is_tangent = true,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                projection_atol_multiplier = 15.0,
                retraction_atol_multiplier = 10.0,
                is_tangent_atol_multiplier = 4*10.0^2,
                retraction_methods = [PolarRetraction(), QRRetraction()],
                inverse_retraction_methods = [PolarInverseRetraction(), QRInverseRetraction()]
            )

            @testset "inner/norm" begin
                v1 = inverse_retract(M, pts[1], pts[2], PolarInverseRetraction())
                v2 = inverse_retract(M, pts[1], pts[3], PolarInverseRetraction())

                @test real(inner(M, pts[1], v1, v2)) ≈ real(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v2)) ≈ -imag(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v1)) ≈ 0

                @test norm(M, pts[1], v1) isa Real
                @test norm(M, pts[1], v1) ≈ sqrt(inner(M, pts[1], v1, v1))
            end
        end
    end

    @testset "Complex" begin
        M = Stiefel(3,2,ℂ)
        @testset "Basics" begin
            @test representation_size(M) == (3,2)
            @test manifold_dimension(M) == 8
            @test !is_manifold_point(M, [1., 0., 0., 0.])
            @test !is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], [0., 0., 1., 0.])
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            @test_throws DomainError is_manifold_point(M, [:a :b; :c :d; :e :f],true)
            @test_throws DomainError is_tangent_vector(M, x, [:a :b; :c :d; :e :f], true)
        end
        types = [
            Matrix{ComplexF64},
        ]
        @testset "Type $T" for T in types
            x = [0.5+0.5im 0.5+0.5im; 0.5+0.5im -0.5-0.5im; 0.0 0.0]
            y = exp(M, x, [0.0 0.0; 0.0 0.0; 1.0 1.0])
            z = exp(M, x, [0.0 0.0; 0.0 0.0; -1.0 1.0])
            pts = convert.(T, [x,y,z])
            v = inverse_retract(M,x,y,PolarInverseRetraction())
            @test !is_manifold_point(M,2*x)
            @test_throws DomainError !is_manifold_point(M,2*x,true)
            @test !is_tangent_vector(M,2*x,v)
            @test_throws DomainError !is_tangent_vector(M,2*x,v,true)
            @test !is_tangent_vector(M,x,y)
            @test_throws DomainError is_tangent_vector(M,x,y,true)
            test_manifold(
                M,
                pts,
                test_exp_log = false,
                default_inverse_retraction_method = PolarInverseRetraction(),
                test_injectivity_radius = false,
                test_is_tangent = true,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                projection_atol_multiplier = 15.0,
                retraction_atol_multiplier = 10.0,
                is_tangent_atol_multiplier = 4*10.0^2,
                retraction_methods = [PolarRetraction(), QRRetraction()],
                inverse_retraction_methods = [PolarInverseRetraction(), QRInverseRetraction()]
            )

            @testset "inner/norm" begin
                v1 = inverse_retract(M, pts[1], pts[2], PolarInverseRetraction())
                v2 = inverse_retract(M, pts[1], pts[3], PolarInverseRetraction())

                @test real(inner(M, pts[1], v1, v2)) ≈ real(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v2)) ≈ -imag(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v1)) ≈ 0

                @test norm(M, pts[1], v1) isa Real
                @test norm(M, pts[1], v1) ≈ sqrt(inner(M, pts[1], v1, v1))
            end
        end
    end

    @testset "Quaternion" begin
        M = Stiefel(3,2,ℍ)
        @testset "Basics" begin
            @test representation_size(M) == (3,2)
            @test manifold_dimension(M) == 18
        end
    end
end


text:
# SD_Error
export keysight_error

const SD_ERROR_OPENING_MODULE                   = Cint(-8000)
const SD_ERROR_CLOSING_MODULE                   = Cint(-8001)
const SD_ERROR_OPENING_HVI                      = Cint(-8002)
const SD_ERROR_CLOSING_HVI                      = Cint(-8003)
const SD_ERROR_MODULE_NOT_OPENED                = Cint(-8004)
const SD_ERROR_MODULE_NOT_OPENED_BY_USER        = Cint(-8005)
const SD_ERROR_MODULE_ALREADY_OPENED            = Cint(-8006)
const SD_ERROR_HVI_NOT_OPENED                   = Cint(-8007)
const SD_ERROR_INVALID_OBJECTID                 = Cint(-8008)
const SD_ERROR_INVALID_MODULEID                 = Cint(-8009)
const SD_ERROR_INVALID_MODULEUSERNAME           = Cint(-8010)
const SD_ERROR_INVALID_HVIID                    = Cint(-8011)
const SD_ERROR_INVALID_OBJECT                   = Cint(-8012)
const SD_ERROR_INVALID_NCHANNEL                 = Cint(-8013)
const SD_ERROR_BUS_DOES_NOT_EXIST               = Cint(-8014)
const SD_ERROR_BITMAP_ASSIGNED_DOES_NOT_EXIST   = Cint(-8015)
const SD_ERROR_BUS_INVALID_SIZE                 = Cint(-8016)
const SD_ERROR_BUS_INVALID_DATA                 = Cint(-8017)
const SD_ERROR_INVALID_VALUE                    = Cint(-8018)
const SD_ERROR_CREATING_WAVE                    = Cint(-8019)
const SD_ERROR_NOT_VALID_PARAMETERS             = Cint(-8020)
const SD_ERROR_AWG                              = Cint(-8021)
const SD_ERROR_DAQ_INVALID_FUNCTIONALITY        = Cint(-8022)
const SD_ERROR_DAQ_POOL_ALREADY_RUNNING         = Cint(-8023)
const SD_ERROR_UNKNOWN                          = Cint(-8024)
const SD_ERROR_INVALID_PARAMETERS               = Cint(-8025)
const SD_ERROR_MODULE_NOT_FOUND                 = Cint(-8026)
const SD_ERROR_DRIVER_RESOURCE_BUSY             = Cint(-8027)
const SD_ERROR_DRIVER_RESOURCE_NOT_READY        = Cint(-8028)
const SD_ERROR_DRIVER_ALLOCATE_BUFFER           = Cint(-8029)
const SD_ERROR_ALLOCATE_BUFFER                  = Cint(-8030)
const SD_ERROR_RESOURCE_NOT_READY               = Cint(-8031)
const SD_ERROR_HARDWARE                         = Cint(-8032)
const SD_ERROR_INVALID_OPERATION                = Cint(-8033)
const SD_ERROR_NO_COMPILED_CODE                 = Cint(-8034)
const SD_ERROR_FW_VERIFICATION                  = Cint(-8035)
const SD_ERROR_COMPATIBILITY                    = Cint(-8036)
const SD_ERROR_INVALID_TYPE                     = Cint(-8037)
const SD_ERROR_DEMO_MODULE                      = Cint(-8038)
const SD_ERROR_INVALID_BUFFER                   = Cint(-8039)
const SD_ERROR_INVALID_INDEX                    = Cint(-8040)
const SD_ERROR_INVALID_NHISTOGRAM               = Cint(-8041)
const SD_ERROR_INVALID_NBINS                    = Cint(-8042)
const SD_ERROR_INVALID_MASK                     = Cint(-8043)
const SD_ERROR_INVALID_WAVEFORM                 = Cint(-8044)
const SD_ERROR_INVALID_STROBE                   = Cint(-8045)
const SD_ERROR_INVALID_STROBE_VALUE             = Cint(-8046)
const SD_ERROR_INVALID_DEBOUNCING               = Cint(-8047)
const SD_ERROR_INVALID_PRESCALER                = Cint(-8048)
const SD_ERROR_INVALID_PORT                     = Cint(-8049)
const SD_ERROR_INVALID_DIRECTION                = Cint(-8050)
const SD_ERROR_INVALID_MODE                     = Cint(-8051)
const SD_ERROR_INVALID_FREQUENCY                = Cint(-8052)
const SD_ERROR_INVALID_IMPEDANCE                = Cint(-8053)
const SD_ERROR_INVALID_GAIN                     = Cint(-8054)
const SD_ERROR_INVALID_FULLSCALE                = Cint(-8055)
const SD_ERROR_INVALID_FILE                     = Cint(-8056)
const SD_ERROR_INVALID_SLOT                     = Cint(-8057)
const SD_ERROR_INVALID_NAME                     = Cint(-8058)
const SD_ERROR_INVALID_SERIAL                   = Cint(-8059)
const SD_ERROR_INVALID_START                    = Cint(-8060)
const SD_ERROR_INVALID_END                      = Cint(-8061)
const SD_ERROR_INVALID_CYCLES                   = Cint(-8062)
const SD_ERROR_HVI_INVALID_NUMBER_MODULES       = Cint(-8063)
const SD_ERROR_DAQ_P2P_ALREADY_RUNNING          = Cint(-8064)

error_dict = Dict(
    SD_ERROR_OPENING_MODULE =>
        "Keysight Error: Opening module",
    SD_ERROR_CLOSING_MODULE =>
        "Keysight Error: Closing module",
    SD_ERROR_OPENING_HVI =>
        "Keysight Error: Opening HVI",
    SD_ERROR_CLOSING_HVI =>
        "Keysight Error: Closing HVI",
    SD_ERROR_MODULE_NOT_OPENED =>
        "Keysight Error: Module not opened",
    SD_ERROR_MODULE_NOT_OPENED_BY_USER =>
        "Keysight Error: Module not opened by user",
    SD_ERROR_MODULE_ALREADY_OPENED =>
        "Keysight Error: Module already opened",
    SD_ERROR_HVI_NOT_OPENED =>
        "Keysight Error: HVI not opened",
    SD_ERROR_INVALID_OBJECTID =>
        "Keysight Error: Invalid ObjectID",
    SD_ERROR_INVALID_MODULEID =>
        "Keysight Error: Invalid ModuleID",
    SD_ERROR_INVALID_MODULEUSERNAME =>
        "Keysight Error: Invalid Module User Name",
    SD_ERROR_INVALID_HVIID =>
        "Keysight Error: Invalid HVI",
    SD_ERROR_INVALID_OBJECT =>
        "Keysight Error: Invalid Object",
    SD_ERROR_INVALID_NCHANNEL =>
        "Keysight Error: Invalid channel number",
    SD_ERROR_BUS_DOES_NOT_EXIST =>
        "Keysight Error: Bus doesn’t exist",
    SD_ERROR_BITMAP_ASSIGNED_DOES_NOT_EXIST =>
        "Keysight Error: Any input assigned to the bitMap does not exist",
    SD_ERROR_BUS_INVALID_SIZE =>
        "Keysight Error: Input size does not fit on this bus",
    SD_ERROR_BUS_INVALID_DATA =>
        "Keysight Error: Input data does not fit on this bus",
    SD_ERROR_INVALID_VALUE =>
        "Keysight Error: Invalid value",
    SD_ERROR_CREATING_WAVE =>
        "Keysight Error: Creating Waveform",
    SD_ERROR_NOT_VALID_PARAMETERS =>
        "Keysight Error: Invalid Parameters",
    SD_ERROR_AWG =>
        "Keysight Error: AWG function failed",
    SD_ERROR_DAQ_INVALID_FUNCTIONALITY =>
        "Keysight Error: Invalid DAQ functionality",
    SD_ERROR_DAQ_POOL_ALREADY_RUNNING =>
        "Keysight Error: DAQ buffer pool is already running",
    SD_ERROR_UNKNOWN =>
        "Keysight Error: Unknown error",
    SD_ERROR_INVALID_PARAMETERS =>
        "Keysight Error: Invalid parameter",
    SD_ERROR_MODULE_NOT_FOUND =>
        "Keysight Error: Module not found",
    SD_ERROR_DRIVER_RESOURCE_BUSY =>
        "Keysight Error: Driver resource busy",
    SD_ERROR_DRIVER_RESOURCE_NOT_READY =>
        "Keysight Error: Driver resource not ready",
    SD_ERROR_DRIVER_ALLOCATE_BUFFER =>
        "Keysight Error: Cannot allocate buffer in driver",
    SD_ERROR_ALLOCATE_BUFFER =>
        "Keysight Error: Cannot allocate buffer",
    SD_ERROR_RESOURCE_NOT_READY =>
        "Keysight Error: Resource not ready",
    SD_ERROR_HARDWARE =>
        "Keysight Error: Hardware error",
    SD_ERROR_INVALID_OPERATION =>
        "Keysight Error: Invalid Operation",
    SD_ERROR_NO_COMPILED_CODE =>
        "Keysight Error: No compiled code in the module",
    SD_ERROR_FW_VERIFICATION =>
        "Keysight Error: Firmware verification failed",
    SD_ERROR_COMPATIBILITY =>
        "Keysight Error: Compatibility error",
    SD_ERROR_INVALID_TYPE =>
        "Keysight Error: Invalid type",
    SD_ERROR_DEMO_MODULE =>
        "Keysight Error: Demo module",
    SD_ERROR_INVALID_BUFFER =>
        "Keysight Error: Invalid buffer",
    SD_ERROR_INVALID_INDEX =>
        "Keysight Error: Invalid index",
    SD_ERROR_INVALID_NHISTOGRAM =>
        "Keysight Error: Invalid histogram number",
    SD_ERROR_INVALID_NBINS =>
        "Keysight Error: Invalid number of bins",
    SD_ERROR_INVALID_MASK =>
        "Keysight Error: Invalid mask",
    SD_ERROR_INVALID_WAVEFORM =>
        "Keysight Error: Invalid waveform",
    SD_ERROR_INVALID_STROBE =>
        "Keysight Error: Invalid strobe",
    SD_ERROR_INVALID_STROBE_VALUE =>
        "Keysight Error: Invalid strobe value",
    SD_ERROR_INVALID_DEBOUNCING =>
        "Keysight Error: Invalid debouncing",
    SD_ERROR_INVALID_PRESCALER =>
        "Keysight Error: Invalid prescaler",
    SD_ERROR_INVALID_PORT =>
        "Keysight Error: Invalid port",
    SD_ERROR_INVALID_DIRECTION =>
        "Keysight Error: Invalid direction",
    SD_ERROR_INVALID_MODE =>
        "Keysight Error: Invalid mode",
    SD_ERROR_INVALID_FREQUENCY =>
        "Keysight Error: Invalid frequency",
    SD_ERROR_INVALID_IMPEDANCE =>
        "Keysight Error: Invalid impedance",
    SD_ERROR_INVALID_GAIN =>
        "Keysight Error: Invalid gain",
    SD_ERROR_INVALID_FULLSCALE =>
        "Keysight Error: Invalid fullscale",
    SD_ERROR_INVALID_FILE =>
        "Keysight Error: Invalid file",
    SD_ERROR_INVALID_SLOT =>
        "Keysight Error: Invalid slot",
    SD_ERROR_INVALID_NAME =>
        "Keysight Error: Invalid product name",
    SD_ERROR_INVALID_SERIAL =>
        "Keysight Error: Invalid serial number",
    SD_ERROR_INVALID_START =>
        "Keysight Error: Invalid start",
    SD_ERROR_INVALID_END =>
        "Keysight Error: Invalid end",
    SD_ERROR_INVALID_CYCLES =>
        "Keysight Error: Invalid number of cycles",
    SD_ERROR_HVI_INVALID_NUMBER_MODULES =>
        "Keysight Error: Invalid number of modules on HVI",
    SD_ERROR_DAQ_P2P_ALREADY_RUNNING =>
        "Keysight Error: DAQ P2P is already running")

keysight_error(errorCode) = error_dict[errorCode]


text:
using StatsBase

"""
    connectID = smoothen_connectID(connectID::Vector{Int32})

Make `connectID` consists of the set of integers `1:length(unique(connectID))``
"""
function smoothen_connectID(connectID::Vector{Int64})
    # Force `connectID` to consist of the integers
    # `1:length(unique(connectID))` without skipping any integers.
    sortindices = sortperm(connectID)
    connectID_sorted = connectID[sortindices]
    connectIDlist = unique(connectID_sorted)
    nperID = counts(connectID_sorted)
    nperID = nperID[nperID .!= 0]
    ncumulative = [0; cumsum(nperID)]
    for nn = 1:length(connectIDlist)
        clusterind = (1:nperID[nn]) .+ ncumulative[nn]
        connectID_sorted[clusterind] .= nn
    end
    connectID[sortindices] = connectID_sorted

    return connectID
end


text:
@testset "ConstantRate" begin
    for rng in [nothing, Random.MersenneTwister()]
        T = 10_000
        freq = rand(1:100)
        dt = rand() * (1/freq)
        rate = freq * dt
        input = isnothing(rng) ? ConstantRate(rate) : ConstantRate(rate; rng = rng)
        # Averaging the spike count of ConstantRate over a fixed window of time should be ConstantRate.rate
        @test isapprox(count(x -> x > 0, [input(t) for t in 1:T]) / T, rate, atol=0.005)
        # Assert frequency constructor matches the rate constructor
        @test ConstantRate(freq, dt).rate == input.rate
    end
end

# A StepCurrent should be zero < threshold time and positive > threshold time
@testset "StepCurrent" begin
    τ = 10
    T = 25
    stepcurrent = StepCurrent(τ)
    prethresh = [stepcurrent(t) for t in 1:τ]
    postthresh = [stepcurrent(t) for t in (τ + 1):T]

    @test all(x -> x == 0, prethresh)
    @test all(x -> x > 0, postthresh)
end

# The spikes of a inhomogeneous Poisson input with a constant λ should be distributed as a Poisson process
@test_skip @testset "PoissonInput" begin
    for rng in [nothing, Random.MersenneTwister()]
        ρ₀ = 0.1
        λ = 0.2
        pI = isnothing(rng) ? PoissonInput(ρ₀, (t; dt) -> λ) : pI = PoissonInput(ρ₀, (t; dt) -> λ; rng = rng)
        lasttime = 0
        outputs = Int[]
        for t in 1:10_000
            if (pI(t) > 0)
                push!(outputs, t - lasttime) # Store the time difference between the last time a spike was observed and now
                lasttime = t # Update the time when last spike was observed
            end
        end
        ds = Exponential(λ)
        ExactOneSampleKSTest(outputs, ds) # Test if sample outputs is from same distribution
    end
end


text:
#md # ```@meta
#md # CurrentModule = DitherPunk
#md # ```
#md #
#md # # DitherPunk.jl
#md # A [dithering / digital halftoning](https://en.wikipedia.org/wiki/Dither) package inspired by [Lucas Pope's Obra Dinn](https://obradinn.com) and [Surma's blogpost](https://surma.dev/things/ditherpunk/) of the same name.
#md #
#md # !!! note
#md #     This package is part of a wider [Julia-based image processing ecosystem](https://github.com/JuliaImages). If you are starting out, then you may benefit from reading about some [fundamental conventions](https://juliaimages.org/latest/quickstart/) that the ecosystem utilizes that are markedly different from how images are typically represented in OpenCV, MATLAB, ImageJ or Python.
#md #
# # Getting started
# We start out by loading an image, in this case the lighthouse
# from [*TestImages.jl*](https://testimages.juliaimages.org).
using DitherPunk
using Images
using TestImages

img = testimage("lighthouse")
img = imresize(img; ratio=1//2)

# To apply binary dithering, we also need to convert the image to grayscale.
img_gray = Gray.(img)

#
# !!! note " Preprocessing"
#     Sharpening the image and adjusting the contrast can emphasize the effect of the algorithms. It is highly recommended to play around with algorithms such as those provided by [ImageContrastAdjustment.jl](https://juliaimages.org/ImageContrastAdjustment.jl/stable/)
#
# ## Binary dithering
# Since we already turned the image to grayscale, we are ready to apply Bayer dithering,
# an [ordered dithering](https://en.wikipedia.org/wiki/Ordered_dithering) algorithm that leads to characteristic cross-hatch patterns.
dither(img_gray, Bayer())

# ### Color spaces
# Depending on the method, dithering in sRGB color space can lead to results that are too bright.
# To obtain a dithered image that more closely matches the human perception of brightness, grayscale images can be converted to linear color space using the boolean keyword argument `to_linear`.
dither(img_gray, Bayer(); to_linear=true)

# ## Separate-space dithering
# All dithering algorithms in DitherPunk can also be applied to color images
# and will automatically apply channel-wise binary dithering.
dither(img, Bayer())

#
# !!! note
#     Because the algorithm is applied once per channel, the output of this algorithm depends on the color type of input image. `RGB` is recommended, but feel free to experiment!
#
# ## Dithering with custom colors
# Let's assume we want to recreate an image by stacking a bunch of Rubik's cubes. Dithering algorithms are perfect for this task!
# We start out by defining a custom color scheme:
white = RGB{Float32}(1, 1, 1)
yellow = RGB{Float32}(1, 1, 0)
green = RGB{Float32}(0, 0.5, 0)
orange = RGB{Float32}(1, 0.5, 0)
red = RGB{Float32}(1, 0, 0)
blue = RGB{Float32}(0, 0, 1)

rubiks_colors = [white, yellow, green, orange, red, blue]

# All methods of type `ErrorDiffusion` and `OrderedDither` support dithering in custom colors.
d = dither(img, FloydSteinberg(), rubiks_colors)

# this looks much better than simply quantizing to the closest color!
d = dither(img, ClosestColor(), rubiks_colors)

# ### ColorSchemes.jl
# Predefined color schemes from [ColorSchemes.jl](https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/#Pre-defined-schemes) can also be used.
using ColorSchemes

dither(img, FloydSteinberg(), ColorSchemes.PuOr_7)

# You can also directly use the corresponding symbol from the
# [ColorSchemes catalogue](https://juliagraphics.github.io/ColorSchemes.jl/stable/catalogue/):
dither(img, FloydSteinberg(), :PuOr_7)

# ### Clustering.jl
# Using [Clustering.jl](https://github.com/JuliaStats/Clustering.jl) allows you to generate
# optimized color schemes. Simply pass the size of the desired color palette:
using Clustering

dither(img, FloydSteinberg(), 8)

# ## UnicodePlots.jl
# Using [UnicodePlots.jl](https://github.com/JuliaPlots/UnicodePlots.jl), it is also possible
# to dither images directly to Braille-characters using `braille`. The interface is the same
# as for binary dithering with `dither`:
using UnicodePlots
img = imresize(img; ratio=1//3)

braille(img, FloydSteinberg())


text:
using Test

# H[1:33:end,30:33:end]
tests_2d = Dict("diffusion_2D_expl.jl"=> [1.4341623497148877e-14 2.021424041084888e-11 4.8022139161215363e-14;
                                          2.394095924752031e-6 0.003374444698589263 8.016498807587078e-6;
                                          0.0006736119165607539 0.594473121056432 0.002255555227680872;
                                          3.1945000985042113e-7 0.00045025860455410284 1.069660825438506e-6],
                "diffusion_2D_impl.jl"=> [1.4341623497148877e-14 2.021424041084888e-11 4.8022139161215363e-14;
                                          2.3940959247520396e-6 0.003374444691354641 8.016498807586814e-6;
                                          0.0006736119165451054 0.6115608556980482 0.0022555552261044605;
                                          3.1945000985042113e-7 0.0004502586045507441 1.069660825438506e-6],
                "diffusion_2D_damp.jl"=> [1.4341623497148877e-14 2.021424041084888e-11 4.8022139161215363e-14;
                                          2.394095924752039e-6 0.0033744446790461404 8.016498807586816e-6;
                                          0.0006736119165187302 0.6106677330399702 0.0022555552234409136;
                                          3.1945000985042113e-7 0.00045025860454511377 1.0696608254385066e-6]
                )

## Test the scripts
# note: this cannot be run in a test-set as the include acts in global scope

do_visu, save_fig = false, false # disable plotting & saving

for (fl, H_ref) in tests_2d
    lx, ly = 10.0, 10.0         # domain size
    nx, ny = 128, 128           # numerical grid resolution
    dx, dy = lx/nx, ly/ny       # grid size

    println("Runninng $fl")
    xco, yco, Hco = include(joinpath(@__DIR__, "../scripts", fl))
    @test xco == LinRange(dx/2, lx-dx/2, nx)
    @test yco == LinRange(dy/2, ly-dy/2, ny)
    @test Hco[1:33:end,30:33:end] ≈ H_ref
end


text:
# Julia 0.6.1 code
BLAS.set_num_threads(8)

function perceptronCycle!(w,g,X,y,eps,sigma_in,sigma_out,rho)
    P=length(X)
    err=0
    for ind = randperm(P)
		eta = sqrt(w'*w*sigma_in^2 + sigma_out^2)*randn()
        h = w'*X[ind]-1 + eta
        if y[ind]*h[1]<0.
            err+=1
            w[:] = (1.0-eps) * w + rho*y[ind]*X[ind]
            w[:] = ((g.*w).>0.).*w
        else
			w[:] = (1.0-eps) * w
		end
    end
    return err
end

function sim(N,P,sigma_w_init,sigma_in,sigma_out,
             rho,eps,f_exc,f_in,p_out,T,stat_interval)
    # g is vector of types of neurons
    g=ones(N)
    g[Int(N*f_exc):end]=-1
    # input patterns
    X=[map(Float64,rand(N).>f_in) for μ=1:P]
    # desired output
    y=2*map(Float64,rand(P).>p_out)-1
    # initial weights
    w=sigma_w_init*randexp(N).*g

    TT = Int(T/stat_interval)
    iter=Array{Float64}(TT)
    err=Array{Float64}(TT)
    I=Array{Float64}(TT)
    IB=Array{Float64}(TT)
    w_norm=Array{Float64}(TT)
    kappa_in=Array{Float64}(TT)
    kappa_out=Array{Float64}(TT)

    for t=1:T
        cycle_error = perceptronCycle!(w,g,X,y,eps,
						   			 sigma_out,sigma_in,rho)/P
        if t % stat_interval == 0
            #Gathering stats
            h=Array{Float64}([(w'*x)[1] for x=X])
            hbar=mean(h)
            gw=g.*w
            C=Array{Float64}([(gw'*x)[1] for x=X])
            kappa_mu=y.*(h-1)
            kout=minimum(kappa_mu)
            wn=norm(w)
            ind = Int(t/stat_interval)
            # saving stats
            err[ind] = cycle_error
            kappa_out[ind] = kout
            kappa_in[ind] = kout/wn
            I[ind] = hbar/(sqrt(N)*std(h))
            IB[ind] = hbar/mean(C)
            w_norm[ind] = wn
            iter[ind] = t
        end
    end

    return iter,err,I,IB,w_norm,kappa_in,kappa_out
end

##
P=900
N=3000

sigma_w_init = 1e-3
sigma_out = 0.15
sigma_in = 0.1
rho = 0.1/N
eps = 5e-7

f_in = 0.5
p_out = 0.5

f_exc = 0.8

T = 10000
stat_interval = 10
TT  = Int(T/stat_interval)
n_sim = 2

err1=Array{Float64}(n_sim,TT)
I1=Array{Float64}(n_sim,TT)
w_norm1=Array{Float64}(n_sim,TT)
IB1=Array{Float64}(n_sim,TT)
kappa_in1=Array{Float64}(n_sim,TT)
kappa_out1=Array{Float64}(n_sim,TT)

for n=1:n_sim
    print(n,"\n")
    t,e,i,ib,wn,kin,kout=sim(N,P,sigma_w_init,
                             sigma_in,sigma_out,
							 rho,eps,
                             f_exc,f_in,p_out,
                             T,stat_interval)
    err1[n,:]=e
    I1[n,:]=i
    IB1[n,:]=ib
    w_norm1[n,:]=wn
    kappa_in1[n,:]=kin
    kappa_out1[n,:]=kout
end
##
using Plots
pyplot()
##
p1 = plot(mean(err1,1)',xscale=:log10,label="Error")
p2 = plot(mean(IB1,1)',xscale=:log10,ylim=(0,0.4),label="IB")
p3 = plot(mean(kappa_in1,1)',xscale=:log10,ylim=(0,0.4),label="κ_in")
p4 = plot(mean(kappa_out1,1)',xscale=:log10,ylim=(0,0.4),label="κ_out")
plot(p1,p2,p3,p4,layout=(2,2))


text:
using BandedMatrices

function example(n)
    t = [ones(n)'*4; ones(n)']
    s = [2ones(n)'; ones(n)']
    tt = Symmetric(BandedMatrix(0=>t[1,:], -1=>t[2,1:end-1]), :L)
    ss = Symmetric(BandedMatrix(0=>s[1,:], -1=>s[2,1:end-1]), :L)
    t, s, tt, ss
end

@testset "bounds($n)" for n = [1,2,3,4,5,7,10,199]
    t, s, = example(n)
    lb00, ub00 = bounds(t, s)
    #lb01, ub01 = bounds(t, -s)
    lb10, ub10 = bounds(-t, s)
    #lb11, ub11 = bounds(-t, -s)
    @test lb00 <= ub00
    #@test lb00 == lb11 && ub00 == ub11
    #@test lb01 == lb10 && ub01 == ub10
    @test lb10 == -ub00 && ub10 == -lb00
end


text:
mutable struct MemoryChunk
    mMemory::vk.VkDeviceMemory
    mVkDevice::vk.VkDevice
    mAllocationOffset::vk.VkDeviceSize
    mOffset::vk.VkDeviceSize
    mSize::vk.VkDeviceSize
    mType::UInt32
    mMappable::Bool
    mDeallocate

    function MemoryChunk(memory::vk.VkDeviceMemory, vkDevice::vk.VkDevice,
                         allocationOffset::vk.VkDeviceSize, offset::vk.VkDeviceSize,
                         size::vk.VkDeviceSize, type::UInt32, mappable::Bool,
                         deallocator = nothing)
        this = new()
        this.mMemory = memory
        this.mVkDevice = vkDevice
        this.mAllocationOffset = allocationOffset
        this.mOffset = offset
        this.mSize = size
        this.mType = type
        this.mMappable = mappable
        this.mDeallocate = deallocator
        return this
    end
end

@class MemoryChunk

function destroy(this::MemoryChunk)
    if (this.mSize == 0)
        return
    end
    if (this.mDeallocate != nothing)
        this.mDeallocate(this)
    else
        vk.vkFreeMemory(this.mVkDevice, this.mMemory, C_NULL)
    end
end

mutable struct MappedMemory
    mOffset::vk.VkDeviceSize
    mSize::vk.VkDeviceSize
    mVkDevice::vk.VkDevice
    mMemory::vk.VkDeviceMemory
    mData::Ptr{Cvoid}

    function MappedMemory(mem::MemoryChunk, dataOffset::vk.VkDeviceSize, size::vk.VkDeviceSize)
        this = new()
        this.mMemory = mem.mMemory
        this.mVkDevice = mem.mVkDevice
        this.mOffset = getOffset(mem) + dataOffset
        this.mSize = min(mem.mSize - dataOffset, size)
        this.mData = VkExt.mapMemory(this.mVkDevice, this.mMemory, this.mOffset, this.mSize, vk.VkFlags(0))
        return this
    end
end
@class MappedMemory

function handle(this::MemoryChunk)::vk.VkDeviceMemory
    return this.mMemory
end

function bindToImage(this::MemoryChunk, image::vk.VkImage)
    if vk.vkBindImageMemory(this.mVkDevice, image, this.mMemory, this.mOffset) != vk.VK_SUCCESS
        error("Failed to bind memory chunk to image!")
    end
end

function bindToBuffer(this::MemoryChunk, buffer::vk.VkBuffer)
    if vk.vkBindBufferMemory(this.mVkDevice, buffer, this.mMemory, this.mOffset) != vk.VK_SUCCESS
        error("Failed to bind memory chunk to buffer!")
    end
end

function getOffset(this::MemoryChunk)::vk.VkDeviceSize
    return this.mOffset
end

function getSize(this::MemoryChunk)::vk.VkDeviceSize
    return this.mSize
end

function getAllocationOffset(this::MemoryChunk)::vk.VkDeviceSize
    return this.mAllocationOffset
end

function isMappable(this::MemoryChunk)::Bool
    return this.mMappable
end

function map(this::MemoryChunk, dataOffset::vk.VkDeviceSize = vk.VkDeviceSize(0), size::vk.VkDeviceSize = vk.VkDeviceSize(vk.VK_WHOLE_SIZE))::MappedMemory
    return MappedMemory(this, dataOffset, size)
end

# When ~MappedMemory is called
function unmap(this::MappedMemory)
    if this.mMemory != C_NULL
        range = [vk.VkMappedMemoryRange(
                    vk.VK_STRUCTURE_TYPE_MAPPED_MEMORY_RANGE, #sType::VkStructureType
                    C_NULL, #pNext::Ptr{Cvoid}
                    this.mMemory, #memory::VkDeviceMemory
                    this.mOffset, #offset::VkDeviceSize
                    vk.VK_WHOLE_SIZE #size::VkDeviceSize
                )]
        vk.vkFlushMappedMemoryRanges(this.mVkDevice, 1, pointer(range))
        vk.vkUnmapMemory(this.mVkDevice, this.mMemory)
    end
end

function getData(this::MappedMemory)::Ptr{Cvoid}
    return this.mData
end

function getSize(this::MappedMemory)::vk.VkDeviceSize
    return this.mSize
end


text:
"""
    read_graph(fname::String, format::String, delim::AbstractChar='')

Reads in a delimited file with 2 (simple) or 3 columns into a `Graph` struct
"""
function read_graph(fname::String, format::String; delim::AbstractChar='\t')::Graph
    data, header = readdlm(fname, delim, Int64, header=true)
    n_nodes, n_edges, = map(x -> parse(Int64, x), header)
    nodes = map(x -> Node(x, 0, 0), 1:n_nodes)
    sorted_data = sortslices(data, dims = 1)
    first_out, end_node = forward_star(n_nodes, sorted_data[:, 1:2])
    if format == "simple"
        edges = map(x -> Edge(x[1], 0, 0), 1:n_edges)
    elseif format == "distance"
        edges = map(x -> Edge(x[1], x[2], 0), zip(1:n_edges, sorted_data[:, 3]))
    elseif format == "capacity"
        edges = map(x -> Edge(x[1], 0, x[2]), zip(1:n_edges, sorted_data[:, 3]))
    elseif format == "assignment"
        edges = map(x -> Edge(x[1], x[2], 0), zip(1:n_edges, sorted_data[:, 3]))
    else
        throw(ArgumentError("Unsupported format. Use either distance, capacity or assignment"))
    end
    return Graph(nodes, edges, first_out, end_node)
end


text:
using Test, SpecialFunctions, Statistics
import TruncatedNormal as TN

@test isnan(TN.tnvar(1, 0))
@test isnan(TN.tnvar(+Inf, 0))
@test isnan(TN.tnvar(0, -Inf))
@test isnan(TN.tnvar(NaN, NaN))
@test isnan(TN.tnvar(0, NaN))
@test isnan(TN.tnvar(NaN, 0))

@test TN.tnvar(-Inf, +Inf) == 1
@test TN.tnvar(+Inf, +Inf) == 0
@test TN.tnvar(-Inf, -Inf) == 0
@test TN.tnvar(0, +Inf) ≈ 1 - 2/π
@test TN.tnvar(-Inf, 0) ≈ 1 - 2/π

for x = -10:10
    @test TN.tnvar(x, +Inf) ≈ 1 + √(2/π) * x / erfcx(+x / √2) - 2/π / erfcx(+x / √2)^2
    @test TN.tnvar(-Inf, x) ≈ 1 - √(2/π) * x / erfcx(-x / √2) - 2/π / erfcx(-x / √2)^2
end

@test TN.tnvar(50, 70) ≈ 0.0003990431868038995479099272265360593305365

@test TN.tnvar(-Inf, Inf, 0, 1) == 1

for x = -10:10
    @test TN.tnvar(x, +Inf, 0, 1) ≈ TN.tnvar(x, +Inf)
    @test TN.tnvar(-Inf, x, 0, 1) ≈ TN.tnvar(-Inf, x)
end

@test TN.tnvar(50, 70, -2, 3) ≈ 0.029373438107168350377591231295634273607812172191712
@test TN.tnvar(-100.0, 0.0, 0.0, 2.0986317998643735) ≈ 1.6004193412141677189841357987638847137391508803335
@test TN.tnvar(0.0, 0.9, 0.0, 0.07132755843183151) ≈ 0.0018487407287725028827020557707636415445504260892486
@test TN.tnvar(-100.0, 100.0, 0.0, 17.185261847875548) ≈ 295.333163899557735486302841237124507431445
@test TN.tnvar(-100.0, 0.5, 0.0, 0.47383322897860064) ≈ 0.145041095812679283837328561547251019229612
@test TN.tnvar(-100.0, 100.0, 0.0, 17.185261847875548) ≈ 295.333163899557735486302841237124507431445
@test TN.tnvar(-10000, 10000, 0, 1) ≈ 1

# https://github.com/JuliaStats/Distributions.jl/issues/827
@test TN.tnvar(999000, 1e6) ≥ 0
@test TN.tnvar(0, 1000, 1000000, 1) ≥ 0

@test_broken TN.tnvar(1e6, Inf) ≈ 9.99999999994000000000050000000e-13
@test_broken TN.tnvar(999000, 1e6) ≈ 1.00200300399898194688784897455e-12
@test_broken TN.tnvar(0, 1000, 1000000, 1) ≈ 1.00200300399898194688784897455909681215272616880064850605775930845141e-12


text:
# sinc interpolation operator

## helper module
module joSincInterp_etc
    using SpecialFunctions
    function kaiser_window(x,r,b)
        return abs(x) <= r ? besseli(0,b*sqrt(1-(x/r)^2))/besseli(0,b) : 0.0
    end
end
using .joSincInterp_etc

export joSincInterp
"""
    julia> joSincInterp(xin,xout;[r=...,][DDT=...,][RDT=...,][name=...])

sinc interpolation matrix for interpolating functions f defined on grid xin to functions defined on grid xout

# Signature

    function joSincInterp(xin::AbstractArray{T,1},xout::AbstractArray{T,1};
        r::I=0,DDT::DataType=T,RDT::DataType=DDT,name::String="joSincInterp")
        where {T<:AbstractFloat,I<:Integer}

# Arguments

- `xin`  - 1D input grid
- `xout` - 1D output grid
- keywords
    - `r`: kaiser window parameter (default: 0, no windowing)
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Notes

- If xout has more than one point, the spacings of xin and xout are normalized to the spacing of xout.

# Examples

% description

    joSincInterp(xin,xout)

examples with DDT/RDT

    % joSincInterp(xin,xout; DDT=Float32)
    % joSincInterp(xin,xout; DDT=Float32,RDT=Float64)

"""
function joSincInterp(xin::AbstractArray{T,1},xout::AbstractArray{T,1};
    r::I=0,DDT::DataType=T,RDT::DataType=DDT,name::String="joSincInterp") where {T<:AbstractFloat,I<:Integer}

    length(xin)>1 ? dx = min(xin[2]-xin[1], xout[2]-xout[1]) : dx = 1
    S = sinc.( (xout .- xin') ./ dx )
    if r > 0
        r_b =[1.24 2.94 4.53 6.31 7.91 9.42 10.95 12.53 14.09 14.18];
        window = [joSincInterp_etc.kaiser_window(xout[i]-xin[j],r,r_b[r]) for i in 1:length(xout), j in 1:length(xin)]
        S = S .* window
    end
    return joMatrix(S,DDT=DDT,RDT=RDT,name=name)

end