text:
export
    say,              # debug output
    @verbose_say,     # define say according to verbose
    timestr,          # time string, useful for file names
    floatstr,         # short float string
    share,            # make a tuple of array to SharedArray
    load_data,        # load synthetic data
    generate_data     # generate synthetic data

""" debug output """
say = println

"""
if :verbose not defined, do nothing;
if verbose true, define say as println;
else mute say.
"""
macro verbose_say()
  return esc(quote
    isdefined(:verbose) && (say = verbose ? println : (sth...) -> nothing)
  end)
end

""" time string, useful for file names, etc. """
timestr() = "$(Dates.format(now(),"yyyymmdd-HHMMSS"))"

""" short float string """
floatstr(num, precision=1000) = "$(floor(num*precision)/precision)"

""" make a tuple of array to SharedArray """
share(t::Tuple) = [convert(SharedArray, i) for i in t]

""" load synthetic data """
function load_data(name="set2")
  file = "data/$name.jld"
  isfile(file) || error("file `$file` not found")
  d = load(file)
  say("dataset `$name` loaded.")

  (d["train_data"], d["train_label"], d["validation_data"], d["validation_label"], d["test_data"], d["test_label"])
end

""" generate synthetic data, column vectors """
function generate_data(;name=timestr(), train=25000, validation=5000, test=20000, d=50, p=2, prob=0.9)
  isdir("data") || mkdir("data")

  n = train + validation + test
  data = rand([-1, 1], d + 1, n)
  data[end, :] = 1

  p_index = sample(1:d, p, replace=false, ordered=true)
  labels = prod(data[p_index, :], 1)[:]

  mutate_index = sample(1:n, floor(Int, (1-prob)*n), replace=false, ordered=true)
  mutate = ones(Int, n)
  mutate[mutate_index] = -1
  labels = labels .* mutate

  train_index = 1 : train
  validation_index = train + 1 : train + validation
  test_index = train + validation + 1 : n

  save("data/$name.jld",  "train_data", data[:, train_index],
                          "validation_data", data[:, validation_index],
                          "test_data", data[:, test_index],
                          "train_label", labels[train_index],
                          "validation_label", labels[validation_index],
                          "test_label", labels[test_index],
                          "p_index", p_index)

  say("data generated in `data/$name.jld`")
end

max_stars_repo_name: innerlee/Halfspace

max_issues_repo_name: innerlee/Halfspace

max_forks_repo_name: innerlee/Halfspace


text:
module Credit

import ..Clarus

function saccr(;params...)
  return Clarus.api_request("Credit","SACCR";params...)
end

end #MODULE-END

ext: jl


text:
# Residual Sum of Squares
function compute_objval_rss(sim_data::Vector{Float64}, exp_data::Vector{Float64})::Float64
    error::Float64 = 0.0

    for i in eachindex(exp_data)
        @inbounds error += (sim_data[i] - exp_data[i])^2
    end

    return error
end

# Cosine similarity
function compute_objval_cos(sim_data::Vector{Float64}, exp_data::Vector{Float64})::Float64

    error::Float64 = 1.0 - dot(sim_data,exp_data)/(norm(sim_data)*norm(exp_data))

    return error
end

function conditions_index(condition_name::String)::Int

    return findfirst(isequal(condition_name),Sim.conditions)
end

function diff_sim_and_exp(sim_matrix::Matrix{Float64},exp_dict::Dict{String,Array{Float64,1}},
                            exp_timepoint::Vector{Float64},conditions::Vector{String};
                            sim_norm_max::Float64,exp_norm_max::Float64)
    sim_result::Vector{Float64} = []
    exp_result::Vector{Float64} = []

    for (i,condition) in enumerate(conditions)
        append!(sim_result,sim_matrix[Int.(exp_timepoint.+1),i])
        append!(exp_result,exp_dict[condition])
    end

    return (sim_result./sim_norm_max, exp_result./exp_norm_max)
end

# Define an objective function to be minimized.
function objective(indiv_gene::Vector{Float64})::Float64

    indiv::Vector{Float64} = decode_gene2val(indiv_gene)

    (p,u0) = update_param(indiv)

    if Sim.simulate!(p,u0) isa Nothing
        error::Vector{Float64} = zeros(18)

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_Akt"),:,[conditions_index("MCF7_EGF"),conditions_index("MCF7_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_Akt")
        error[1] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_Akt"), Int.(exp_t.+1), conditions_index("MCF7_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_Akt")]["MCF7_EGF"]
        )
        error[2] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_Akt"), Int.(exp_t.+1), conditions_index("MCF7_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_Akt")]["MCF7_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_ERK"),:,[conditions_index("MCF7_EGF"),conditions_index("MCF7_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_ERK")
        error[3] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERK"), Int.(exp_t.+1), conditions_index("MCF7_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERK")]["MCF7_EGF"]
        )
        error[4] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERK"), Int.(exp_t.+1), conditions_index("MCF7_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERK")]["MCF7_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_cFos"),:,[conditions_index("MCF7_EGF"),conditions_index("MCF7_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_cFos")
        error[5] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("MCF7_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["MCF7_EGF"]
        )
        error[6] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("MCF7_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["MCF7_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_Akt"),:,[conditions_index("BT474_EGF"),conditions_index("BT474_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_Akt")
        error[7] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_Akt"), Int.(exp_t.+1), conditions_index("BT474_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_Akt")]["BT474_EGF"]
        )
        error[8] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_Akt"), Int.(exp_t.+1), conditions_index("BT474_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_Akt")]["BT474_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_ERK"),:,[conditions_index("BT474_EGF"),conditions_index("BT474_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_ERK")
        error[9] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERK"), Int.(exp_t.+1), conditions_index("BT474_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERK")]["BT474_EGF"]
        )
        error[10] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERK"), Int.(exp_t.+1), conditions_index("BT474_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERK")]["BT474_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_cFos"),:,[conditions_index("BT474_EGF"),conditions_index("BT474_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_cFos")
        error[11] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("BT474_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["BT474_EGF"]
        )
        error[12] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("BT474_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["BT474_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_Akt"),:,[conditions_index("MDAMB231_EGF"),conditions_index("MDAMB231_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_Akt")
        error[13] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_Akt"), Int.(exp_t.+1), conditions_index("MDAMB231_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_Akt")]["MDAMB231_EGF"]
        )
        error[14] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_Akt"), Int.(exp_t.+1), conditions_index("MDAMB231_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_Akt")]["MDAMB231_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_ERK"),:,[conditions_index("MDAMB231_EGF"),conditions_index("MDAMB231_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_ERK")
        error[15] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERK"), Int.(exp_t.+1), conditions_index("MDAMB231_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERK")]["MDAMB231_EGF"]
        )
        error[16] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERK"), Int.(exp_t.+1), conditions_index("MDAMB231_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERK")]["MDAMB231_HRG"]
        )

        sim_norm_max = maximum(
            Sim.simulations[observables_index("Phosphorylated_cFos"),:,[conditions_index("MDAMB231_EGF"),conditions_index("MDAMB231_HRG")]]
        )
        exp_t = Exp.get_timepoint("Phosphorylated_cFos")
        error[17] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("MDAMB231_EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["MDAMB231_EGF"]
        )
        error[18] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("MDAMB231_HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["MDAMB231_HRG"]
        )

        return sum(error)
    else
        return NaN
    end
end


text:
# Julia wrapper for header: /usr/include/expat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Void XML_SetElementDeclHandler (XML_Parser, XML_ElementDeclHandler) libexpat
@c Void XML_SetAttlistDeclHandler (XML_Parser, XML_AttlistDeclHandler) libexpat
@c Void XML_SetXmlDeclHandler (XML_Parser, XML_XmlDeclHandler) libexpat
@c XML_Parser XML_ParserCreate (Ptr{XML_Char},) libexpat
@c XML_Parser XML_ParserCreateNS (Ptr{XML_Char}, XML_Char) libexpat
@c XML_Parser XML_ParserCreate_MM (Ptr{XML_Char}, Ptr{XML_Memory_Handling_Suite}, Ptr{XML_Char}) libexpat
@c XML_Bool XML_ParserReset (XML_Parser, Ptr{XML_Char}) libexpat
@c Void XML_SetEntityDeclHandler (XML_Parser, XML_EntityDeclHandler) libexpat
@c Void XML_SetElementHandler (XML_Parser, XML_StartElementHandler, XML_EndElementHandler) libexpat
@c Void XML_SetStartElementHandler (XML_Parser, XML_StartElementHandler) libexpat
@c Void XML_SetEndElementHandler (XML_Parser, XML_EndElementHandler) libexpat
@c Void XML_SetCharacterDataHandler (XML_Parser, XML_CharacterDataHandler) libexpat
@c Void XML_SetProcessingInstructionHandler (XML_Parser, XML_ProcessingInstructionHandler) libexpat
@c Void XML_SetCommentHandler (XML_Parser, XML_CommentHandler) libexpat
@c Void XML_SetCdataSectionHandler (XML_Parser, XML_StartCdataSectionHandler, XML_EndCdataSectionHandler) libexpat
@c Void XML_SetStartCdataSectionHandler (XML_Parser, XML_StartCdataSectionHandler) libexpat
@c Void XML_SetEndCdataSectionHandler (XML_Parser, XML_EndCdataSectionHandler) libexpat
@c Void XML_SetDefaultHandler (XML_Parser, XML_DefaultHandler) libexpat
@c Void XML_SetDefaultHandlerExpand (XML_Parser, XML_DefaultHandler) libexpat
@c Void XML_SetDoctypeDeclHandler (XML_Parser, XML_StartDoctypeDeclHandler, XML_EndDoctypeDeclHandler) libexpat
@c Void XML_SetStartDoctypeDeclHandler (XML_Parser, XML_StartDoctypeDeclHandler) libexpat
@c Void XML_SetEndDoctypeDeclHandler (XML_Parser, XML_EndDoctypeDeclHandler) libexpat
@c Void XML_SetUnparsedEntityDeclHandler (XML_Parser, XML_UnparsedEntityDeclHandler) libexpat
@c Void XML_SetNotationDeclHandler (XML_Parser, XML_NotationDeclHandler) libexpat
@c Void XML_SetNamespaceDeclHandler (XML_Parser, XML_StartNamespaceDeclHandler, XML_EndNamespaceDeclHandler) libexpat
@c Void XML_SetStartNamespaceDeclHandler (XML_Parser, XML_StartNamespaceDeclHandler) libexpat
@c Void XML_SetEndNamespaceDeclHandler (XML_Parser, XML_EndNamespaceDeclHandler) libexpat
@c Void XML_SetNotStandaloneHandler (XML_Parser, XML_NotStandaloneHandler) libexpat
@c Void XML_SetExternalEntityRefHandler (XML_Parser, XML_ExternalEntityRefHandler) libexpat
@c Void XML_SetExternalEntityRefHandlerArg (XML_Parser, Ptr{Void}) libexpat
@c Void XML_SetSkippedEntityHandler (XML_Parser, XML_SkippedEntityHandler) libexpat
@c Void XML_SetUnknownEncodingHandler (XML_Parser, XML_UnknownEncodingHandler, Ptr{Void}) libexpat
@c Void XML_DefaultCurrent (XML_Parser,) libexpat
@c Void XML_SetReturnNSTriplet (XML_Parser, Int32) libexpat
@c Void XML_SetUserData (XML_Parser, Ptr{Void}) libexpat
@c Int32 XML_SetEncoding (XML_Parser, Ptr{XML_Char}) libexpat
@c Void XML_UseParserAsHandlerArg (XML_Parser,) libexpat
@c Int32 XML_UseForeignDTD (XML_Parser, XML_Bool) libexpat
@c Int32 XML_SetBase (XML_Parser, Ptr{XML_Char}) libexpat
@c Ptr{XML_Char} XML_GetBase (XML_Parser,) libexpat
@c Int32 XML_GetSpecifiedAttributeCount (XML_Parser,) libexpat
@c Int32 XML_GetIdAttributeIndex (XML_Parser,) libexpat
@c Int32 XML_Parse (XML_Parser, Ptr{UInt8}, Int32, Int32) libexpat
@c Ptr{Void} XML_GetBuffer (XML_Parser, Int32) libexpat
@c Int32 XML_ParseBuffer (XML_Parser, Int32, Int32) libexpat
@c Int32 XML_StopParser (XML_Parser, XML_Bool) libexpat
@c Int32 XML_ResumeParser (XML_Parser,) libexpat
@c Void XML_GetParsingStatus (XML_Parser, Ptr{XML_ParsingStatus}) libexpat
@c XML_Parser XML_ExternalEntityParserCreate (XML_Parser, Ptr{XML_Char}, Ptr{XML_Char}) libexpat
@c Int32 XML_SetParamEntityParsing (XML_Parser, Void) libexpat
@c Int32 XML_SetHashSalt (XML_Parser, UInt32) libexpat
@c Int32 XML_GetErrorCode (XML_Parser,) libexpat
@c XML_Size XML_GetCurrentLineNumber (XML_Parser,) libexpat
@c XML_Size XML_GetCurrentColumnNumber (XML_Parser,) libexpat
@c XML_Index XML_GetCurrentByteIndex (XML_Parser,) libexpat
@c Int32 XML_GetCurrentByteCount (XML_Parser,) libexpat
@c Ptr{UInt8} XML_GetInputContext (XML_Parser, Ptr{Int32}, Ptr{Int32}) libexpat
@c Void XML_FreeContentModel (XML_Parser, Ptr{XML_Content}) libexpat
@c Ptr{Void} XML_MemMalloc (XML_Parser, size_t) libexpat
@c Ptr{Void} XML_MemRealloc (XML_Parser, Ptr{Void}, size_t) libexpat
@c Void XML_MemFree (XML_Parser, Ptr{Void}) libexpat
@c Void XML_ParserFree (XML_Parser,) libexpat
@c Ptr{XML_LChar} XML_ErrorString (Void,) libexpat
@c Ptr{XML_LChar} XML_ExpatVersion () libexpat
@c XML_Expat_Version XML_ExpatVersionInfo () libexpat
@c Ptr{XML_Feature} XML_GetFeatureList () libexpat


text:
using Keldysh, LinearAlgebra, Test

@testset "storage" begin
  let norb=2, N=10
    X = rand(norb, norb, N, N)
    Y = rand(norb, norb, N, N)
    Z = X .+ 2 .* Y .- 3 .* X

    A = Keldysh.GenericStorage(X)
    B = Keldysh.GenericStorage(Y)
    C = A + 2*B - 3*A

    D = similar(A)
    for (i,j) in Iterators.product(1:N,1:N)
      D[i,j] = Z[:,:,i,j]
    end

    E = zero(A)
    for (k,l,i,j) in Iterators.product(1:norb, 1:norb, 1:N, 1:N)
      E[k,l,i,j] = Z[k,l,i,j]
    end

    @test size(A) == (norb, norb, N, N)

    # test getindex
    for (i,j) in Iterators.product(1:N,1:N)
      @test C[i,j] == Z[:,:,i,j]
      @test D[i,j] == Z[:,:,i,j]
      @test E[i,j] == Z[:,:,i,j]
    end
    for (k,l,i,j) in Iterators.product(1:norb, 1:norb, 1:N, 1:N)
      @test C[k,l,i,j] == Z[k,l,i,j]
      @test D[k,l,i,j] == Z[k,l,i,j]
      @test E[k,l,i,j] == Z[k,l,i,j]
    end
  end

  let norb=2, N=10
    get_rand_herm = () -> begin
      X = rand(ComplexF64, norb, norb, N, N)
      X = im * (X + conj(permutedims(X, (2,1,4,3))))
      return X
    end

    X = get_rand_herm()
    Y = get_rand_herm()

    Z = X .+ 2 .* Y .- 3 .* X

    A = Keldysh.AntiHermitianStorage(X)
    B = Keldysh.AntiHermitianStorage(Y)
    C = A + 2*B - 3*A

    Dup = similar(A)
    Dlo = similar(A)
    for (i,j) in Iterators.product(1:N, 1:N)
      i >= j && (Dup[i,j] = Z[:,:,i,j])
      i <= j && (Dlo[i,j] = Z[:,:,i,j])
    end

    Eup = zero(A)
    Elo = zero(A)
    for (k,l,i,j) in Iterators.product(1:norb, 1:norb, 1:N, 1:N)
      i >= j && (Eup[k,l,i,j] = Z[k,l,i,j])
      i <= j && (Elo[k,l,i,j] = Z[k,l,i,j])
    end

    for (i,j) in Iterators.product(1:N, 1:N)
      @test C[i,j] == Z[:,:,i,j]
      @test Dup[i,j] == Z[:,:,i,j]
      @test Dlo[i,j] == Z[:,:,i,j]
      @test Eup[i,j] == Z[:,:,i,j]
      @test Elo[i,j] == Z[:,:,i,j]
    end

    for (k,l,i,j) in Iterators.product(1:norb, 1:norb, 1:N, 1:N)
      @test C[k,l,i,j] == Z[k,l,i,j]
      @test Dup[k,l,i,j] == Z[k,l,i,j]
      @test Dlo[k,l,i,j] == Z[k,l,i,j]
      @test Eup[k,l,i,j] == Z[k,l,i,j]
      @test Elo[k,l,i,j] == Z[k,l,i,j]
    end
  end

  let norb=2, N=10
    get_rand_herm_toeplitz = () -> begin
      v = rand(ComplexF64, norb, norb, N)
      v[:,:,1] = im*(v[:,:,1] + v[:,:,1]')
      X = zeros(ComplexF64, norb, norb, N, N)
      for (i,j) in Iterators.product(1:N, 1:N)
        X[:,:,i,j] = i >= j ? v[:,:,i-j+1] : -v[:,:,j-i+1]'
      end
      return X
    end

    X = get_rand_herm_toeplitz()
    Y = get_rand_herm_toeplitz()

    Z = X .+ 2 .* Y .- 3 .* X

    A = Keldysh.AntiHermitianToeplitzStorage(X)
    B = Keldysh.AntiHermitianToeplitzStorage(Y)
    C = A + 2*B - 3*A

    Dup = similar(A)
    Dlo = similar(A)
    for i in 1:N
      Dup[i,1] = Z[:,:,i,1]
      Dlo[1,i] = Z[:,:,1,i]
    end

    Eup = zero(A)
    Elo = zero(A)
    for (k,l,i) in Iterators.product(1:norb, 1:norb, 1:N)
      Eup[k,l,i,1] = Z[k,l,i,1]
      Elo[k,l,1,i] = Z[k,l,1,i]
    end

    for (i,j) in Iterators.product(1:N, 1:N)
      @test C[i,j] == Z[:,:,i,j]
      @test Dup[i,j] == Z[:,:,i,j]
      @test Dlo[i,j] == Z[:,:,i,j]
      @test Eup[i,j] == Z[:,:,i,j]
      @test Elo[i,j] == Z[:,:,i,j]
    end

    for (k,l,i,j) in Iterators.product(1:norb, 1:norb, 1:N, 1:N)
      @test C[k,l,i,j] == Z[k,l,i,j]
      @test Dup[k,l,i,j] == Z[k,l,i,j]
      @test Dlo[k,l,i,j] == Z[k,l,i,j]
      @test Eup[k,l,i,j] == Z[k,l,i,j]
      @test Elo[k,l,i,j] == Z[k,l,i,j]
    end
  end

  let norb=2, N=10
    get_rand_periodic = () -> begin
      v = rand(norb, norb, N)
      X = zeros(norb, norb, N, N)
      for (i,j) in Iterators.product(1:N, 1:N)
        X[:,:,i,j] = i >= j ? v[:,:,i-j+1] : v[:,:,i-j+N]
      end
      return X
    end

    X = get_rand_periodic()
    Y = get_rand_periodic()
    Z = X .+ 2 .* Y .- 3 .* X

    A = Keldysh.PeriodicStorage(X[:,:,:,1])
    B = Keldysh.PeriodicStorage(Y[:,:,:,1])
    C = A + 2*B - 3*A

    Dup = similar(A)
    for i in 1:N
      Dup[i,1] = Z[:,:,i,1]
    end

    Eup = zero(A)
    for (k,l,i) in Iterators.product(1:norb, 1:norb, 1:N)
      Eup[k,l,i,1] = Z[k,l,i,1]
    end

    for (i,j) in Iterators.product(1:N, 1:N)
      @test C[i,j] == Z[:,:,i,j]
      @test Dup[i,j] == Z[:,:,i,j]
      @test Eup[i,j] == Z[:,:,i,j]
    end

    for (k,l,i,j) in Iterators.product(1:norb, 1:norb, 1:N, 1:N)
      @test C[k,l,i,j] == Z[k,l,i,j]
      @test Dup[k,l,i,j] == Z[k,l,i,j]
      @test Eup[k,l,i,j] == Z[k,l,i,j]
    end

  end
end


text:
using OMEinsum
using Graphs
using Test, Random
using SparseArrays
using OMEinsumContractionOrders
using KaHyPar

@testset "graph coarse graining" begin
    Random.seed!(2)
    adj = zeros(5, 6)
    for ind in [[1,1], [1,3], [2,2], [3,2], [4,4], [4,5], [4,3], [5, 4]]
        adj[ind...] = 1
    end
    parts = [[1,3], [2], [4]]
    incidence_list = OMEinsumContractionOrders.get_coarse_grained_graph(sparse(adj), parts)
    @test incidence_list.v2e[1] == [1,2,3]
    @test incidence_list.v2e[2] == [2]
    @test incidence_list.v2e[3] == [3,4,5]
    @test incidence_list.openedges == [4]

    @test length(OMEinsumContractionOrders._connected_components(adj, parts[1])) == 2

    res = OMEinsumContractionOrders.coarse_grained_optimize(adj, parts, ones(6), GreedyMethod(OMEinsum.MinSpaceOut(), 10))
    @test res == OMEinsum.ContractionTree(OMEinsum.ContractionTree(1,2), 3)
    @test OMEinsumContractionOrders.map_tree_to_parts(res, [[[1,2], 3], [7,6], [9, [4,1]]]) == [[[[1,2], 3], [7,6]], [9, [4,1]]]
end

@testset "kahypar" begin
    Random.seed!(2)
    function random_regular_eincode(n, k)
        g = Graphs.random_regular_graph(n, k)
        ixs = [minmax(e.src,e.dst) for e in Graphs.edges(g)]
        return EinCode((ixs..., [(i,) for i in Graphs.vertices(g)]...), ())
    end

    g = random_regular_graph(220, 3)
    rows = Int[]
    cols = Int[]
    for (i,edge) in enumerate(edges(g))
        push!(rows, edge.src, edge.dst)
        push!(cols, i, i)
    end
    graph = sparse(rows, cols, ones(Int, length(rows)))
    sc_target = 28.0
    log2_sizes = fill(1, size(graph, 2))
    b = KaHyParBipartite(sc_target=sc_target, imbalances=[0.0:0.02:0.8...])
    group1, group2 = OMEinsumContractionOrders.bipartite_sc(b, graph, collect(1:size(graph, 1)), log2_sizes)
    @test OMEinsumContractionOrders.group_sc(graph, group1, log2_sizes) <= sc_target
    @test OMEinsumContractionOrders.group_sc(graph, group2, log2_sizes) <= sc_target
    sc_target = 27.0
    group11, group12 = OMEinsumContractionOrders.bipartite_sc(b, graph, group1, log2_sizes)
    @test OMEinsumContractionOrders.group_sc(graph, group11, log2_sizes) <= sc_target
    @test OMEinsumContractionOrders.group_sc(graph, group12, log2_sizes) <= sc_target

    code = random_regular_eincode(220, 3)
    res = optimize_kahypar(code,uniformsize(code, 2); max_group_size=50, sc_target=30)
    tc, sc = OMEinsum.timespace_complexity(res, uniformsize(code, 2))
    @test sc <= 30

    # contraction test
    code = random_regular_eincode(50, 3)
    codeg = optimize_kahypar(code, uniformsize(code, 2); max_group_size=10, sc_target=12)
    codek = optimize_greedy(code, uniformsize(code, 2))
    tc, sc = OMEinsum.timespace_complexity(codek, uniformsize(code, 2))
    @test sc <= 12
    xs = [[2*randn(2, 2) for i=1:75]..., [randn(2) for i=1:50]...]
    resg = codeg(xs...)
    resk = codek(xs...)
    @test resg ≈ resk

    Random.seed!(2)
    code = random_regular_eincode(220, 3)
    codeg_auto = optimize_kahypar_auto(code, uniformsize(code, 2), greedy_method=OMEinsum.MinSpaceOut())
    tc, sc = OMEinsum.timespace_complexity(codeg_auto, uniformsize(code, 2))
    @test sc <= 30
end

@testset "regression test" begin
    code = ein"i->"
    optcode = optimize_kahypar(code, Dict('i'=>4), sc_target=10, max_group_size=10)
    @test optcode isa NestedEinsum
    x = randn(4)
    @test optcode(x) ≈ code(x)

    code = ein"i,j->"
    optcode = optimize_kahypar(code, Dict('i'=>4, 'j'=>4), sc_target=10, max_group_size=10)
    @test optcode isa NestedEinsum
    x = randn(4)
    y = randn(4)
    @test optcode(x, y) ≈ code(x, y)

    code = ein"ij,jk,kl->ijl"
    optcode = optimize_kahypar(code, Dict('i'=>4, 'j'=>4, 'k'=>4, 'l'=>4), sc_target=4, max_group_size=2)
    @test optcode isa NestedEinsum
    a, b, c = [rand(4,4) for i=1:4]
    @test optcode(a, b, c) ≈ code(a, b, c)

    code = ein"ij,jk,kl->ijl"
    optcode = optimize_kahypar(code, Dict('i'=>3, 'j'=>3, 'k'=>3, 'l'=>3), sc_target=4, max_group_size=2)
    @test optcode isa NestedEinsum
    a, b, c = [rand(3,3) for i=1:4]
    @test optcode(a, b, c) ≈ code(a, b, c)
end


text:
using DataFrames
using Test

m = page_model()
include("../src/components/NonMarketDamages.jl")

nonmarketdamages = addnonmarketdamages(m)

set_param!(m, :NonMarketDamages, :rtl_realizedtemperature, readpagedata(m, "test/validationdata/rtl_realizedtemperature.csv"))
set_param!(m, :NonMarketDamages, :rcons_per_cap_MarketRemainConsumption, readpagedata(m,
"test/validationdata/rcons_per_cap_MarketRemainConsumption.csv"))
set_param!(m, :NonMarketDamages, :rgdp_per_cap_MarketRemainGDP, readpagedata(m,
"test/validationdata/rgdp_per_cap_MarketRemainGDP.csv"))
set_param!(m, :NonMarketDamages, :atl_adjustedtolerableleveloftemprise, readpagedata(m,
"test/validationdata/atl_adjustedtolerableleveloftemprise_nonmarket.csv"))
set_param!(m, :NonMarketDamages, :imp_actualreduction, readpagedata(m,
"test/validationdata/imp_actualreduction_nonmarket.csv"))
set_param!(m, :NonMarketDamages, :isatg_impactfxnsaturation, 28.333333333333336)

p = load_parameters(m)
p["y_year_0"] = 2008.
p["y_year"] = Mimi.dim_keys(m.md, :time)
set_leftover_params!(m, p)

run(m)

rcons_per_cap = m[:NonMarketDamages, :rcons_per_cap_NonMarketRemainConsumption]
rcons_per_cap_compare = readpagedata(m, "test/validationdata/rcons_per_cap_NonMarketRemainConsumption.csv")
@test rcons_per_cap ≈ rcons_per_cap_compare rtol=1e-2

rgdp_per_cap = m[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP]
rgdp_per_cap_compare = readpagedata(m, "test/validationdata/rgdp_per_cap_NonMarketRemainGDP.csv")
@test rgdp_per_cap ≈ rgdp_per_cap_compare rtol=1e-2


text:
#=
# 3D Poisson, Dirichlet bc
# Uses aMat target
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("p3damat");

useLog("p3damatlog", level=3)

generateFor("../targets/target_amat_cg.jl")

domain(3)
functionSpace(order=2)

mesh(HEXMESH, elsperdim=4)

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, "0")

# Write the weak form
coefficient("f", "-29*pi*pi*sin(2*pi*x)*sin(3*pi*y)*sin(4*pi*z)")
weakForm(u, "-dot(grad(u),grad(v)) - f*v")

solve(u);

finalize_finch()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module JuKeBOX

# Write your package code here.
using ComradeBase
using ConcreteStructs: @concrete
using Elliptic
using DocStringExtensions
using Elliptic
using LinearAlgebra
using Polyester
using StaticArrays
using StructArrays

export Observer, BAM, bam, raytrace, intensity_point

abstract type AccretionModel <: ComradeBase.AbstractPolarizedModel end
ComradeBase.visanalytic(::Type{<:AccretionModel}) = ComradeBase.NotAnalytic()
ComradeBase.imanalytic(::Type{<:AccretionModel}) = ComradeBase.IsAnalytic()

struct SimpleModel{B,G,O} <: AccretionModel
    acc::B
    g::G
    o::O
end
ComradeBase.isprimitive(::Type{<:SimpleModel}) = ComradeBase.IsPrimitive()

@inline function ComradeBase.intensity_point(s::AccretionModel, α, β)
    raytrace(-α, β, s.g, s.o, s.acc)
end

include("model_helpers.jl")
include("metric/kerr.jl")
include("bam.jl")
end


text:
struct RGscale{T<:Number}
    IR::T
    UV::T
    kspan::Tuple{T,T}
end
RGscale(UV,IR)=RGscale(IR,UV,(UV,IR))

mutable struct ODEConfig
    atol::AbstractFloat
    rtol::AbstractFloat
    dtmax::AbstractFloat
    maxiters::Int
    progress::Bool
    adaptive::Bool
    dense::Bool
    save_on::Bool
    save_start::Bool
    save_end::Bool
end

function ODEConfig(;
    atol = 1e-14,
    rtol = 1e-14,
    dtmax = 0.1,
    maxiters = 10^5,
    progress = true,
    adaptive = true,
    dense = true,
    save_on = true,
    save_start = true,
    save_end = true,
)
    ODEConfig(
        atol,
        rtol,
        dtmax,
        maxiters,
        progress,
        adaptive,
        dense,
        save_on,
        save_start,
        save_end,
    )
end


text:
export
    LineSegment1D

"""
A line segment in one dimension with a ≤ b
"""
immutable LineSegment1D
    a::Float64
    b::Float64
    LineSegment1D(a::Real, b::Real) = new(convert(Float64, min(a,b)), convert(Float64, max(a,b)))
end

Base.:+(seg::LineSegment1D, v::Real) = LineSegment1D(seg.a + v, seg.b + v)
Base.:-(seg::LineSegment1D, v::Real) = LineSegment1D(seg.a - v, seg.b - v)

"""
The distance between the line segment and the point P
"""
function get_distance(seg::LineSegment1D, P::Real)::Float64
    if P < seg.a
        return seg.a - P
    elseif P > seg.b
        return P - seg.b
    else
        return 0.0
    end
end

Base.contains(P::LineSegment1D, v::Float64) = P.a ≤ v ≤ P.b
Base.contains(P::LineSegment1D, Q::LineSegment1D) = Q.a ≥ P.a && Q.b ≤ P.b

intersects(P::LineSegment1D, Q::LineSegment1D) = P.b ≥ Q.a && Q.b ≥ P.a


text:
#Import modules
#import Pkg; Pkg.add("DataFrames")
#import Pkg; Pkg.add("JSONTables")
#import Pkg; Pkg.add("CSV")

#Using the following modules
using JSON
using DataFrames, JSONTables
using CSV

#Parse file
parsedFile = JSON.parsefile("contacts.json")
#println(parsedFile)

#Create new dictionary
userDatabase = Dict()
#Format is as such
# userDatabase = {
#   userid => {
#   "Emails" => [email1, email2],
#   "Phones" => [phone1, phone2, phone3],
#   "OrderIds" => [order1, order2],
#   "totalcontacts" => totalcontacts
#       }
#   }

function newUser(target)
    println()
    print("Adding new user ")
    print(target["OrderId"])
    userDatabase[target["OrderId"]] = Dict(
        "Emails" => [target["Email"]],
        "Phones" => [target["Phone"]],
        "OrderIds" => [target["OrderId"]],
        "TicketIds" => [target["Id"]],
        "TotalContacts" => target["Contacts"]
    )

end

function addToUser(uid, target)
    println()
    print("Updating user ")
    print(target["OrderId"])
    #Emails
    if !(target["Email"] in userDatabase[uid]["Emails"])
        push!(userDatabase[uid]["Emails"],target["Email"])
    elseif (isempty(userDatabase[uid]["Emails"]))#If not email record, start emails
        userDatabase[uid]["Emails"] = [target["Email"]]
    end

    #Phones
    if !(target["Phone"] in userDatabase[uid]["Phones"])
        push!(userDatabase[uid]["Phones"],target["Phone"])
    elseif (isempty(userDatabase[uid]["Phones"]))   #If no phone record, start array
        userDatabase[uid]["Phones"] = [target["Phone"]]
    end

    #Order Id
    if !(target["OrderId"] in userDatabase[uid]["OrderIds"])
        push!(userDatabase[uid]["OrderIds"],target["OrderId"])
    elseif (isempty(userDatabase[uid]["OrderIds"]))   #If no order record, start array
        userDatabase[uid]["OrderIds"] = [target["OrderId"]]
    end

    userDatabase[uid]["TotalContacts"]+= target["Contacts"]

end

#Search Criteria, user id refers to the individual user, target refers to the matching record
function searchCriterias(uid, target)
    #Boolean for found
    recordFound = false
    if (target["Email"] in userDatabase[uid]["Emails"])
        recordFound = true
    elseif (target["Phone"] in userDatabase[uid]["Phones"])
        recordFound = true
    elseif (target["OrderId"] in userDatabase[uid]["OrderIds"])
        recordFound = true
    end

    return recordFound
end

#Find record
function findRecord(target, searchType="Normal")
    #variable to start search
    searchBegin = true
    userFound = false
    #Get all users
    users = keys(userDatabase)
    #Loop for each user
    for (user) in users
        #Start search in target user
        findUser = searchCriterias(user, target)
        #if search ends, break out of loop
        if (findUser)
            if (searchType=="Normal")
                addToUser(user,target)
            else
                return user
            end
            userFound = true
            break
        end
    end
    if (!userFound)
        if (searchType=="Normal")
            newUser(target)
        end
    end
end

#Loop for each contact in json
for (contact) in parsedFile
    #print(contact);
    findRecord(contact)
end

println("Saving to list")
#print("Total Users: ")
#print(size(userDatabase))

#Create databank to export
exportIds = []
exportDatabank = []
for (contact) in parsedFile
    userRecords = findRecord(contact,"Look")
    ticketTrace = join(userDatabase[userRecords]["TicketIds"], "-")
    totalcontacts = userDatabase[userRecords]["TotalContacts"]
    output = ticketTrace * ", " * totalcontacts
    push!(exportIds,contact)
    push!(exportDatabank,output)
end

#Writing CSV
# new file created
touch("miaResults.csv")

# file handling in write mode
efg = open("miaResults.csv", "w")

#Store them into data frame
df = DataFrame(ticket_id = exportIds,
                ticket_trace_contact = exportDatabank
             )
println()
print("Exporting df")

CSV.write("miaResults.csv", df)

max_stars_repo_path: 01 Multi-Channel Contacts/01-mia.jl

max_issues_repo_path: 01 Multi-Channel Contacts/01-mia.jl

max_forks_repo_path: 01 Multi-Channel Contacts/01-mia.jl


text:
"""
Port for REST API on Tablo
"""
const port = ":18080/plex/"
const plex = string("http://",TabloIP,port)
const ch_ids = string(plex,"ch_ids")
const ch_info = string(plex,"ch_info?id=")
const rec_ids = string(plex,"rec_ids")
const rec_info = string(plex,"rec_info?id=")


text:
# https://github.com/bokeh/Bokeh.jl

# make255(x) = round(Int, 255 * x)

# function bokehcolor(c::Colorant)
#   @sprintf("rgba(%d, %d, %d, %1.3f)", [make255(f(c)) for f in [red,green,blue]]..., alpha(c))
# end
# bokehcolor(cs::ColorScheme) = bokehcolor(getColor(cs))

const _glyphtypes = Dict(
    :ellipse    => :Circle,
    :rect       => :Square,
    :diamond    => :Diamond,
    :utriangle  => :Triangle,
    :dtriangle  => :InvertedTriangle,
    # :pentagon   =>
    # :hexagon    =>
    # :heptagon   =>
    # :octagon    =>
    :cross      => :Cross,
    :xcross     => :X,
    :star5      => :Asterisk,
  )

function bokeh_glyph_type(d::Dict)
  lt = d[:linetype]
  mt = d[:markershape]
  if lt == :scatter && mt == :none
    mt = :ellipse
  end

  # if we have a marker, use that
  if lt == :scatter || mt != :none
    return _glyphtypes[mt]
  end

  # otherwise return a line
  return :Line
end

function get_stroke_vector(linestyle::Symbol)
  dash = 12
  dot = 3
  gap = 2
  linestyle == :solid && return Int[]
  linestyle == :dash && return Int[dash, gap]
  linestyle == :dot && return Int[dot, gap]
  linestyle == :dashdot && return Int[dash, gap, dot, gap]
  linestyle == :dashdotdot && return Int[dash, gap, dot, gap, dot, gap]
  error("unsupported linestyle: ", linestyle)
end

# ---------------------------------------------------------------------------

function _create_plot(pkg::BokehPackage; kw...)
  d = Dict(kw)

  # dumpdict(d, "plot", true)

  # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
  # TODO: initialize the plot... title, xlabel, bgcolor, etc

  datacolumns = Bokeh.BokehDataSet[]
  tools = Bokeh.tools()
  filename = tempname() * ".html"
  title = d[:title]
  w, h = d[:size]
  xaxis_type = d[:xscale] == :log10 ? :log : :auto
  yaxis_type = d[:yscale] == :log10 ? :log : :auto
  # legend = d[:legend] ? xxxx : nothing
  legend = nothing
  extra_args = Dict()  # TODO: we'll put extra settings (xlim, etc) here
  bplt = Bokeh.Plot(datacolumns, tools, filename, title, w, h, xaxis_type, yaxis_type, legend) #, extra_args)

  Plot(bplt, pkg, 0, d, Dict[])
end

function _add_series(::BokehPackage, plt::Plot; kw...)
  d = Dict(kw)

  # dumpdict(d, "plot!", true)

  bdata = Dict{Symbol, Vector}(:x => collect(d[:x]), :y => collect(d[:y]))

  glyph = Bokeh.Bokehjs.Glyph(
      glyphtype = bokeh_glyph_type(d),
      linecolor = webcolor(d[:linecolor]),  # shape's stroke or line color
      linewidth = d[:linewidth],          # shape's stroke width or line width
      fillcolor = webcolor(d[:markercolor]),
      size      = ceil(Int, d[:markersize] * 2.5),  # magic number 2.5 to keep in same scale as other backends
      dash      = get_stroke_vector(d[:linestyle])
    )

  legend = nothing  # TODO
  push!(plt.o.datacolumns, Bokeh.BokehDataSet(bdata, glyph, legend))

  push!(plt.seriesargs, d)
  plt
end

# ----------------------------------------------------------------

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot(plt::Plot{BokehPackage}, d::Dict)
end

function _update_plot_pos_size(plt::PlottingObject{BokehPackage}, d::Dict)
end

# ----------------------------------------------------------------

# accessors for x/y data

function Base.getindex(plt::Plot{BokehPackage}, i::Int)
  series = plt.o.datacolumns[i].data
  series[:x], series[:y]
end

function Base.setindex!(plt::Plot{BokehPackage}, xy::Tuple, i::Integer)
  series = plt.o.datacolumns[i].data
  series[:x], series[:y] = xy
  plt
end

# ----------------------------------------------------------------

function _add_annotations{X,Y,V}(plt::Plot{BokehPackage}, anns::AVec{@compat(Tuple{X,Y,V})})
  for ann in anns
    # TODO: add the annotation to the plot
  end
end

# ----------------------------------------------------------------

function _create_subplot(subplt::Subplot{BokehPackage}, isbefore::Bool)
  # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example

end

function _expand_limits(lims, plt::Plot{BokehPackage}, isx::Bool)
  # TODO: call expand limits for each plot data
end

function _remove_axis(plt::Plot{BokehPackage}, isx::Bool)
  # TODO: if plot is inner subplot, might need to remove ticks or axis labels
end

# ----------------------------------------------------------------

function Base.writemime(io::IO, ::MIME"image/png", plt::PlottingObject{BokehPackage})
  # TODO: write a png to io
  warn("mime png not implemented")
end

function Base.display(::PlotsDisplay, plt::Plot{BokehPackage})
  Bokeh.showplot(plt.o)
end

function Base.display(::PlotsDisplay, plt::Subplot{BokehPackage})
  # TODO: display/show the subplot
end


text:
@testset "sources" begin

@test diagnostic("abc", 0, "bad") == "bad at (0,0)\n[Before start]\n^\n"
@test diagnostic("abc", 1, "bad") == "bad at (1,1)\nabc\n^\n"
@test diagnostic("abc", 2, "bad") == "bad at (1,2)\nabc\n ^\n"
@test diagnostic("abc", 3, "bad") == "bad at (1,3)\nabc\n  ^\n"
@test diagnostic("abc", 4, "bad") == "bad at (2,0)\n[After end]\n^\n"

@test diagnostic("l1\nl2", 2, "bad") == "bad at (1,2)\nl1\n ^\n"
@test diagnostic("l1\nl2", 3, "bad") == "bad at (1,3)\nl1\n  ^\n"
@test diagnostic("l1\nl2", 4, "bad") == "bad at (2,1)\nl2\n^\n"

@test diagnostic(LineSource("abc"), LineIter(0, 0), "bad") == "bad at (0,0)\n[Not available]\n^\n"
@test diagnostic(LineSource("abc"), LineIter(1, 0), "bad") == "bad at (1,0)\n[Not available]\n^\n"
@test diagnostic(LineSource("abc"), LineIter(1, 1), "bad") == "bad at (1,1)\nabc\n^\n"
@test diagnostic(LineSource("abc"), LineIter(1, 2), "bad") == "bad at (1,2)\nabc\n ^\n"
@test diagnostic(LineSource("abc"), LineIter(1, 3), "bad") == "bad at (1,3)\nabc\n  ^\n"
@test diagnostic(LineSource("abc"), LineIter(1, 4), "bad") == "bad at (1,4)\n[Not available]\n   ^\n"
@test diagnostic(LineSource("abc"), LineIter(2, 0), "bad") == "bad at (2,0)\n[Not available]\n^\n"
@test diagnostic(LineSource("abc"), LineIter(2, 1), "bad") == "bad at (2,1)\n[Not available]\n^\n"

@test diagnostic(LineSource("l1\nl2"), LineIter(1, 2), "bad") == "bad at (1,2)\nl1\n ^\n"
@test diagnostic(LineSource("l1\nl2"), LineIter(1, 3), "bad") == "bad at (1,3)\nl1\n  ^\n"
@test diagnostic(LineSource("l1\nl2"), LineIter(2, 1), "bad") == "bad at (2,1)\nl2\n^\n"

line = Trace(p"(.|\n)+"[0:end] + Eos())
@test parse_one("abc\n", line) == ["abc\n"]
@test parse_one(LineSource("abc\n"), line) == ["abc\n"]

println("sources ok")

end


text:
mutable struct vanderPol <: SpecialDiffSys
    μ::Float64
    noise::AbstractArray
end

name(::vanderPol) = "van der Pol equation"
dimension(::vanderPol) = 2

vanderPol(; μ=0.3, noise=zeros(2)) = vanderPol(μ, noise)

function vec_field(vdp::vanderPol, x::AbstractArray)
    u, v = x
    μ = vdp.μ
    du = v
    dv = -u-μ*v*(u^2-1.0)
    return [du, dv]
end

function jacobian(vdp::vanderPol, x::AbstractArray)
    u, v = x
    μ = vdp.μ
    J = [0 1; -1-2μ*u*v -μ*(u^2-1)]
    return J
end


text:
using GeometryBasics

function vertex_interp(iso, p1::Point3{T}, p2::Point3{T}, valp1, valp2,
                       tol = 1e-5) where {T}
    if abs(iso - valp1) < tol
        return p1
    elseif abs(iso - valp2) < tol
        return p2
    elseif abs(valp1 - valp2) < tol
        return (p1 + p2) / 2
    else
        return p1 + T((iso - valp1) / (valp2 - valp1)) * (p2 - p1)
    end
end

@inline value(values::Array, points, i, j, k) = @inbounds values[i, j, k]
@inline value(values::Function, points, i, j, k) = @inbounds values(points[i, j, k])

function marching_cubes(values::Union{Function,Array}, points::Array{Point3{PT},3}, iso,
                        cube_size = 1) where {PT<:AbstractFloat}
    vertices              = Vector{Point3{PT}}()
    indices               = Vector{NTuple{3,Int}}()
    vertex_list           = zeros(Point3{PT}, 12)
    ind_pushed            = false
    edge_table, tri_table = get_edge_tri_table()
    i_max, j_max, k_max   = size(points)

    val = (i, j, k) -> value(values, points, i, j, k)
    @inbounds for k in 1:cube_size:k_max-1,
                  j in 1:cube_size:j_max-1,
                  i in 1:cube_size:i_max-1

        cube_index = val(i, j, k) < iso ? 1 : 0
        val(i + cube_size, j, k) < iso && (cube_index |= 2)
        val(i + cube_size, j + cube_size, k) < iso && (cube_index |= 4)
        val(i, j + cube_size, k) < iso && (cube_index |= 8)
        val(i, j, k + cube_size) < iso && (cube_index |= 16)
        val(i + cube_size, j, k + cube_size) < iso && (cube_index |= 32)
        val(i + cube_size, j + cube_size, k + cube_size) < iso && (cube_index |= 64)
        val(i, j + cube_size, k + cube_size) < iso && (cube_index |= 128)
        cube_index += 1
        edge_entry = edge_table[cube_index]

        if edge_entry == 0
            continue
        end
        if edge_entry & 1 == 1
            vertex_list[1] = vertex_interp(iso, points[i, j, k], points[i+cube_size, j, k],
                                           val(i, j, k), val(i + cube_size, j, k))
        end
        if edge_entry & 2 == 2
            vertex_list[2] = vertex_interp(iso, points[i+cube_size, j, k],
                                           points[i+cube_size, j+cube_size, k],
                                           val(i + cube_size, j, k),
                                           val(i + cube_size, j + cube_size, k))
        end
        if edge_entry & 4 == 4
            vertex_list[3] = vertex_interp(iso, points[i+cube_size, j+cube_size, k],
                                           points[i, j+cube_size, k],
                                           val(i + cube_size, j + cube_size, k),
                                           val(i, j + cube_size, k))
        end
        if edge_entry & 8 == 8
            vertex_list[4] = vertex_interp(iso, points[i, j+cube_size, k], points[i, j, k],
                                           val(i, j + cube_size, k), val(i, j, k))
        end
        if edge_entry & 16 == 16
            vertex_list[5] = vertex_interp(iso, points[i, j, k+cube_size],
                                           points[i+cube_size, j, k+cube_size],
                                           val(i, j, k + cube_size),
                                           val(i + cube_size, j, k + cube_size))
        end
        if edge_entry & 32 == 32
            vertex_list[6] = vertex_interp(iso, points[i+cube_size, j, k+cube_size],
                                           points[i+cube_size, j+cube_size, k+cube_size],
                                           val(i + cube_size, j, k + cube_size),
                                           val(i + cube_size, j + cube_size, k + cube_size))
        end
        if edge_entry & 64 == 64
            vertex_list[7] = vertex_interp(iso,
                                           points[i+cube_size, j+cube_size, k+cube_size],
                                           points[i, j+cube_size, k+cube_size],
                                           val(i + cube_size, j + cube_size, k + cube_size),
                                           val(i, j + cube_size, k + cube_size))
        end
        if edge_entry & 128 == 128
            vertex_list[8] = vertex_interp(iso, points[i, j+cube_size, k+cube_size],
                                           points[i, j, k+cube_size],
                                           val(i, j + cube_size, k + cube_size),
                                           val(i, j, k + cube_size))
        end
        if edge_entry & 256 == 256
            vertex_list[9] = vertex_interp(iso, points[i, j, k], points[i, j, k+cube_size],
                                           val(i, j, k), val(i, j, k + cube_size))
        end
        if edge_entry & 512 == 512
            vertex_list[10] = vertex_interp(iso, points[i+cube_size, j, k],
                                            points[i+cube_size, j, k+cube_size],
                                            val(i + cube_size, j, k),
                                            val(i + cube_size, j, k + cube_size))
        end
        if edge_entry & 1024 == 1024
            vertex_list[11] = vertex_interp(iso, points[i+cube_size, j+cube_size, k],
                                            points[i+cube_size, j+cube_size, k+cube_size],
                                            val(i + cube_size, j + cube_size, k),
                                            val(i + cube_size, j + cube_size,
                                                k + cube_size))
        end
        if edge_entry & 2048 == 2048
            vertex_list[12] = vertex_interp(iso, points[i, j+cube_size, k],
                                            points[i, j+cube_size, k+cube_size],
                                            val(i, j + cube_size, k),
                                            val(i, j + cube_size, k + cube_size))
        end

        table_entry = tri_table[cube_index]
        n = 1
        while table_entry[n] != -1
            vert = vertex_list[table_entry[n]+1]
            push!(vertices, vert)
            # if !ind_pushed
            push!(indices, (i, j, k))
            # ind_pushed = true
            n += 1
        end
        # ind_pushed = false
    end
    return vertices, indices
end

function get_edge_tri_table()
    edge_table = UInt16[0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905,
                        0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a,
                        0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a,
                        0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435,
                        0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
                        0x3a0, 0x2a9, 0x1a3, 0xaa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5,
                        0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460, 0x569, 0x663,
                        0x76a, 0x66, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a,
                        0x963, 0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5,
                        0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
                        0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c, 0xe5c, 0xf55,
                        0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3,
                        0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca,
                        0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf,
                        0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
                        0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55,
                        0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3,
                        0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0xff, 0x1f6, 0x6fa,
                        0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f,
                        0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460,
                        0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac,
                        0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, 0xd30, 0xc39, 0xf33,
                        0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636,
                        0x13a, 0x33, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f,
                        0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190,
                        0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c,
                        0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0]
    tri_table = Array{Array{Int8,1},1}([[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1,
                                         -1],
                                        [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1,
                                         -1],
                                        [3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1,
                                         -1],
                                        [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1,
                                         -1],
                                        [4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1,
                                         -1],
                                        [9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1,
                                         -1],
                                        [10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1,
                                         -1],
                                        [5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1,
                                         -1],
                                        [5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1,
                                         -1],
                                        [8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1,
                                         -1],
                                        [2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1,
                                         -1],
                                        [2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1,
                                         -1],
                                        [11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1,
                                         -1],
                                        [5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0,
                                         -1],
                                        [11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0,
                                         -1],
                                        [11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1,
                                         -1],
                                        [2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1,
                                         -1],
                                        [6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1,
                                         -1],
                                        [3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1,
                                         -1],
                                        [6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1,
                                         -1],
                                        [6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1,
                                         -1],
                                        [8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1,
                                         -1],
                                        [7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1],
                                        [3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1,
                                         -1],
                                        [0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1,
                                         -1],
                                        [9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6,
                                         -1],
                                        [8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1,
                                         -1],
                                        [5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11,
                                         -1],
                                        [0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1],
                                        [6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1,
                                         -1],
                                        [10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1,
                                         -1],
                                        [1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1,
                                         -1],
                                        [0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1,
                                         -1],
                                        [3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1,
                                         -1],
                                        [6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1,
                                         -1],
                                        [9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1,
                                         -1],
                                        [8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1],
                                        [3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1,
                                         -1],
                                        [10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1,
                                         -1],
                                        [10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1,
                                         -1],
                                        [2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1],
                                        [7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1,
                                         -1],
                                        [2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7,
                                         -1],
                                        [1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11,
                                         -1],
                                        [11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1,
                                         -1],
                                        [8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1],
                                        [0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1,
                                         -1],
                                        [7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1,
                                         -1],
                                        [7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1,
                                         -1],
                                        [10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1,
                                         -1],
                                        [0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1,
                                         -1],
                                        [7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1,
                                         -1],
                                        [6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1,
                                         -1],
                                        [4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1,
                                         -1],
                                        [10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3,
                                         -1],
                                        [8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1,
                                         -1],
                                        [1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1,
                                         -1],
                                        [10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1],
                                        [10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1,
                                         -1],
                                        [9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1,
                                         -1],
                                        [7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1,
                                         -1],
                                        [3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1],
                                        [7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1,
                                         -1],
                                        [3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1,
                                         -1],
                                        [6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1],
                                        [9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1,
                                         -1],
                                        [1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1],
                                        [4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10,
                                         -1],
                                        [7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1,
                                         -1],
                                        [6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1,
                                         -1],
                                        [0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1,
                                         -1],
                                        [6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1,
                                         -1],
                                        [0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10,
                                         -1],
                                        [11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5,
                                         -1],
                                        [6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1,
                                         -1],
                                        [5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1,
                                         -1],
                                        [9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1],
                                        [1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1],
                                        [10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1,
                                         -1],
                                        [0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1,
                                         -1],
                                        [11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1,
                                         -1],
                                        [9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1,
                                         -1],
                                        [7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1],
                                        [2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1,
                                         -1],
                                        [9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1,
                                         -1],
                                        [9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1],
                                        [1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1,
                                         -1],
                                        [0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1,
                                         -1],
                                        [10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4,
                                         -1],
                                        [2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1,
                                         -1],
                                        [0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11,
                                         -1],
                                        [0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1],
                                        [9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1,
                                         -1],
                                        [5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1],
                                        [5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1,
                                         -1],
                                        [8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1,
                                         -1],
                                        [9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1,
                                         -1],
                                        [1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1,
                                         -1],
                                        [3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4,
                                         -1],
                                        [4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1,
                                         -1],
                                        [9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3,
                                         -1],
                                        [11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1,
                                         -1],
                                        [2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1,
                                         -1],
                                        [9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1],
                                        [3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10,
                                         -1],
                                        [1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1,
                                         -1],
                                        [4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1,
                                         -1],
                                        [4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1,
                                         -1],
                                        [0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1,
                                         -1],
                                        [1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1],
                                        [0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1],
                                        [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                         -1, -1, -1]])
    return edge_table, tri_table
end


text: using Test, TestPackage0


text:
using SafeTestsets: @safetestset

# r"Main\..*\."

@safetestset "no type param" begin
    io = IOBuffer()
    struct Foo
        x::Int64
    end

    show(io, Foo(3))
    @show String(take!(io)) # because of safetestset it will be like "Main.##257.Foo(3)""
end

@safetestset "type param" begin
    using ShowTools: ShowTools
    io = IOBuffer()
    struct Foo{T}
        x::T
    end

    Base.show(io::Base.IO, x::Foo) = ShowTools.show_default_no_typeparams(io, x)
    show(io, Foo(3))
    @show String(take!(io)) # because of safetestset it will be like "Main.##123.Foo(3)""
    #Foo(3)
end

@safetestset "type param" begin
    using ShowTools: ShowTools
    io = IOBuffer()
    struct Foo{S, T}
        x::T
    end

    Base.show(io::Base.IO, x::Foo) = ShowTools.show_default_no_typeparams(io, x)
    show(io, Foo{Some{Float64}, Int64}(3))
    @show String(take!(io))
    # maybe Foo{Some{Float64}, _}(3)
end

@safetestset "type param" begin
    using ShowTools: ShowTools
    io = IOBuffer()
    struct Foo{S, T}
        x::T
    end

    # Base.show(io::Base.IO, x::Foo) = ShowTools.show_default_no_typeparams(io, x)
    show(io, Foo{:bitsbits, Int64}(3))
    @show String(take!(io))
    # maybe Foo{:bitsbits, _}(3)
end


text:
module LatticeBasisReduction

export gram_schmidt, LLL, diophantine, gauss, sda

include("misc.jl")
include("gram_schmidt.jl")
include("diophantine.jl")
include("gauss.jl")
include("sda.jl")
include("lll.jl")

end


text:
struct Config{
    #= Params =#
    MATMUL_SHAPE,               # MNK, overall shape of the MATMUL operation
    BLOCK_SHAPE,                # MNK, shape of each CTA tile
    WARPS_PER_BLOCK,            # scalar, number of warps per CTA

    MEM_A_WARP,                 # MK, shape of each warp tile during memory operations involving matrix A
    MEM_A_THREAD,               # MK, shape of each thread tile during memory operations involving matrix A

    MEM_B_WARP,                 # KN, shape of each warp tile during memory operations involving matrix B
    MEM_B_THREAD,               # KN, shape of each thread tile during memory operations involving matrix B

    MEM_CD_WARP,                # MN, shape of each warp tile during memory operations involving matrix C or D
    MEM_CD_THREAD,              # MN, shape of each thread tile during memory operations involving matrix C or D

    COMPUTE_WARP,               # MNK, shape of each warp tile during the inner loop computations
    COMPUTE_OP_SHAPE,           # MNK, shape of the operation used in the inner loop

    #= Layouts =#
    GLOBAL_A_LAYOUT,            # layout of the A matrix in global memory
    GLOBAL_B_LAYOUT,            # layout of the B matrix in global memory
    GLOBAL_C_LAYOUT,            # layout of the C matrix in global memory
    GLOBAL_D_LAYOUT,            # layout of the D matrix in global memory

    SHARED_A_LAYOUT,            # layout of the A matrix in shared memory
    SHARED_B_LAYOUT,            # layout of the B matrix in shared memory
    SHARED_C_LAYOUT,            # layout of the C matrix in shared memory
    SHARED_D_LAYOUT,            # layout of the D matrix in shared memory

    #= Operator =#
    OPERATOR,                   # which operator to use in the inner loop

    #= Is A & B stored in Column major order? This determines the iteration order of the parallellisation =#
    IS_A_COL_MAJOR,
    IS_B_COL_MAJOR
   }
end

@inline function Base.getproperty(conf::Type{Config{MATMUL_SHAPE, BLOCK_SHAPE, WARPS_PER_BLOCK, MEM_A_WARP, MEM_A_THREAD, MEM_B_WARP, MEM_B_THREAD, MEM_CD_WARP, MEM_CD_THREAD, COMPUTE_WARP, COMPUTE_OP_SHAPE, GLOBAL_A_LAYOUT, GLOBAL_B_LAYOUT, GLOBAL_C_LAYOUT, GLOBAL_D_LAYOUT, SHARED_A_LAYOUT, SHARED_B_LAYOUT, SHARED_C_LAYOUT, SHARED_D_LAYOUT, OPERATOR, IS_A_COL_MAJOR, IS_B_COL_MAJOR}}, sym::Symbol) where {MATMUL_SHAPE, BLOCK_SHAPE, WARPS_PER_BLOCK, MEM_A_WARP, MEM_A_THREAD, MEM_B_WARP, MEM_B_THREAD, MEM_CD_WARP, MEM_CD_THREAD, COMPUTE_WARP, COMPUTE_OP_SHAPE, GLOBAL_A_LAYOUT, GLOBAL_B_LAYOUT, GLOBAL_C_LAYOUT, GLOBAL_D_LAYOUT, SHARED_A_LAYOUT, SHARED_B_LAYOUT, SHARED_C_LAYOUT, SHARED_D_LAYOUT, OPERATOR, IS_A_COL_MAJOR, IS_B_COL_MAJOR}
    if sym == :launch_args
        (threads = WARPS_PER_BLOCK * 32,
         blocks = (MATMUL_SHAPE.M ÷ BLOCK_SHAPE.M, MATMUL_SHAPE.N ÷ BLOCK_SHAPE.N),
         shmem = 64 * 1024)

    # convenience accessors for typevars
    elseif sym == :matmul_shape
        MATMUL_SHAPE
    elseif sym == :block_shape
        BLOCK_SHAPE
    elseif sym == :warps_per_block
        WARPS_PER_BLOCK
    elseif sym == :mem_a_warp
        MEM_A_WARP
    elseif sym == :mem_a_thread
        MEM_A_THREAD
    elseif sym == :mem_b_warp
        MEM_B_WARP
    elseif sym == :mem_b_thread
        MEM_B_THREAD
    elseif sym == :mem_cd_warp
        MEM_CD_WARP
    elseif sym == :mem_cd_thread
        MEM_CD_THREAD
    elseif sym == :compute_warp
        COMPUTE_WARP
    elseif sym == :compute_op_shape
        COMPUTE_OP_SHAPE
    elseif sym == :global_a_layout
        GLOBAL_A_LAYOUT
    elseif sym == :global_b_layout
        GLOBAL_B_LAYOUT
    elseif sym == :global_c_layout
        GLOBAL_C_LAYOUT
    elseif sym == :global_d_layout
        GLOBAL_D_LAYOUT
    elseif sym == :shared_a_layout
        SHARED_A_LAYOUT
    elseif sym == :shared_b_layout
        SHARED_B_LAYOUT
    elseif sym == :shared_c_layout
        SHARED_C_LAYOUT
    elseif sym == :shared_d_layout
        SHARED_D_LAYOUT
    elseif sym == :operator
        OPERATOR
    elseif sym == :is_a_col_major
        IS_A_COL_MAJOR
    elseif sym == :is_b_col_major
        IS_B_COL_MAJOR

    # fallback
    else
        getfield(conf, sym)
    end
end

function heuristic_block_shape(shared_a_layout, shared_b_layout, shared_c_layout, shared_d_layout)
    # Determining the tile size of each block is a little trickier.
    # We apply the following heuristics:
    # 1) Ideally, the block shape in the M and N dimensions should be square or
    #    nearly-square to maximise data reuse. More specifically, tiling sizes
    #    of the form (M = k, N = k) and (M = 2k, N = k) work well.
    # 2) The tile size should be as large as possible.
    # 3) The size in the M and N dimension is limited by the fact that a tile of
    #    the C (and D) matrix of that size must fit in shared memory.
    # 4) The size in the K dimension is limited by the fact that both a M x K tile
    #    of A and a K x N tile of B must fit in shared memory, at the same time.

    num_bytes_A(M, N, K) = prod(Layout.physical_size(shared_a_layout, (M = M, K = K))) * sizeof(Layout.eltype(shared_a_layout))
    num_bytes_B(M, N, K) = prod(Layout.physical_size(shared_b_layout, (K = K, N = N))) * sizeof(Layout.eltype(shared_b_layout))
    num_bytes_C(M, N, K) = prod(Layout.physical_size(shared_c_layout, (M = M, N = N))) * sizeof(Layout.eltype(shared_c_layout))
    num_bytes_D(M, N, K) = prod(Layout.physical_size(shared_d_layout, (M = M, N = N))) * sizeof(Layout.eltype(shared_d_layout))

    next_MN(M, N, K) = M == N ? (2 * M, N, K) : (M, 2 * N, K)
    next_K(M, N, K) = (M, N, 2 * K)

    cur = 1, 1, 1     # M, N, K
    nxt = next_MN(cur...)

    while (max(num_bytes_C(nxt...), num_bytes_D(nxt...)) <= 64 * 1024)
        cur = nxt
        nxt = next_MN(cur...)
    end

    nxt = next_K(cur...)

    while (num_bytes_A(nxt...) + num_bytes_B(nxt...) <= 64 * 1024)
        cur = nxt
        nxt = next_K(cur...)
    end

    return (M = cur[1], N = cur[2], K = cur[3])
end

# Helper function that returns the logical size of a set of adjacent elements, taking care not
# to make the size larger than the parent tile
function adjacent_elements(num, parent_size, is_col_major)
    p = Tuple(parent_size)

    if is_col_major
        t = (min(num, p[1]), num ÷ min(num, p[1]))
    else
        t = (num ÷ min(num, p[2]), min(num, p[2]))
    end

    return typeof(parent_size)(t)
end

function get_config(; gemm_shape, operator, global_a_layout, global_c_layout, kwargs...)
    params = Dict(kwargs)

    # Use some simple heuristics to get sensible defaults for parameters the user does not specify.

    # Get the global layouts for B & D.
    # Fallback to the layouts of A & C, respectively.
    global_b_layout = get(params, :global_b_layout, global_a_layout)
    global_d_layout = get(params, :global_d_layout, global_c_layout)

    # Get the shared layouts for A, B, C, D.
    # For A & B, add padding to reduce bank conflicts, but preserve 128-bit (16 byte) alignment.
    shared_a_layout = get(params, :shared_a_layout,
                          Layout.Padded{global_a_layout, 16 ÷ sizeof(Layout.eltype(global_a_layout))})
    shared_b_layout = get(params, :shared_b_layout,
                          Layout.Padded{global_b_layout, 16 ÷ sizeof(Layout.eltype(global_b_layout))})
    shared_c_layout = get(params, :shared_c_layout, global_c_layout)
    shared_d_layout = get(params, :shared_d_layout, global_d_layout)

    # Apply heuristic for block shape
    block_shape = get(params, :block_shape,
        heuristic_block_shape(shared_a_layout, shared_b_layout, shared_c_layout, shared_d_layout))

    # 8 warps in a 4 x 2 arrangement usually works well
    warps_per_block = get(params, :warps_per_block, 8)
    op_shape = Operator.shape(operator)
    compute_warp = get(params, :compute_warp,
                       (M = block_shape.M ÷ 4, N = block_shape.N ÷ 2, K = op_shape.K))

    # Is the layout col-major or not? This is needed to find good values for mem_a_warp, mem_b_warp, etc.
    # TODO: Let the layouts handle this?
    is_a_col_major = get(params, :is_a_col_major, true)
    is_b_col_major = get(params, :is_b_col_major, true)
    is_cd_col_major = get(params, :is_cd_col_major, true)

    # Heuristics for memory tiling sizes:
    # 1) The tiles should encompass 128 bits (16 bytes) to enable vectorisation.
    # 2) The tiles should be as small as possible (i.e. each thread exactly 128 bits) to enable coalescing.

    mem_a_warp = get(params, :mem_a_warp,
        adjacent_elements(32 * 16 ÷ sizeof(Layout.eltype(global_a_layout)), (M = block_shape.M, K = block_shape.K), is_a_col_major))
    mem_b_warp = get(params, :mem_b_warp,
        adjacent_elements(32 * 16 ÷ sizeof(Layout.eltype(global_b_layout)), (K = block_shape.K, N = block_shape.N), is_b_col_major))
    mem_cd_warp = get(params, :mem_cd_warp,
        adjacent_elements(32 * 16 ÷ sizeof(Layout.eltype(global_c_layout)), (M = block_shape.M, N = block_shape.N), is_cd_col_major))

    mem_a_thread = get(params, :mem_a_thread,
        adjacent_elements(16 ÷ sizeof(Layout.eltype(global_a_layout)), (M = block_shape.M, K = block_shape.K), is_a_col_major))
    mem_b_thread = get(params, :mem_b_thread,
        adjacent_elements(16 ÷ sizeof(Layout.eltype(global_b_layout)), (K = block_shape.K, N = block_shape.N), is_b_col_major))
    mem_cd_thread = get(params, :mem_cd_thread,
        adjacent_elements(16 ÷ sizeof(Layout.eltype(global_c_layout)), (M = block_shape.M, N = block_shape.N), is_cd_col_major))

    return Config{
        #= Params =#
        gemm_shape,
        block_shape,
        warps_per_block,
        mem_a_warp,
        mem_a_thread,
        mem_b_warp,
        mem_b_thread,
        mem_cd_warp,
        mem_cd_thread,
        compute_warp,
        op_shape,

        #= Layouts =#
        global_a_layout,
        global_b_layout,
        global_c_layout,
        global_d_layout,

        shared_a_layout,
        shared_b_layout,
        shared_c_layout,
        shared_d_layout,

        #= Operators =#
        operator,

        #= Is A & B Col Major? =#
        is_a_col_major,
        is_b_col_major,
    }
end


text:
using Printf

function distcheck(f::Function, rep::Int=10000, Δ::Int=3)
    smpl = f(rep)
    counts = Dict(k => count(smpl .== k) for k in unique(smpl))
    expected = rep / length(counts)
    lbound = expected * (1 - 0.01Δ)
    ubound = expected * (1 + 0.01Δ)
    noobs = count(x -> !(lbound ≤ x ≤ ubound), values(counts))
    if noobs > 0 warn(@sprintf "%2.4f%% values out of bounds" noobs / rep) end
    return counts
end

# Dice5 check
distcheck(x -> rand(1:5, x))
# Dice7 check
distcheck(dice7)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    LineBufferedIO(dest::IO)

A thread safe line buffered IO wrapper which buffers writes until a full line
(delimited by '\\n') is received. Full lines are written to the downstream
`dest` `IO`.

To ensure that the tail of the stream is written (even without a trailing
'\\n'), be sure to call `close()`.
"""
struct LineBufferedIO{DestIO<:IO} <: IO
    buf::IOBuffer
    buf_lock::ReentrantLock
    dest::DestIO
end

function LineBufferedIO(dest::IO)
    LineBufferedIO(IOBuffer(), ReentrantLock(), dest)
end

function Base.unsafe_write(io::LineBufferedIO, p::Ptr{UInt8}, n::UInt)
    first_newline = nothing
    last_newline = nothing
    for i = 1:n
        if unsafe_load(p, i) == UInt8('\n')
            first_newline = i
            break
        end
    end
    if !isnothing(first_newline)
        for i = n:-1:1
            if unsafe_load(p, i) == UInt8('\n')
                last_newline = i
                break
            end
        end
    end
    firstline = nothing
    unbuf_start = nothing
    unbuf_end = 0
    lock(io.buf_lock) do
        # Invariant: `buf` buffers the prefix of the current line
        buf = io.buf
        if isnothing(first_newline)
            # No newline: just buffer the whole thing
            unsafe_write(buf, p, n)
            unbuf_start = 1
            unbuf_end = 0
        else
            if position(buf) != 0
                unsafe_write(buf, p, first_newline)
                # At this point we have one full line in the buffer
                firstline = take!(buf)
                unbuf_start = first_newline+1
            else
                unbuf_start = 1
            end
            if last_newline != n
                # Also record any trailing chars into the buffer
                unsafe_write(buf, p+last_newline, n-last_newline)
            end
            unbuf_end = last_newline
        end
    end
    if !isnothing(firstline)
        write(io.dest, firstline)
    end
    linestart = unbuf_start
    for i = unbuf_start:unbuf_end
        if unsafe_load(p, i) == UInt8('\n')
            unsafe_write(io.dest, p + linestart-1, i - linestart + 1)
            linestart = i+1
        end
    end
    return n
end

function Base.write(io::LineBufferedIO, c::UInt8)
    Base.unsafe_write(io, Ref(c), 1)
end

function Base.close(io::LineBufferedIO)
    local remaining_line
    lock(io.buf_lock) do
        remaining_line = take!(io.buf)
        close(io.buf)
    end
    if !isempty(remaining_line)
        write(io.dest, remaining_line)
    end
end

Base.isopen(io::LineBufferedIO) = lock(()->isopen(io.buf), io.buf_lock)