text:
module R_constraints
    include(dirname(@__FILE__)*"\\Utils_cons.jl")
    include(dirname(dirname(@__FILE__))*"\\Moments.jl")
    using .Utils_cons ; using .Moments
    const uc = Utils_cons ; const mom = Moments

    export make_mon_expo_keys,
           make_PSD_con,
           make_ord4_con,
           make_loc_cons_S1,
           make_loc_cons_S2,
           make_loc_cons_S3,
           make_G_con

    make_mon_expo_keys(n,t::Int) = Moments.make_mon_expo(n,t*2)

    """L([x,y]ₜ[x,y]ₜᵀ) ⪰ 0"""
    function make_PSD_con(d,t,Lx)
        MB = mom.get_ℝ_block_diag(d,t)
        Bs = [keys(MB)...]
        PSD_con = Dict(zip(Bs, map(b -> uc.idx2var(Lx, MB[b]),Bs)))
        return PSD_con
    end

    """L(xxᵀ⊗yyᵀ)"""
    make_ord4_con(d,Lx) = uc.idx2var(Lx,mom.make_xxᵀ⨂yyᵀ(d))

    """ L(g⋅η) ⪰ 0 for η ∈ even-degree-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ
                        g = √(ρₘₐₓ) - xᵢ² , √(ρₘₐₓ) - yᵢ² and i ∈ [d] """
    function make_loc_con(eᵢ,sqr_ρ,B,Lx)
        Lxᵢ = uc.idx2var(Lx, B + repeat([2*eᵢ], size(B)...)) #  L(xᵢ²⋅η)
        Lρ = sqr_ρ*uc.idx2var(Lx, B)  # L(√ρₘₐₓ⋅η)
        return  Lρ - Lxᵢ
    end

    function make_loc_cons_S1(ρ,d,t,Lx)
        n     = sum(d)
        MB    = mom.get_ℝ_block_diag(d,t.- 1)
        sqr_ρ = sqrt(maximum(ρ))

        loc_con = Dict()
        for b in keys(MB)
            for k in 1:d[1]
                loc_con[(b,"x²_$k")] = make_loc_con(uc.eᵢ(n,k),sqr_ρ,MB[b],Lx) # Constraint:  L((√(ρₘₐₓ) - xᵢ²))⋅η) ⪰ 0
            end
            for k in 1:d[2]
                loc_con[(b,"y²_$k")] = make_loc_con(uc.eᵢ(n,k+d[1]),sqr_ρ,MB[b],Lx) # Constraint:  L((√(ρₘₐₓ) - yᵢ²))⋅η) ⪰ 0
            end
        end
        return loc_con
    end

    """loc_con: L((√(Tr(ρ)) - ∑xᵢ²)⋅η) ⪰ 0 ∀ η ∈ even-degree-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ
            g₂: L((∑xᵢ² -  ∑yᵢ²)⋅η) = 0 ∀ η ∈ even-degree-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ """
    function make_loc_cons_S2(ρ,d,t,Lx)
        d₁,d₂ = d ; n = sum(d) ; stρ = sqrt(sum([ρ[i,i] for i in 1:n]))
        MB = mom.get_ℝ_block_diag(d,t.- 1)
        loc_con = Dict() ; g₂ = Dict()
        for b in keys(MB)
            xRterm = sum([uc.idx2var(Lx,MB[b] .+ 2* [uc.eᵢ(n,k)])    for k in 1:d₁]) # L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            yRterm = sum([uc.idx2var(Lx,MB[b] .+ 2* [uc.eᵢ(n,k+d₁)]) for k in 1:d₂]) # L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            loc_con[b,"x"] = stρ*uc.idx2var(Lx,MB[b])-xRterm # √Tr(ρ)⋅L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            # loc_con[b,"y"] = stρ*uc.idx2var(Lx,MB[b])-yRterm # √Tr(ρ)⋅L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            g₂[b] = xRterm - yRterm
        end
        return loc_con, g₂
    end

    """loc_con: L((Tr(ρ) - ∑xᵢ²)⋅η) ⪰ 0 ∀ η ∈ even-degree-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ
     ideal_con: L((1 -  ∑yᵢ²)⋅η) = 0 ∀ η ∈ even-degree-principle-submatrices of [x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ """
    function make_loc_cons_S3(ρ,d,t,Lx)
        d₁,d₂ = d ; n = sum(d) ; trρ = sum([ρ[i,i] for i in 1:size(ρ)[1]])
        MB    = mom.get_ℝ_block_diag(d,t.- 1)
        loc_con = Dict() ; ideal_con = Dict()
        for b in keys(MB)
            MBLx   = uc.idx2var(Lx, MB[b]) # L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ)
            xRterm = sum([uc.idx2var(Lx, MB[b] .+ [2*uc.eᵢ(n,k)])    for k in 1:d₁]) # L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            yRterm = sum([uc.idx2var(Lx, MB[b] .+ [2*uc.eᵢ(n,k+d₁)]) for k in 1:d₂]) # L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ )
            loc_con[b] = trρ*MBLx-xRterm #  tr(ρ)L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑xᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ ) ⪰ 0
            ideal_con[b] = MBLx - yRterm #  L([x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ) - L((∑yᵢ²)⋅[x,y]ₜ₋₁[x,y]ₜ₋₁ᵀ ) = 0
        end
        return loc_con, ideal_con
    end

    """ρ⊗L(η) - L((xxᵀ⊗yyᵀ) ⊗ η ) ⪰ 0 ∀ η ∈ even-degree-principle-submatrices of ([x,y]₌ₜ₋₂[x,y]₌ₜ₋₂ᵀ) """
    function make_G_con(ρ,d,t,Lx)
        MB = mom.get_ℝ_block_diag(d,t.- 2) ; xxᵀ⨂yyᵀ = mom.make_xxᵀ⨂yyᵀ(d)

        G_con = Dict()
        for b in keys(MB)
            Lterm   = kron(ρ, uc.idx2var(Lx,MB[b]))  # ρ⊗L([x, y]ₜ₋₂[x, y]ᵀₜ₋₂)
            Rterm   = uc.idx2var(Lx, mom.var_kron(xxᵀ⨂yyᵀ,MB[b]))  # L( (xxᵀ⊗yyᵀ) ⊗ ([x, y]ₜ₋₂[x, y]ᵀₜ₋₂) )
            G_con[b] = Lterm - Rterm # ρ⊗L([x,y]ₜ₋₂[x,y]ᵀₜ₋₂) - L( (xxᵀ⊗yyᵀ) ⊗ ([x,y]ₜ₋₂[x,y]ᵀₜ₋₂) )
        end
        return G_con
    end
end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
using BQCESubroutine
using BQCESubroutine.Utilities
using LinearAlgebra
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using Test
using StaticArrays
using BenchmarkTools
using YaoArrayRegister
using LoopVectorization
using BQCESubroutine: subspace_mul!, subspace_mul_generic!, subspace_mul4x4!
using PaddedMatrices

function subspace_smm!(st::AbstractMatrix, indices, U::AbstractMatrix, subspace)
    @inbounds @views for k in subspace
        idx = k .+ indices
        for b in 1:size(st, 1)
            st[b, idx] = U * st[b, idx]
        end
    end
    return st
end

n = 20
locs = Locations((3,4,5))
subspace = bsubspace(n, locs)
comspace = bcomspace(n, locs)
U = rand(ComplexF64, 8, 8)
static_indices = @SVector [i+1 for i in comspace]
stride_indices = StrideArray(MVector(static_indices));
st=rand(ComplexF64, 100, 1 << n)
subspace_smm!(copy(st), static_indices, U, subspace) ≈ subspace_mul_generic!(copy(st), stride_indices, U, subspace)

subspace_smm!(copy(st), static_indices, U, subspace) ≈ subspace_mul!(copy(st), comspace, U, subspace)
@benchmark subspace_smm!(st, $static_indices, $U, $subspace) setup=(st=rand(ComplexF64, 100, 1 << $n))
@benchmark subspace_mul_generic!(st, $stride_indices, $U, $subspace) setup=(st=rand(ComplexF64, 100, 1 << $n))
@benchmark subspace_mul!(st, $comspace, $U, $subspace) setup=(st=rand(ComplexF64, 100, 1 << $n))

locs = Locations((3,5,15))
broutine!(copy(st), U, locs) ≈ transpose(instruct!(transpose(copy(st)), U, Tuple(locs)))
@benchmark broutine!($(copy(st)), $U, $locs)
@benchmark instruct!(st, $U, $(Tuple(locs))) setup=(st=transpose(rand(ComplexF64, 100, 1 << $n)))

n = 20
locs = Locations((3,4,5,6))
subspace = bsubspace(n, locs)
comspace = bcomspace(n, locs)
U = rand(ComplexF64, 16, 16)
static_indices = @SVector [i+1 for i in comspace]
stride_indices = StrideArray(MVector(static_indices));
st=rand(ComplexF64, 1 << n)
broutine!(copy(st), U, locs) ≈ instruct!(copy(st), U, Tuple(locs))
@benchmark broutine!(st, $U, $locs) setup=(st=rand(ComplexF64, 1 << $n))
@benchmark instruct!(st, $U, $(Tuple(locs))) setup=(st=rand(ComplexF64, 1 << $n))

ext: jl


text:
export ROCRANDError

struct ROCRANDError <: Exception
    code::rocrand_status
    msg::AbstractString
end
Base.show(io::IO, err::ROCRANDError) = print(io, "ROCRANDError(code $(err.code), $(err.msg))")

function ROCRANDError(code::rocrand_status)
    msg = status_message(code)
    return ROCRANDError(code, msg)
end

function status_message(status)
    if status == ROCRAND_STATUS_SUCCESS
        return "the operation completed successfully"
    elseif status == ROCRAND_STATUS_VERSION_MISMATCH
        return "header file and linked library version do not match"
    elseif status == ROCRAND_STATUS_NOT_CREATED
        return "generator was not created"
    elseif status == ROCRAND_STATUS_ALLOCATION_FAILED
        return "memory allocation failed"
    elseif status == ROCRAND_STATUS_TYPE_ERROR
        return "wrong generator type"
    elseif status == ROCRAND_STATUS_OUT_OF_RANGE
        return "argument out of range"
    elseif status == ROCRAND_STATUS_LENGTH_NOT_MULTIPLE
        return "requested size is not a multiple of quasirandom generator's dimension"
    elseif status == ROCRAND_STATUS_DOUBLE_PRECISION_REQUIRED
        return "GPU does not support double precision"
    elseif status == ROCRAND_STATUS_LAUNCH_FAILURE
        return "kernel launch failure"
    elseif status == ROCRAND_STATUS_INTERNAL_ERROR
        return "internal library error"
    else
        return "unknown error"
    end
end

macro check(rand_func)
    quote
        local err::rocrand_status
        err = $(esc(rand_func::Expr))
        if err != ROCRAND_STATUS_SUCCESS
            throw(ROCRANDError(err))
        end
        err
    end
end


text:
using AppleAccelerateLinAlgWrapper
using Test

@testset "AppleAccelerateLinAlgWrapper.jl" begin

    for T ∈ (Float32,Float64)
        A = rand(T,10,10);
        B = rand(T,10,10);
        @test AppleAccelerateLinAlgWrapper.gemm(A,B) ≈ A*B
        @test AppleAccelerateLinAlgWrapper.rdiv(A,B) ≈ A/B
        @test AppleAccelerateLinAlgWrapper.ldiv(A,B) ≈ A\B
    end
end


text:
using BivMatFun;
using BenchmarkTools;
using Printf;
using LinearAlgebra;
using DelimitedFiles;
using MAT;

include("experiments_common.jl")

function run_test()

    success = true;

    matrices = [ "jordbloc1", "grcar", "smoke", "kahan2", "lesp", "sampling", "grcar-randn" ];

    m = 64;
    n = m;
    a = .5;

    fncts = [
    Dict("name" => "sqrt", "f" => (x,y,i,j) -> sqrt(x+y)),
    Dict("name" => "invsqrt", "f" => (x,y,i,j)  -> 1 ./ sqrt(x+y)),
    Dict("name" => "phi1", "f" => (x,y,i,j) -> expm1(x+y) ./ (x+y)),
    Dict("name" => "expa", "f" => (x,y,i,j) -> exp((x+y).^a))
    ]

    shift = 1.00;

    lesp = n -> convert(Matrix{ComplexF64}, diagm(-1 => 1 ./ (2:n), 0 => -(5 : 2 : 2*n+3), 1 => 2 : n));
    smoke = n -> diagm(0 => exp.(2im * pi * (1 : n) ./ n), 1 => complex(ones(n-1)), -n+1 => [ 1.0 + 0im ]);

    # err_fun2m, err_diag, err_diagm, time_fun2m, time_diag, time_diagm, nblocksA, nblocksB, d_fun2m, d_fundiag, errest
    data = zeros(length(fncts), length(matrices), 11);

    for i = 1 : length(fncts)
        for j = 1 : length(matrices);
            matrix_name = matrices[j];

            fname = fncts[i]["name"];
            f = fncts[i]["f"];

            vars = matread("$(data_directory)/experiment2_$(matrix_name)_$(fname).mat")
            A = vars["A"]; B = vars["B"]; C = vars["C"]; Y = vars["Y"]; errest = vars["errest"]

            Y2, infom = fun2m(f, A, B, C, method = BivMatFun.Diag, bs = max(m,n));

            if ! is_ci_test()
                tdiagm = @benchmark fun2m($f, $A, $B, $C, method = $BivMatFun.Diag, bs = $max($m, $n))
                tdiagm = median(tdiagm.times) / 1e9;
            else
                tdiagm = 0.0
            end

            XD = BivMatFun.diag_fun(f, A, B, C);
            if ! is_ci_test()
                tdiag = @benchmark BivMatFun.diag_fun($f, $A, $B, $C)
                tdiag = median(tdiag.times) / 1e9;
            else
                tdiag = 0.0
            end

            X, info = fun2m(f, A, B, C, method = BivMatFun.Diag);
            if ! is_ci_test()
                tsp = @benchmark fun2m($f, $A, $B, $C, method = BivMatFun.Diag);
                tsp = median(tsp.times) / 1e9;
            else
                tsp = 0.0
            end

            #% err_fun2m, err_diag, time_fun2m, time_diag, time_diagm, nblocksA, nblocksB, d_fun2m, d_fundiagm, errest
            data[i, j, 1] = norm(X - Y) / norm(Y);
            data[i, j, 2] = norm(XD - Y) / norm(Y);
            data[i, j, 3] = norm(Y2 - Y) / norm(Y);
            data[i, j, 4] = tsp;
            data[i, j, 5] = tdiag;
            data[i, j, 6] = tdiagm;
            data[i, j, 7] = info.nblocksA;
            data[i, j, 8] = info.nblocksB;
            data[i, j, 9] = info.digits;
            data[i, j, 10] = infom.digits;
            data[i, j, 11] = errest;

            # In these tests, it is expected that fun2m with Diag should work in all cases, and that
            # the method using Taylor may fail. Hence, we only check if the residue is sufficiently
            # for the case with Diag; we allow some tolerance, because the error estimate is rather
            # rough.
            success &= data[i, j, 1] < errest * 1e2;

            # writedlm("experiment2.dat", data, '\t');

            @printf("%s: %s: |X - Y|/|Y|: %e (nblocks = %d, %d, time = %fs, digits = %d); errest = %e\n",
            fname, matrix_name, norm(X - Y) / norm(Y), info.nblocksA, info.nblocksB, tsp, info.digits, errest);
            @printf("%s: %s: |Y2 - Y| / |Y|: %e, time = %fs\n", fname, matrix_name, norm(Y2 - Y) / norm(Y), tdiagm)
            @printf("%s: %s: |XD - Y|/|Y|: %e, time = %fs\n", fname, matrix_name, norm(XD - Y) / norm(Y), tdiag);
        end

        # We format the data in a way that is easier to handle and ready to include
        # in the paper TeX files.
        writedlm("exp2-$i.dat", data[i,:,:], '\t')
    end

    return success;

end

run_test();

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
function xfer(s, xfer_exp)
    xfer_size = 10^xfer_exp
    xfer_block = 10^(xfer_exp - 4)

    @sync begin
        @async begin
            # read in chunks of xfer_block
            bread = 0
            while bread < xfer_size
                data = read(s, Uint8, xfer_block)
                @assert length(data) == xfer_block
                bread = bread + xfer_block
            end

#            println("process $(myid()) received $bread bytes")

        end

        @async begin
            # write in chunks of xfer_block
            data = fill!(zeros(Uint8, xfer_block), int8(65))
            bwritten = 0
            while bwritten < xfer_size
                write(s, data)
                bwritten = bwritten + xfer_block
            end
#            println("process $(myid()) sent $bwritten bytes")
        end
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#TODO: I want strings, so then I can do a string macro to return the colortype
crayons = Dict(
    "Almond"=> "#EFDECD",
    "Antique Brass"=> "#CD9575",
    "Apricot"=> "#FDD9B5",
    "Aquamarine"=> "#78DBE2",
    "Asparagus"=> "#87A96B",
    "Atomic Tangerine"=> "#FFA474",
    "Banana Mania"=> "#FAE7B5",
    "Beaver"=> "#9F8170",
    "Bittersweet"=> "#FD7C6E",
    "Black"=> "#000000",
    "Blue"=> "#1F75FE",
    "Blue Bell"=> "#A2A2D0",
    "Blue Green"=> "#0D98BA",
    "Blue Violet"=> "#7366BD",
    "Blush"=> "#DE5D83",
    "Brick Red"=> "#CB4154",
    "Brown"=> "#B4674D",
    "Burnt Orange"=> "#FF7F49",
    "Burnt Sienna"=> "#EA7E5D",
    "Cadet Blue"=> "#B0B7C6",
    "Canary"=> "#FFFF99",
    "Caribbean Green"=> "#00CC99",
    "Carnation Pink"=> "#FFAACC",
    "Cerise"=> "#DD4492",
    "Cerulean"=> "#1DACD6",
    "Chestnut"=> "#BC5D58",
    "Copper"=> "#DD9475",
    "Cornflower"=> "#9ACEEB",
    "Cotton Candy"=> "#FFBCD9",
    "Dandelion"=> "#FDDB6D",
    "Denim"=> "#2B6CC4",
    "Desert Sand"=> "#EFCDB8",
    "Eggplant"=> "#6E5160",
    "Electric Lime"=> "#CEFF1D",
    "Fern"=> "#71BC78",
    "Forest Green"=> "#6DAE81",
    "Fuchsia"=> "#C364C5",
    "Fuzzy Wuzzy"=> "#CC6666",
    "Gold"=> "#E7C697",
    "Goldenrod"=> "#FCD975",
    "Granny Smith Apple"=> "#A8E4A0",
    "Gray"=> "#95918C",
    "Green"=> "#1CAC78",
    "Green Yellow"=> "#F0E891",
    "Hot Magenta"=> "#FF1DCE",
    "Inchworm"=> "#B2EC5D",
    "Indigo"=> "#5D76CB",
    "Jazzberry Jam"=> "#CA3767",
    "Jungle Green"=> "#3BB08F",
    "Laser Lemon"=> "#FEFE22",
    "Lavender"=> "#FCB4D5",
    "Macaroni and Cheese"=> "#FFBD88",
    "Magenta"=> "#F664AF",
    "Mahogany"=> "#CD4A4C",
    "Manatee"=> "#979AAA",
    "Mango Tango"=> "#FF8243",
    "Maroon"=> "#C8385A",
    "Mauvelous"=> "#EF98AA",
    "Melon"=> "#FDBCB4",
    "Midnight Blue"=> "#1A4876",
    "Mountain Meadow"=> "#30BA8F",
    "Navy Blue"=> "#1974D2",
    "Neon Carrot"=> "#FFA343",
    "Olive Green"=> "#BAB86C",
    "Orange"=> "#FF7538",
    "Orchid"=> "#E6A8D7",
    "Outer Space"=> "#414A4C",
    "Outrageous Orange"=> "#FF6E4A",
    "Pacific Blue"=> "#1CA9C9",
    "Peach"=> "#FFCFAB",
    "Periwinkle"=> "#C5D0E6",
    "Piggy Pink"=> "#FDDDE6",
    "Pine Green"=> "#158078",
    "Pink Flamingo"=> "#FC74FD",
    "Pink Sherbert"=> "#F78FA7",
    "Plum"=> "#8E4585",
    "Purple Heart"=> "#7442C8",
    "Purple Mountains Majesty"=> "#9D81BA",
    "Purple Pizzazz"=> "#FE4EDA",
    "Radical Red"=> "#FF496C",
    "Raw Sienna"=> "#D68A59",
    "Razzle Dazzle Rose"=> "#FF48D0",
    "Razzmatazz"=> "#E3256B",
    "Red"=> "#EE204D",
    "Red Orange"=> "#FF5349",
    "Red Violet"=> "#C0448F",
    "Robin's Egg Blue"=> "#1FCECB",
    "Royal Purple"=> "#7851A9",
    "Salmon"=> "#FF9BAA",
    "Scarlet"=> "#FC2847",
    "Screamin' Green"=> "#76FF7A",
    "Sea Green"=> "#93DFB8",
    "Sepia"=> "#A5694F",
    "Shadow"=> "#8A795D",
    "Shamrock"=> "#45CEA2",
    "Shocking Pink"=> "#FB7EFD",
    "Silver"=> "#CDC5C2",
    "Sky Blue"=> "#80DAEB",
    "Spring Green"=> "#ECEABE",
    "Sunglow"=> "#FFCF48",
    "Sunset Orange"=> "#FD5E53",
    "Tan"=> "#FAA76C",
    "Tickle Me Pink"=> "#FC89AC",
    "Timberwolf"=> "#DBD7D2",
    "Tropical Rain Forest"=> "#17806D",
    "Tumbleweed"=> "#DEAA88",
    "Turquoise Blue"=> "#77DDE7",
    "Unmellow Yellow"=> "#FFFF66",
    "Violet (Purple)"=> "#926EAE",
    "Violet Red"=> "#F75394",
    "Vivid Tangerine"=> "#FFA089",
    "Vivid Violet"=> "#8F509D",
    "White"=> "#FFFFFF",
    "Wild Blue Yonder"=> "#A2ADD0",
    "Wild Strawberry"=> "#FF43A4",
    "Wild Watermelon"=> "#FC6C85",
    "Wisteria"=> "#CDA4DE",
    "Yellow"=> "#FCE883",
    "Yellow Green"=> "#C5E384",
    "Yellow Orange"=> "#FFAE42"
)

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
abstract type AbstractStructFRSpace <: KitBase.AbstractStructPhysicalSpace end
abstract type AbstractUnstructFRSpace <: KitBase.AbstractUnstructPhysicalSpace end

"""
    struct FRPSpace1D{
        A,
        I<:Integer,
        B<:AbstractVector{<:AbstractFloat},
        C<:AbstractMatrix{<:AbstractFloat},
    } <: AbstractStructFRSpace
        base::A
        deg::I
        J::B
        np::I
        xpl::B
        xpg::C
        wp::B
        dl::C
        ll::B
        lr::B
        dhl::B
        dhr::B
    end

1D physical space for flux reconstruction method

"""
struct FRPSpace1D{
    A,
    I<:Integer,
    B<:AbstractVector{<:AbstractFloat},
    C<:AbstractMatrix{<:AbstractFloat},
} <: AbstractStructFRSpace
    base::A

    deg::I
    J::B
    np::I

    xpl::B
    xpg::C
    wp::B

    dl::C
    ll::B
    lr::B
    dll::B
    dlr::B
    dhl::B
    dhr::B
    V::C
    iV::C
end

function FRPSpace1D(x0::Real, x1::Real, nx::Integer, deg::Integer, ng = 0::Integer)
    ps = PSpace1D(x0, x1, nx, ng)
    J = [ps.dx[i] / 2 for i in eachindex(ps.dx)]

    r = legendre_point(deg) .|> eltype(ps.x)
    xi = push!(ps.x - 0.5 * ps.dx, ps.x[end] + 0.5 * ps.dx[end]) .|> eltype(ps.x)
    xp = global_sp(xi, r)
    wp = gausslegendre(deg + 1)[2]

    ll = lagrange_point(r, -1.0)
    lr = lagrange_point(r, 1.0)
    lpdm = ∂lagrange(r)

    V = vandermonde_matrix(deg, r) |> Array
    iV = inv(V)
    dVf = ∂vandermonde_matrix(deg, [-1.0, 1.0])
    ∂lf = zeros(2, deg + 1)
    for i = 1:2
        ∂lf[i, :] .= V' \ dVf[i, :]
    end
    dll = ∂lf[1, :]
    dlr = ∂lf[2, :]

    dhl, dhr = ∂radau(deg, r)

    return FRPSpace1D{typeof(ps),typeof(deg),typeof(J),typeof(xp)}(
        ps,
        deg,
        J,
        deg + 1,
        r,
        xp,
        wp,
        lpdm,
        ll,
        lr,
        dll,
        dlr,
        dhl,
        dhr,
        V,
        iV,
    )
end

"""

2D physical space for flux reconstruction method

"""
struct FRPSpace2D{
    A,
    I<:Integer,
    B,
    C<:AbstractVector{<:AbstractFloat},
    D<:AbstractArray{<:AbstractFloat,5},
    E<:AbstractMatrix{<:AbstractFloat},
} <: AbstractStructFRSpace
    base::A

    deg::I
    J::B
    iJ::B
    Ji::B
    np::I

    xpl::C
    xpg::D
    wp::E

    dl::E
    ll::C
    lr::C
    dll::C
    dlr::C
    dhl::C
    dhr::C
    V::E
    iV::E
end

function FRPSpace2D(base::AbstractPhysicalSpace2D, deg::Integer)
    # solution points in 1D standard element
    r = legendre_point(deg) .|> eltype(base.x)

    # Jacobian and its inverse
    J = rs_jacobi(r, base.vertices)

    iJ = deepcopy(J)
    for i in axes(iJ, 1), j in axes(iJ, 2)
        for k = 1:deg+1, l = 1:deg+1
            iJ[i, j][k, l] .= inv(J[i, j][k, l])
        end
    end

    # interface flux points
    ri = zeros(eltype(base.x), 4, deg + 1)
    ri[1, :] .= r
    ri[2, :] .= 1.0
    ri[3, :] .= r[end:-1:1]
    ri[4, :] .= 0.0

    si = zeros(eltype(base.x), 4, deg + 1)
    si[1, :] .= 0.0
    si[2, :] .= r
    si[3, :] .= 1.0
    si[4, :] .= r[end:-1:1]

    # Jacobian of flux points
    Ji = rs_jacobi(ri, si, base.vertices)

    # solution points in global coordinates
    xpg = OffsetArray{eltype(base.x)}(
        undef,
        axes(base.x, 1),
        axes(base.y, 2),
        1:deg+1,
        1:deg+1,
        1:2,
    )
    for i in axes(xpg, 1), j in axes(xpg, 2), k = 1:deg+1, l = 1:deg+1
        @. xpg[i, j, k, l, :] =
            (r[k] - 1.0) * (r[l] - 1.0) / 4 * base.vertices[i, j, 1, :] +
            (r[k] + 1.0) * (1.0 - r[l]) / 4 * base.vertices[i, j, 2, :] +
            (r[k] + 1.0) * (r[l] + 1.0) / 4 * base.vertices[i, j, 3, :] +
            (1.0 - r[k]) * (r[l] + 1.0) / 4 * base.vertices[i, j, 4, :]
    end

    # quadrature weights
    w = gausslegendre(deg + 1)[2] .|> eltype(base.x)
    wp = [w[i] * w[j] for i = 1:deg+1, j = 1:deg+1]

    # Lagrange polynomials
    ll, lr, lpdm = standard_lagrange(r)

    V = vandermonde_matrix(deg, r)
    dVf = ∂vandermonde_matrix(deg, [-1.0, 1.0])
    ∂lf = zeros(eltype(base.x), 2, deg + 1)
    for i = 1:2
        ∂lf[i, :] .= V' \ dVf[i, :]
    end
    dll = ∂lf[1, :]
    dlr = ∂lf[2, :]

    dhl, dhr = ∂radau(deg, r)

    # 2D Vandermonde matrix
    r2d = zeros(eltype(r), length(r), length(r))
    for j in axes(r2d, 2)
        r2d[:, j] .= r
    end
    s2d = deepcopy(r2d)
    for i in axes(s2d, 1)
        s2d[i, :] .= r
    end
    V = vandermonde_matrix(Quad, deg, r2d[:], s2d[:])
    iV = inv(V)

    return FRPSpace2D{typeof(base),typeof(deg),typeof(J),typeof(r),typeof(xpg),typeof(wp)}(
        base,
        deg,
        J,
        iJ,
        Ji,
        (deg + 1)^2,
        r,
        xpg,
        wp,
        lpdm,
        ll,
        lr,
        dll,
        dlr,
        dhl,
        dhr,
        V,
        iV,
    )
end

function FRPSpace2D(
    x0::Real,
    x1::Real,
    nx::Integer,
    y0::Real,
    y1::Real,
    ny::Integer,
    deg::Integer,
    ngx = 0::Integer,
    ngy = 0::Integer,
)
    ps = PSpace2D(x0, x1, nx, y0, y1, ny, ngx, ngy)
    return FRPSpace2D(ps, deg)
end

"""
Unstructued physical space for flux reconstruction method

"""
struct UnstructFRPSpace{
    A,
    H,
    G<:Integer,
    B<:AbstractMatrix{<:AbstractFloat},
    F<:AbstractArray{<:AbstractFloat,3},
    E<:AbstractVector{<:AbstractFloat},
    I<:AbstractArray{<:AbstractFloat,4},
    J,
} <: AbstractUnstructFRSpace
    #--- general ---#
    base::A # basic unstructured mesh info that contains:
    #=
    cells::A # all information: cell, line, vertex
    points::B # locations of vertex points
    cellid::C # node indices of elements
    cellType::D # inner/boundary cell
    cellNeighbors::C # neighboring cells id
    cellFaces::C # cell edges id
    cellCenter::B # cell center location
    cellArea::E # cell size
    cellNormals::F # cell unit normal vectors
    facePoints::C # ids of two points at edge
    faceCells::C # ids of two cells around edge
    faceCenter::B # edge center location
    faceType::D # inner/boundary face
    faceArea::E # face area
    =#

    #--- FR specific ---#
    J::H # Jacobi
    deg::G # polynomial degree
    np::G # number of solution points
    xpl::B # local coordinates of solution points
    xpg::F # global coordinates of solution points
    wp::E # weights of solution points
    xfl::F # local coordinates of flux points
    xfg::I # global coordinates of flux points
    wf::B # weights of flux points
    V::B # Vandermonde matrix
    ψf::F # Vandermonde matrix along faces
    Vr::B # ∂V/∂r
    Vs::B # ∂V/∂s
    ∂l::F # ∇l
    lf::F # Lagrange polynomials along faces
    ϕ::F # correction field
    fpn::J # adjacent flux points index in neighbor cell
end

function TriFRPSpace(file::T, deg::Integer) where {T<:AbstractString}
    ps = UnstructPSpace(file)

    J = rs_jacobi(ps.cellid, ps.points)
    np = (deg + 1) * (deg + 2) ÷ 2
    xpl, wp = tri_quadrature(deg)
    V = vandermonde_matrix(Tri, deg, xpl[:, 1], xpl[:, 2])
    Vr, Vs = ∂vandermonde_matrix(Tri, deg, xpl[:, 1], xpl[:, 2])
    ∂l = ∂lagrange(V, Vr, Vs)
    ϕ = correction_field(deg, V)

    xfl, wf = triface_quadrature(deg)
    ψf = zeros(3, deg + 1, np)
    for i = 1:3
        ψf[i, :, :] .= vandermonde_matrix(Tri, deg, xfl[i, :, 1], xfl[i, :, 2])
    end

    lf = zeros(3, deg + 1, np)
    for i = 1:3, j = 1:deg+1
        lf[i, j, :] .= V' \ ψf[i, j, :]
    end

    xpg = global_sp(ps.points, ps.cellid, deg)
    xfg = global_fp(ps.points, ps.cellid, deg)
    ncell = size(ps.cellid, 1)
    fpn = [neighbor_fpidx([i, j, k], ps, xfg) for i = 1:ncell, j = 1:3, k = 1:deg+1]

    return UnstructFRPSpace(
        ps,
        J,
        deg,
        np,
        xpl,
        xpg,
        wp,
        xfl,
        xfg,
        wf,
        V,
        ψf,
        Vr,
        Vs,
        ∂l,
        lf,
        ϕ,
        fpn,
    )
end


text:
module SpringCollab2020ChangeThemeTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/ChangeThemeTrigger" ChangeThemeTrigger(x::Integer, y::Integer,
    width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, enable::Bool=true)

const placements = Ahorn.PlacementDict(
    "Change Theme Trigger (Spring Collab 2020)" => Ahorn.EntityPlacement(
        ChangeThemeTrigger,
        "rectangle",
    ),
)

end


text:
module ConformalEdge

export DMap, trace_fieldline, SOL, fluxmap

using LinearAlgebra, ForwardDiff, ConformalMaps, FFTW

include("./dmap.jl")
include("./fieldlines.jl")
include("./sol.jl")

end