text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants

using LaTeXStrings
using Plots
# gr()
# plotlyjs()
pyplot()

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

nr=128
rmax = 500 * AU
nz=128
zmax = 200 * AU

# Get the interpolator object for this particular setup.
interp_S = DiskTracer.model.get_S_nu_grid(pars, mol, mol.nu_0, nr, rmax, nz, zmax)
interp_U = DiskTracer.model.get_Upsilon_nu_grid(pars, mol, mol.nu_0, nr, rmax, nz, zmax)
interp_D2 = DiskTracer.model.get_DeltaV2_grid(pars, mol, nr, rmax, nz, zmax)

# Make a 2D plot of the source function and Upsilon
# density structure
function plot_S(pars::AbstractParameters)

    nr = 128
    nz = 128
    rs = linspace(0.1, 600 * AU, nr)
    zs = linspace(0.1, 300 * AU, nz)

    rr = rs./AU
    zz = zs./AU

    Ss = Array{Float64}(nz, nr)

    for i=1:nz
        for j=1:nr
            Ss[i,j] = interp_S(rs[j], zs[i])
        end
    end

    p = heatmap(Ss, xlabel=L"$r$ [AU]", ylabel=L"$z$ [AU]") #, size=(800,800))
    savefig(p, "source_function_interp.png")

end

function plot_Upsilon(pars::AbstractParameters)

    nr = 128
    nz = 128
    rs = linspace(0.1, 600 * AU, nr)
    zs = linspace(0.1, 300 * AU, nz)

    Us = Array{Float64}(nz, nr)

    for i=1:nz
        for j=1:nr
            Us[i,j] = interp_U(rs[j], zs[i])
        end
    end

    p = heatmap(Us, xlabel=L"$r$ [AU]", ylabel=L"$z$ [AU]") #, size=(800,800))
    savefig(p, "upsilon_interp.png")

end

function plot_DeltaV2(pars::AbstractParameters)

    nr = 128
    nz = 128
    rs = linspace(0.1, 600 * AU, nr)
    zs = linspace(0.1, 300 * AU, nz)

    D2s = Array{Float64}(nz, nr)

    for i=1:nz
        for j=1:nr
            D2s[i,j] = interp_D2(rs[j], zs[i])
        end
    end

    p = heatmap(D2s, xlabel=L"$r$ [AU]", ylabel=L"$z$ [AU]") #, size=(800,800))
    savefig(p, "DeltaV2_interp.png")

end

plot_S(pars)
plot_Upsilon(pars)
plot_DeltaV2(pars)


text:
using StructLH
using Random, Test

# Construct ROV object with child objects
struct RTP2
    r1
    r2
end

struct RTP3
    r4
    r2
end

struct RTP1
    x :: RTP2
    y :: RTP3
end

mutable struct ROV
    x
    y
    z
end

function string_sum(sV :: Vector{String})
    outStr = "";
    for s in sV
        outStr = outStr * s;
    end
    return outStr
end

# Because `hasmethod(sum, Tuple{Vector{Any}})== true`, but `sum(["a", "b"])` errors, we cannot use sum for the test
function float_sum(sV :: Vector{F1}) where F1 <: Real
    return sum(sV);
end

function retrieve_test()
    @testset "Retrieve property" begin
        x = RTP2(1, 2);
        y = RTP3(4, 5);
        z = RTP1(x, y);

        @test has_property(z, :r1)
        @test !has_property(z, :r12)
        r1 = retrieve_property(z, :r1);
        @test r1 == x.r1
        r2 = retrieve_property(z, :r2);
        @test r2 == x.r2
        r3 = retrieve_property(z, :r3);
        @test isnothing(r3)
        r4 = retrieve_property(z, :x);
        @test isa(r4, RTP2)

        r2a = retrieve_child_property(z, :x, :r2);
        @test r2a == x.r2
        r2b = retrieve_child_property(z, :y, :r2);
        @test r2b == y.r2
        # Try an object that is in a different child
        @test isnothing(retrieve_child_property(z, :x, :r4))
    end
end

function reduce_ov_test()
    @testset "Reduce object vector" begin
        rov1 = ROV(rand(3,2), 1.23, "a");
        rov2 = ROV(rand(3,2), 3.2, "b");
        rov3 = ROV(rand(3,2), 0.6, "c");
        objV = [rov1, rov2, rov3];
        oOut = reduce_object_vector(objV, float_sum);

        @test oOut.y ≈ rov1.y .+ rov2.y .+ rov3.y
        @test oOut.x ≈ rov1.x .+ rov2.x .+ rov3.x

        oOut2 = reduce_object_vector(objV, string_sum, fieldTypes = [String]);
        @test oOut2.z == "abc"
    end
end

mutable struct STD1
    x
end

mutable struct STD2
    x
    y
end

Base.isequal(s1 :: STD1, s2 :: STD1) = isequal(s1.x, s2.x);
Base.isequal(s1 :: STD2, s2 :: STD2) =
    isequal(s1.x, s2.x)  &&  isequal(s1.y, s2.y);

function struct2dict_test()
    @test isequal(struct2dict(1.2), NodeInfo(1.2))
    @test isequal(struct2dict([1, 2]), NodeInfo{Array{Int64, 1}}((2,), [1, 2]))

    @test isequal(struct2dict("abc"), NodeInfo("abc"))
    @test isequal(struct2dict(["abc", "def"]),
        NodeInfo{Array{String,1}}((2,), ["abc", "def"]))

    @test isequal(struct2dict(STD1([1, 2])),
        Dict{Symbol, Any}([:x => NodeInfo{Array{Int64,1}}((2,), [1, 2])]))

    s1 = STD2(1.2, STD1([1,2]));
    d1 = struct2dict(s1);
    @test isequal(d1[:x], NodeInfo(1.2))
    @test isequal(d1[:y], struct2dict(STD1([1,2])))

    s11 = STD2(0.0, STD1([1]));
    dict2struct!(s11, d1);
    @test isequal(s11.x, s1.x)
    @test isequal(s11, s1)
    d11 = struct2dict(s11);
    @test isequal(d1, d11)
end

@testset "All" begin
    retrieve_test()
    reduce_ov_test()
    struct2dict_test()
    include("merge_test.jl");
    include("apply_fct_test.jl");
    include("reduce_test.jl");
    include("helpers_test.jl");
end

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

ext: jl


text:
module EchoviewEcs

export load

MYFLOAT = Float32

"""
    FLOAT_LABELS

List of known calibration fields that contain floating point
values. Labels NOT in this list are assumed to be `String`s.

"""
FLOAT_LABELS = ["AbsorptionCoefficient",
     "EK60SaCorrection",
     "Ek60TransducerGain",
     "Frequency",
     "MajorAxis3dbBeamAngle",
     "MajorAxisAngleOffset",
     "MajorAxisAngleSensitivity",
     "MinorAxis3dbBeamAngle",
     "MinorAxisAngleOffset",
     "MinorAxisAngleSensitivity",
     "SoundSpeed",
     "TransmittedPower",
     "TransmittedPulseLength",
     "TwoWayBeamAngle"]

function load(lines::Vector{String})
    fileset = Dict()
    n = []
    m = fileset
    for line in lines
        if occursin("#", line)
            line = split(line, "#")[1]
        end
        if startswith(strip(line),"SourceCal")
            m = copy(fileset)
            push!(n,m)
        end
        if occursin("=", line) && m != nothing
            d = split(line, "=")
            k = strip(d[1])
            v =strip(d[2])
            if k in FLOAT_LABELS
                v = Base.parse(MYFLOAT,v)
            end
            m[k] =v
        end
    end

    return n
end

load(f::IOStream) = load(readlines(f))

"""
    load(filename::AbstractString)

Loads an Echoview calibration supplement file, returning
a `Vector` of `Dict`. Each `Dict` represents a transducer calibration
with keys and values containing the various calibration settings.

See [Echoview calibration supplement files](http://support.echoview.com/WebHelp/Reference/File_formats/Echoview_calibration_supplement_files.html).

"""
function load(filename::AbstractString)
    open(filename) do f
        load(f)
    end
end

end # module


text:
using ImageFiltering, Statistics, Test
using ImageFiltering: IdentityUnitRange

@testset "mapwindow" begin
    function groundtruth(f, A, window::Tuple)
        Aex = copy(A)
        hshift = map(x->x>>1+1, window)
        for Ishift in CartesianIndices(window)
            for I in CartesianIndices(size(A))
                Aex[I] = f(Aex[I], A[map(d->clamp(I[d]+Ishift[d]-hshift[d], 1, size(A,d)), 1:ndims(A))...])
            end
        end
        Aex
    end
    groundtruth(f, A, window) = groundtruth(f, A, (window...,))
    # 1d case
    A = [0.1,0.3,0.2,0.3,0.4]
    mm = mapwindow(extrema, A, 1)
    @test mm == [(a,a) for a in A]
    mm = mapwindow(extrema, A, 2)
    @test mm == [(0.1,0.1),(0.1,0.3),(0.2,0.3),(0.2,0.3),(0.3,0.4)]
    mm = mapwindow(extrema, A, 3)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.2,0.3),(0.2,0.4),(0.3,0.4)]
    mm = mapwindow(extrema, A, 4)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.1,0.3),(0.2,0.4),(0.2,0.4)]
    mm = mapwindow(extrema, A, 5)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.1,0.4),(0.2,0.4),(0.2,0.4)]
    mm = mapwindow(extrema, A, 6)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.1,0.4),(0.1,0.4),(0.2,0.4)]
    mm = mapwindow(extrema, A, 7)
    @test mm == [(0.1,0.3),(0.1,0.4),(0.1,0.4),(0.1,0.4),(0.2,0.4)]
    A = [0.1,0.3,0.5,0.4,0.2]
    mm = mapwindow(extrema, A, 1)
    @test mm == [(a,a) for a in A]
    mm = mapwindow(extrema, A, 2)
    @test mm == [(0.1,0.1),(0.1,0.3),(0.3,0.5),(0.4,0.5),(0.2,0.4)]
    mm = mapwindow(extrema, A, 3)
    @test mm == [(0.1,0.3),(0.1,0.5),(0.3,0.5),(0.2,0.5),(0.2,0.4)]
    mm = mapwindow(extrema, A, 4)
    @test mm == [(0.1,0.3),(0.1,0.5),(0.1,0.5),(0.2,0.5),(0.2,0.5)]
    mm = mapwindow(extrema, A, 5)
    @test mm == [(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.2,0.5),(0.2,0.5)]
    mm = mapwindow(extrema, A, 6)
    @test mm == [(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.2,0.5)]
    mm = mapwindow(extrema, A, 7)
    @test mm == [(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.2,0.5)]
    # 2d case
    A = rand(5,5)/10
    A[2,2] = 0.8
    A[4,4] = 0.6
    for w in ((2,2), (2,3), (3,2), (3,3), (2,5))
        mm = mapwindow(extrema, A, w)
        maxval = last.(mm)
        Amax = groundtruth(max, A, w)
        @test maxval == Amax
        minval = first.(mm)
        Amin = groundtruth(min, A, w)
        @test minval == Amin
    end
    # 3d case
    A = rand(5,5,5)/10
    A[2,2,2] = 0.7
    A[4,4,2] = 0.4
    A[2,2,4] = 0.5
    for w in ((2,2,2), (2,3,2), (3,2,2), (2,2,3), (3,3,3), (2,5,3))
        mm = mapwindow(extrema, A, w)
        maxval = last.(mm)
        Amax = groundtruth(max, A, w)
        @test maxval == Amax
        minval = first.(mm)
        Amin = groundtruth(min, A, w)
        @test minval == Amin
    end

    # median
    for f in (median, median!)
        a = [1,1,1,2,2,2]
        @test mapwindow(f, a, -1:1) == a
        @test mapwindow(f, a, (-1:1,)) == a
        @test mapwindow(f, a, -2:2) == a
        @test mapwindow(f, a, -3:3) == a
        b = [1,100,1,2,-1000,2]
        @test mapwindow(f, b, -1:1) == [1,1,2,1,2,2]
        @test mapwindow(f, b, -2:2) == a

        A = [1 5 -2 3 7;
             2 0 3  4 4;
             3 3 6  2 5;
             1 -3 5 3 0]
        @test mapwindow(f, A, (3,3)) == [1 1 3 3 4;
                                         2 3 3 4 4;
                                         2 3 3 4 4;
                                         1 3 3 3 2]
    end

    # resolve_f
    replace_me(x) = 1
    ImageFiltering.MapWindow.replace_function(::typeof(replace_me)) = x -> 2
    @test mapwindow!(replace_me, randn(3), randn(3), (1,)) == [2,2,2]

    function groundtruth(f, A, window::Tuple; indices=nothing)
        mapwindow(f,A,window)[indices...]
    end
    for (f,img, window, imginds) ∈ [
            (mean   , randn(10)  , (1,)      , (1:2:10,)),
            (median!, randn(10)  , (-1:1,)   , (1:2:8,)),
            (mean   , randn(10)  , (-1:1,)   , (1:2:8,)),
            (mean   , randn(10,5), (-1:1,0:0), (1:2:8,1:3)),
            (mean   , randn(10,5), (-1:1,0:0), (Base.OneTo(2),1:3)),
        ]

        expected = groundtruth(f,img,window,indices=imginds)
        @test expected == @inferred mapwindow(f,img,window,indices=imginds)
        out = similar(expected)
        @test expected == @inferred mapwindow!(f,out,img,window,indices=imginds)
    end
    for (inds, kw) ∈ [((Base.OneTo(3),), (border="replicate", indices=2:2:7)),
                        ((2:7,), (indices=2:7,)),
                        ((Base.OneTo(10),), ())
                       ]
        @test inds == axes(mapwindow(mean, randn(10), (3,); kw...))
    end

    img_48 = 10*collect(1:10)
    @test mapwindow(first, img_48, (1,), border=Inner()) == img_48
    res_48 = mapwindow(first, img_48, (0:1,), border=Inner())
    @test axes(res_48) === (IdentityUnitRange(1:9),)
    @test res_48 == img_48[axes(res_48)...]
    inds_48 = 2:2:8
    @test mapwindow(first, img_48, (0:2,),
                    border=Inner(),
                    indices=inds_48) == img_48[inds_48]

    @testset "desugaring window argument #58" begin
        img58 = rand(10)
        canonical_window = (-1:1,)
        truth = mapwindow(median, img58, canonical_window)
        for window in [3, (3,), [3], [-1:1], -1:1]
            @test truth == mapwindow(median, img58, window)
            out = similar(img58)
            @test truth == mapwindow!(median,out, img58, window)
        end
    end

    @testset ">3D mapwindow #105" begin
        img105 = ones(5,5,5,5,5)
        out105 = mapwindow(sum, img105, [1,1,1,1,3]; border = Fill(0))
        foo = centered([1])
        bar = centered([1,1,1])
        @test out105 == imfilter(img105, kernelfactors((foo, foo, foo, foo, bar)), Fill(0))
        out1052 = mapwindow(sum, img105, [1,1,1,1,3])
        @test all(isequal(3.0), out1052)
    end
end


text:
module PINNS

using ArgParse
using BioSequences
using RCall
using BioBridgeR.APE
using GeneticVariation: NG86, pairwise_dNdS

export cmd

function parse_command_line()
    s = ArgParseSettings()
    s.add_help = true

    @add_arg_table s begin
        "nullsim"
            help = "Simulate null distribution."
            action = :command
        "process"
            help = "Calculate statistics for real data and null distributions."
            action = :command
    end

    @add_arg_table s["nullsim"] begin
        "inputfile"
            help = "An input fasta file."
            arg_type = String
            required = true
        "reps"
            help = "Number of simulations to produce."
            arg_type = Int64
            required = true
        "seed"
            help = "A seed to set for RNG"
            arg_type = Int64
            default = rand(1:1000000)
        "--treemethod", "-t"
            help = "Method used to construct the starting tree - 'nj' or 'upgma'."
            arg_type = String
            default = "nj"
        "--submodel", "-s"
            help = "The substitution model used in model fitting and ancestral sequence estimation - 'JC' or 'F81'."
            arg_type = String
            default = "JC"
        "--ancmethod", "-a"
            help = "the criterion to assign the internal nodes during ancestral sequence reconstruction. - 'marginal', 'ml', or 'bayes'."
            arg_type = String
            default = "ml"
    end

    @add_arg_table s["process"] begin
        "realfile"
            help = "The FASTA file containing your actual sequence alignment."
            arg_type = String
        "simfile"
            help = "The FASTA file containing the similated sequence alignments."
            arg_type = String
    end

    return parse_args(s)
end

function cmd()
    arguments = parse_command_line()
    if arguments["%COMMAND%"] == "nullsim"
        args = arguments["nullsim"]
        nullsim(args["inputfile"], args["seed"], args["reps"], args["submodel"], args["treemethod"], args["ancmethod"])
    end
    if arguments["%COMMAND%"] == "process"
        args = arguments["process"]
        process(args["realfile"], args["simfile"])
    end
end

include("simulate.jl")
include("process.jl")

end # Module PINNS.


text:
@testset "cache per configuration" begin
    # cache key should be same for the same configurations
    let
        k1 = JET.get_cache_key(JETAnalyzer())

        k2 = JET.get_cache_key(JETAnalyzer())

        @test k1 == k2
    end

    # cache key should be different for different configurations
    let
        analyzer1 = JETAnalyzer(; max_methods=3)
        k1 = JET.get_cache_key(analyzer1)

        analyzer2 = JETAnalyzer(; max_methods=4)
        k2 = JET.get_cache_key(analyzer2)

        @test k1 ≠ k2
    end

    # configurations other than `InferenceParams` and `ReportPass`
    # shouldn't affect the cache key identity
    let
        analyzer1 = JETAnalyzer(; toplevel_logger=nothing)
        k1 = JET.get_cache_key(analyzer1)

        analyzer2 = JETAnalyzer(; toplevel_logger=IOBuffer())
        k2 = JET.get_cache_key(analyzer2)

        @test k1 == k2
    end

    # cache key should be different for different report passes
    let
        analyzer1 = JETAnalyzer(; report_pass=JET.BasicPass())
        k1 = JET.get_cache_key(analyzer1)

        analyzer2 = JETAnalyzer(; report_pass=JET.SoundPass())
        k2 = JET.get_cache_key(analyzer2)

        @test k1 ≠ k2
    end

    # end to end test
    let
        m = Module()
        @eval m begin
            foo(a::Val{1}) = 1
            foo(a::Val{2}) = 2
            foo(a::Val{3}) = 3
            foo(a::Val{4}) = undefvar
        end

        # run first analysis and cache
        result = @eval m $report_call((Int,); max_methods=3) do a
            foo(Val(a))
        end
        @test isempty(get_reports(result))

        # should use the cached result
        result = @eval m $report_call((Int,); max_methods=3) do a
            foo(Val(a))
        end
        @test isempty(get_reports(result))

        # should re-run analysis, and should get a report
        result = @eval m $report_call((Int,); max_methods=4) do a
            foo(Val(a))
        end
        @test any(get_reports(result)) do r
            isa(r, GlobalUndefVarErrorReport) &&
            r.name === :undefvar
        end

        # should run the cached previous result
        result = @eval m $report_call((Int,); max_methods=4) do a
            foo(Val(a))
        end
        @test any(get_reports(result)) do r
            isa(r, GlobalUndefVarErrorReport) &&
            r.name === :undefvar
        end
    end
end

@testset "report no method matching" begin
    # if there is no method matching case, it should be reported
    let
        # NOTE: we can't just wrap them into `let`, closures can't be inferred correctly
        m = Module()
        result = Core.eval(m, quote
            foo(a::Integer) = :Integer
            $report_call((AbstractString,)) do a
                foo(a)
            end
        end)
        @test length(get_reports(result)) === 1
        report = first(get_reports(result))
        @test report isa NoMethodErrorReport
        @test report.t === Tuple{typeof(m.foo), AbstractString}
    end

    let
        result = report_call(()->sum([]))
        @test length(get_reports(result)) === 1
        report = first(get_reports(result))
        @test report isa SeriousExceptionReport
        @test report.err isa MethodError
        @test report.err.f === zero
    end

    # if there is no method matching case in union-split, it should be reported
    let
        m = Module()
        result = Core.eval(m, quote
            foo(a::Integer) = :Integer
            foo(a::AbstractString) = "AbstractString"

            $report_call(a->foo(a), (Union{Nothing,Int},))
        end)

        @test length(get_reports(result)) === 1
        report = first(get_reports(result))
        @test report isa NoMethodErrorReport
        @test report.t == Any[Tuple{typeof(m.foo), Nothing}]
    end
end

@testset "report local undefined variables" begin
    let
        result = report_call((Bool,)) do b
            if b
                bar = rand(Int)
                return bar
            end
            return bar # undefined in this pass
        end
        @test length(get_reports(result)) === 1
        r = first(get_reports(result))
        @test r isa LocalUndefVarErrorReport
        @test r.name === :bar
        @test last(r.vst).line == (@__LINE__)-6
    end

    # deeper level
    let
        m = @fixturedef begin
            function foo(b)
                if b
                    bar = rand(Int)
                    return bar
                end
                return bar # undefined in this pass
            end
            baz(a) = foo(a)
        end

        result = Core.eval(m, :($report_call(baz, (Bool,))))
        @test length(get_reports(result)) === 1
        r = first(get_reports(result))
        @test r isa LocalUndefVarErrorReport
        @test r.name === :bar
        @test last(r.vst).line == (@__LINE__)-10

        # works when cached
        result = Core.eval(m, :($report_call(baz, (Bool,))))
        @test length(get_reports(result)) === 1
        r = first(get_reports(result))
        @test r isa LocalUndefVarErrorReport
        @test r.name === :bar
        @test last(r.vst).line == (@__LINE__)-18
    end

    # try to exclude false negatives as possible (by collecting reports in after-optimization pass)
    let
        result = report_call((Bool,)) do b
            if b
                bar = rand()
            end

            return if b
                return bar # this shouldn't be reported
            else
                return nothing
            end
        end
        @test isempty(get_reports(result))
    end

    let
        result = report_call((Bool,)) do b
            if b
                bar = rand()
            end

            return if b
                return nothing
            else
                # ideally we want to have report for this pass, but tons of work will be
                # needed to report this pass
                return bar
            end
        end
        @test_broken length(get_reports(result)) === 1 &&
            first(get_reports(result)) isa LocalUndefVarErrorReport &&
            first(get_reports(result)).name === :bar
    end

    let
        result = report_call((Int,)) do a
            function inner(n)
                if n > 0
                   a = n
                end
            end
            inner(rand(Int))
            return a
        end
        @test isempty(get_reports(result))
    end

    let # should work for top-level analysis
        res = @analyze_toplevel begin
            foo = let
                if rand(Bool)
                    bar = rand(Int)
                else
                    bar # undefined in this pass
                end
            end
        end
        @test length(res.inference_error_reports) === 1 &&
              first(res.inference_error_reports) isa LocalUndefVarErrorReport &&
              first(res.inference_error_reports).name === :bar
    end
end

@testset "report undefined (global) variables" begin
    let
        result = report_call(()->foo)
        @test length(get_reports(result)) === 1
        @test first(get_reports(result)) isa GlobalUndefVarErrorReport
        @test first(get_reports(result)).name === :foo
    end

    # deeper level
    let
        m = @fixturedef begin
            foo(bar) = bar + baz
            qux(a) = foo(a)
        end

        result = Core.eval(m, :($report_call(qux, (Int,))))
        @test length(get_reports(result)) === 1
        @test first(get_reports(result)) isa GlobalUndefVarErrorReport
        @test first(get_reports(result)).name === :baz

        # works when cached
        result = Core.eval(m, :($report_call(qux, (Int,))))
        @test length(get_reports(result)) === 1
        @test first(get_reports(result)) isa GlobalUndefVarErrorReport
        @test first(get_reports(result)).name === :baz
    end
end

@testset "report non-boolean condition error" begin
    # simple case
    let
        result = report_call((Int,)) do a
            a ? a : nothing
        end
        @test length(get_reports(result)) === 1
        er = first(get_reports(result))
        @test er isa NonBooleanCondErrorReport
        @test er.t === Int
    end

    # don't report when a type can be `Bool` (Bool ⊑ type)
    let
        result = report_call((Integer,)) do a
            a ? a : nothing
        end
        @test isempty(get_reports(result))
    end

    # report union split case
    let
        result = report_call((Union{Nothing,Bool},)) do a
            a ? a : false
        end
        @test length(get_reports(result)) === 1
        let r = first(get_reports(result))
            @test r isa NonBooleanCondErrorReport
            @test r.t == [Nothing]
            @test occursin("for 1 of union split cases", r.msg)
        end
    end

    let
        result = report_call() do
            anyary = Any[1,2,3]
            first(anyary) ? first(anyary) : nothing
        end
        @test isempty(get_reports(result)) # very untyped, we can't report on this ...
    end
end

@testset "UndefKeywordError" begin
    isa2(t) = x -> isa(x, t)
    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            foo(a; #= can be undef =# kw) =  a, kw
            $report_call(foo, (Any,))
        end)
        @test !isempty(get_reports(result))
        @test any(get_reports(result)) do r
            r isa SeriousExceptionReport || return false
            err = r.err
            err isa UndefKeywordError && err.var === :kw
        end
        # there shouldn't be duplicated report for the `throw` call
        @test !any(isa2(UncaughtExceptionReport), get_reports(result))
    end
end

@testset "DivideError" begin
    isa2(t) = x -> isa(x, t)
    let
        apply(f, args...) = f(args...)

        result = report_call() do
            apply(div, 1, 0)
        end
        @test !isempty(get_reports(result))
        @test any(isa2(DivideErrorReport), get_reports(result))

        result = report_call() do
            apply(rem, 1, 0)
        end
        @test !isempty(get_reports(result))
        @test any(isa2(DivideErrorReport), get_reports(result))

        # JET analysis isn't sound
        result = report_call((Int,Int)) do a, b
            apply(div, a, b)
        end
        @test isempty(get_reports(result))
    end
end

@testset "report `throw` calls" begin
    # simplest case
    let
        result = report_call(()->throw("foo"))
        @test !isempty(get_reports(result))
        @test first(get_reports(result)) isa UncaughtExceptionReport
    end

    # throws in deep level
    let
        foo(a) = throw(a)
        result = report_call(()->foo("foo"))
        @test !isempty(get_reports(result))
        @test first(get_reports(result)) isa UncaughtExceptionReport
    end

    # don't report possibly false negative `throw`s
    let
        foo(a) = a ≤ 0 ? throw("a is $(a)") : a
        result = report_call(foo, (Int,))
        @test isempty(get_reports(result))
    end

    # constant prop sometimes helps exclude false negatives
    let
        foo(a) = a ≤ 0 ? throw("a is $(a)") : a
        result = report_call(()->foo(0))
        @test !isempty(get_reports(result))
        @test first(get_reports(result)) isa UncaughtExceptionReport
    end

    # report even if there're other "critical" error exist
    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            foo(a) = sum(a) # should be reported
            bar(a) = throw(a) # shouldn't be reported first
            $report_call((Bool, String)) do b, s
                b && foo(s)
                bar(s)
            end
        end)
        @test length(get_reports(result)) === 3
        test_sum_over_string(get_reports(result))
    end

    # end to end
    let
        # this should report `throw(ArgumentError("Sampler for this object is not defined")`
        result = report_call(rand, (Char,))
        @test !isempty(get_reports(result))
        @test first(get_reports(result)) isa UncaughtExceptionReport

        # this should not report `throw(DomainError(x, "sin(x) is only defined for finite x."))`
        result = report_call(sin, (Int,))
        @test isempty(get_reports(result))

        # again, constant prop sometimes can exclude false negatives
        result = report_call(()->sin(Inf))
        @test !isempty(get_reports(result))
        @test first(get_reports(result)) isa UncaughtExceptionReport
    end
end

@testset "keyword argument methods" begin
    result = Core.eval(gen_virtual_module(), quote
        f(a; b = nothing, c = nothing) = return
        $report_call((Any,)) do b
            f(1; b)
        end
    end)
    @test isempty(get_reports(result))
end

@testset "don't early escape if type grows up to `Any`" begin
    vmod = gen_virtual_module()
    res = @analyze_toplevel context = vmod virtualize = false begin
        abstract type Foo end
        struct Foo1 <: Foo
            bar
        end
        struct Foo2 <: Foo
            baz # typo
        end
        struct Foo3 <: Foo
            qux # typo
        end
        bar(foo::Foo) = foo.bar

        f = rand(Bool) ? Foo1(1) : rand(Bool) ? Foo2(1) : Foo3(1)
        bar(f)
    end

    # `bar(::Foo1)` is valid but its return type is `Any`, and so we can't collect possible
    # error points for `bar(::Foo2)` and `bar(::Foo3)` if we bail out on `Any`-grew return type
    @test length(res.inference_error_reports) === 2
    @test any(res.inference_error_reports) do er
        return er isa NoFieldErrorReport &&
            er.typ === vmod.Foo2 &&
            er.name === :bar
    end
    @test any(res.inference_error_reports) do er
        return er isa NoFieldErrorReport &&
            er.typ === vmod.Foo3 &&
            er.name === :bar
    end
end

@testset "abstract_invoke" begin
    # non-`Type` `argtypes`
    result = report_call() do
        invoke(sin, :this_should_be_type, 1.0)
    end
    @test length(get_reports(result)) == 1
    r = first(get_reports(result))
    @test isa(r, InvalidInvokeErrorReport)

    # invalid `argtype`
    result = report_call() do
        Base.@invoke sin(1.0::Int)
    end
    @test length(get_reports(result)) == 1
    r = first(get_reports(result))
    @test isa(r, InvalidInvokeErrorReport)

    # don't report errors collected in `invoke`d functions
    result = @eval Module() begin
        foo(i::Integer) = throw(string(i))
        $report_call((Int,)) do a
            Base.@invoke foo(a::Integer)
        end
    end
    @test !isempty(get_reports(result))
    @test !any(r->isa(r, InvalidInvokeErrorReport), get_reports(result))

    #== LINE SENSITIVITY START ===#
    BAR_LINE = (@__LINE__) + 3
    BAZ_LINE = (@__LINE__) + 5
    result = @eval begin
        bar(a) = return a + undefvar
        function baz(a)
            a += 1
            Base.@invoke bar(a::Int) # `invoke` is valid, but error should happne within `bar`
        end
        $report_call(baz, (Any,))
    end
    #== LINE SENSITIVITY END ===#
    @test !isempty(get_reports(result))
    @test !any(r->isa(r, InvalidInvokeErrorReport), get_reports(result))
    # virtual stack trace should include frames of both `bar` and `baz`
    # we already `@assert` it within `typeinf` when `JET_DEV_MODE` is enabled,
    # but test it explicitly here just to make sure it's working
    @test all(get_reports(result)) do r
        any(r.vst) do vf
            vf.file === Symbol(@__FILE__) &&
            vf.line == BAR_LINE && # `bar`
            vf.linfo.def.name === :bar
        end || return false
        any(r.vst) do vf
            vf.file === Symbol(@__FILE__) &&
            vf.line == BAZ_LINE && # `baz`
            vf.linfo.def.name === :baz
        end || return false
        return true
    end
end

@testset "staged programming" begin
    m = @fixturedef begin
        @generated function foo(a)
            if a <: Integer
                return :(a)
            elseif a <: Number
                return :(undefvar) # report me this case
            end
            throw("invalid argument")
        end

        bar(args...) = foo(args...)
    end

    let # successful code generation, valid code
        result = report_call(m.foo, (Int,))
        @test isempty(get_reports(result))
    end

    let # successful code generation, invalid code
        result = report_call(m.foo, (Float64,))
        @test length(get_reports(result)) == 1
        r = first(get_reports(result))
        @test isa(r, GlobalUndefVarErrorReport)
        @test r.name === :undefvar
    end

    let # unsuccessful code generation
        result = report_call(m.foo, (String,))
        @test length(get_reports(result)) == 1
        r = first(get_reports(result))
        @test isa(r, GeneratorErrorReport) && r.err == "invalid argument"
    end

    let # should work if cached
        result = report_call(m.bar, (String,))
        @test length(get_reports(result)) == 1
        r = first(get_reports(result))
        @test isa(r, GeneratorErrorReport) && r.err == "invalid argument"
    end
end

@testset "report invalid builtin call" begin
    result = report_call((Int, Type{Int}, Any)) do a, b, c
        isa(a, b, c)
    end
    @test length(get_reports(result)) === 1
    report = first(get_reports(result))
    @test report isa InvalidBuiltinCallErrorReport &&
        widenconst.(report.argtypes) == [Int, Type{Int}, Any]

    @testset "constant propagation" begin
        m = gen_virtual_module()
        Core.eval(m, quote
            struct T
                v
            end
            access_field(t, sym) = getfield(t, sym)
        end)

        result = Core.eval(m, quote
            $report_call(t->access_field(t,:v), (T,))
        end)
        @test isempty(get_reports(result))

        result = Core.eval(m, quote
            $report_call(t->access_field(t,:w), (T,))
        end)
        @test length(get_reports(result)) === 1
        er = first(get_reports(result))
        @test er isa NoFieldErrorReport
        @test er.typ === m.T
        @test er.name === :w

        result = Core.eval(m, quote
            $report_call(t->access_field(t,:v), (T,))
        end)
        @test isempty(get_reports(result))
    end
end

@testset "malformed getfield" begin
    let
        # shouldn't error
        result = report_call((Any,)) do a
            getfield(a)
        end
        @test length(get_reports(result)) == 1
        @test first(get_reports(result)) isa InvalidBuiltinCallErrorReport
    end
end

@testset "special case `return_type`" begin
    # don't report invalid method calls simulated in `return_type_tfunc`
    let
        result = report_call(()->CC.return_type(sum, Tuple{String}))
        @test isempty(get_reports(result))
    end

    # report invalid call of `return_type` itself
    let
        result = report_call(()->CC.return_type(sum))
        @test length(get_reports(result)) == 1
        @test isa(first(get_reports(result)), InvalidReturnTypeCall)
    end

    # end to end
    let
        # this shouldn't report "no matching method found for call signature: Base.iterate(itr::DataType)",
        # which otherwise will be caught in `abstract_cal` in `return_type_tfunc`
        result = report_call(() -> Dict('a' => 1, :b => 2))
        @test isempty(get_reports(result))
    end
end

@testset "target_modules" begin
    # from `PrintConfig` docstring

    M = Module()
    @eval M begin
        function foo(a)
            r1 = sum(a)       # => Base: MethodError(+(::Char, ::Char)), MethodError(zero(::Type{Any}))
            r2 = undefsum(a)  # => @__MODULE__: UndefVarError(:undefsum)
            return r1, r2
        end
    end

    let
        result = @report_call M.foo("julia")
        @test length(get_reports(result)) == 3
        @test any(get_reports(result)) do report
            isa(report, GlobalUndefVarErrorReport) && report.name === :undefsum
        end
    end

    let
        result = @report_call target_modules=(M,) M.foo("julia")
        report = only(get_reports(result))
        @test isa(report, GlobalUndefVarErrorReport) && report.name === :undefsum
    end
end

@testset "BasicPass" begin
    @testset "basicfilter" begin
        # skip errors on abstract dispatch
        let # https://github.com/aviatesk/JET.jl/issues/154
            res = @analyze_toplevel analyze_from_definitions=true begin
                struct Foo
                    x::AbstractVector{<:AbstractString}
                end
            end
            @test isempty(res.inference_error_reports)
        end

        # should still report anything within entry frame
        let
            res = @eval Module() begin
                foo(a::Int) = "hello"
                $report_call((AbstractString,)) do a # this abstract call isn't concrete dispatch
                    foo(a)
                end
            end
            @test !isempty(get_reports(res))
            @test any(r->isa(r,NoMethodErrorReport), get_reports(res))
        end

        # skip errors on abstract entry frame entered by `analyze_from_definitions!`
        let
            res = @analyze_toplevel analyze_from_definitions=true begin
                struct Foo end
                function isfoo(x)
                    # ==(::Missing, ::Foo) -> Missing will lead to `NonBooleanCondErrorReport` otherwise
                    return x == Foo() ? :foo : :bar
                end
            end
            @test !any(res.inference_error_reports) do r
                isa(r, NonBooleanCondErrorReport)
            end
        end
    end
end

@testset "SoundPass" begin
    let # `:basicfilter`
        m = Module()
        # `==(::Missing, ::Any) -> Missing` will be reported
        @eval m foo(a, b) = a == b ? 0 : 1
        result = @eval m $report_call((Any,Symbol)) do a, b
            foo(a, b)
        end
        @test isempty(get_reports(result))

        result = report_call((Any,Symbol); mode = :sound) do a, b
            a == b ? 0 : 1
        end
        @test any(get_reports(result)) do r
            isa(r, NonBooleanCondErrorReport) &&
            r.t == Type[Missing]
        end
    end

    let # `!(t ⊑ Bool)` -> error
        basic = report_call((Integer,)) do cond
            cond ? 0 : 1
        end
        @test isempty(get_reports(basic))

        sound = report_call((Integer,); mode=:sound) do cond
            cond ? 0 : 1
        end
        @test any(get_reports(sound)) do r
            isa(r, NonBooleanCondErrorReport) &&
            r.t == Integer
        end
    end
end

@testset "TypoPass" begin
    @test_call mode=:typo sum("julia") # don't report NoMethodError, etc.

    @testset "global undef var" begin
        let
            result = report_call(()->foo; mode=:typo)
            @test length(get_reports(result)) === 1
            @test first(get_reports(result)) isa GlobalUndefVarErrorReport
            @test first(get_reports(result)).name === :foo
        end

        let # deeper level
            m = @fixturedef begin
                foo(bar) = bar + baz
                qux(a) = foo(a)
            end

            result = Core.eval(m, :($report_call(qux, (Int,); mode=:typo)))
            @test length(get_reports(result)) === 1
            @test first(get_reports(result)) isa GlobalUndefVarErrorReport
            @test first(get_reports(result)).name === :baz

            # works when cached
            result = Core.eval(m, :($report_call(qux, (Int,); mode=:typo)))
            @test length(get_reports(result)) === 1
            @test first(get_reports(result)) isa GlobalUndefVarErrorReport
            @test first(get_reports(result)).name === :baz
        end
    end

    @testset "no field" begin
        m = @fixturedef begin
            struct T
                v
            end
            access_field(t, sym) = getfield(t, sym)
        end

        result = Core.eval(m, :($report_call(t->access_field(t,:w), (T,); mode=:typo)))
        @test length(get_reports(result)) === 1
        er = first(get_reports(result))
        @test er isa NoFieldErrorReport
        @test er.typ === m.T
        @test er.name === :w
    end
end


text:
import Base: copyto!, copy, deepcopy_internal
# include("karray.jl") ## _unsafe_copy!, KnetArray

# Generalizing low level copy using linear indexing to/from gpu arrays:
# copyto!{T}(dest::Array{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer)
# Note that this is an unsafe operation, no argument or bounds checking performed.
# Defined in Base:
# _unsafe_copy!{T}(dest::Ptr{T}, src::Ptr{T}, n) at array.jl:73
# _unsafe_copy!{T}(dest::Array{T,N}, doffs, src::Array{T,N}, soffs, n) at array.jl:79

const KorA{T} = Union{KnetArray{T},Array{T}}

function _copyto!(dest::KorA{T}, doffs::Integer, src::KorA{T}, soffs::Integer, n::Integer) where {T}
    if n == 0; return dest; end
    if n < 0; throw(ArgumentError()); end
    if soffs < 1 || doffs < 1 || soffs+n-1 > length(src) || doffs+n-1 > length(dest)
        throw(BoundsError())
    end
    _unsafe_copy!(dest, doffs, src, soffs, n)
end

function _copyto!(dest::KorA{T}, src::KorA{T}) where {T}
    if length(dest) < length(src); throw(BoundsError()); end
    copyto!(dest, 1, src, 1, length(src))
end

# Avoid Array->Array to prevent base conflict.
copyto!(dest::KnetArray{T}, doffs::Integer, src::KnetArray{T}, soffs::Integer, n::Integer) where {T} = _copyto!(dest, doffs, src, soffs, n)
copyto!(dest::KnetArray{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer) where {T} = _copyto!(dest, doffs, src, soffs, n)
copyto!(dest::Array{T}, doffs::Integer, src::KnetArray{T}, soffs::Integer, n::Integer) where {T} = _copyto!(dest, doffs, src, soffs, n)
copyto!(dest::KnetArray{T}, src::KnetArray{T}) where {T} = _copyto!(dest, src)
copyto!(dest::KnetArray{T}, src::Array{T}) where {T} = _copyto!(dest, src)
copyto!(dest::Array{T}, src::KnetArray{T}) where {T} = _copyto!(dest, src)

function copy(a::KnetArray)
    _unsafe_copy!(similar(a),1,a,1,length(a))
end

# This will make deepcopy work properly
deepcopy_internal(x::KnetArray, s::IdDict)=if haskey(s,x); s[x]; else; copy(x); end


text:
# This file was generated, do not modify it. # hide
purchase = caravan.Purchase
vals     = unique(purchase)


text:
@inline function matmul_axes(C, A, B)
    M = indices((C,A), (StaticInt{1}(),StaticInt{1}()))
    K = indices((A,B), (StaticInt{2}(),StaticInt{1}()))
    N = indices((C,B), (StaticInt{2}(),StaticInt{2}()))
    M, K, N
end

@inline _select(::StaticInt{M}, ::StaticInt{M}) where {M} = StaticInt{M}()
@noinline _select(::StaticInt{M}, ::StaticInt{N}) where {M,N} = throw("$M ≠ $N")
@inline _select(::StaticInt{M}, _) where {M} = StaticInt{M}()
@inline _select(_, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline _select(x, _) = x
@inline function matmul_sizes(C,A,B)
    MC, NC = size(C)
    MA, KA = size(A)
    KB, NB = size(B)
    @assert ((MC == MA) & (KA == KB) & (NC == NB)) "Size mismatch."
    (_select(MA, MC), _select(KA, KB), _select(NB, NC))
end
# @inline function matmul_sizes(C, A, B)
#     M, K, N = matmul_axes(C, A, B)
#     static_length(M), static_length(K), static_length(N)
# end

@inline function loopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{0}, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ
    end
    nothing
end
@inline function loopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{1}, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] += Cₘₙ
    end
    nothing
end
@inline function loopmul!(C, A, B, ::StaticInt{1}, β, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ + β * C[m,n]
    end
    nothing
end
@inline function loopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{0}, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ -= A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ
    end
    nothing
end
@inline function loopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{1}, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ -= A[m,k] * B[k,n]
        end
        C[m,n] += Cₘₙ
    end
    nothing
end
@inline function loopmul!(C, A, B, ::StaticInt{-1}, β, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ -= A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ + β * C[m,n]
    end
    nothing
end
@inline function loopmul!(C, A, B, α, ::StaticInt{0}, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = α * Cₘₙ
    end
    nothing
end

@inline function loopmul!(C, A, B, α, ::StaticInt{1}, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] += α * Cₘₙ
    end
    nothing
end
@inline function loopmul!(C, A, B, α, β, (M, K, N))
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = α * Cₘₙ + β * C[m,n]
    end
    nothing
end

@inline function packamul!(
    C, Ãₚ, A, B,
    ::StaticInt{1}, ::StaticInt{0},
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ += Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] = Cₘₙ
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    ::StaticInt{1}, ::StaticInt{1},
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ += Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] += Cₘₙ
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    ::StaticInt{1}, β,
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ += Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] = Cₘₙ + β * C[m,n]
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    ::StaticInt{-1}, ::StaticInt{0},
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ -= Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] = Cₘₙ
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    ::StaticInt{-1}, ::StaticInt{1},
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ -= Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] += Cₘₙ
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    ::StaticInt{-1}, β,
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ -= Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] = Cₘₙ + β * C[m,n]
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    α, ::StaticInt{0},
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ += Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] = α * Cₘₙ
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    α, ::StaticInt{1},
    M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ += Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] += α * Cₘₙ
    end
end
@inline function packamul!(
    C, Ãₚ, A, B,
    α, β, M, K, N
)
    # @avx for n ∈ N, m ∈ M
    @avx for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Aₘₖ = A[m,k]
            Cₘₙ += Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] = α * Cₘₙ + β * C[m,n]
   end
end
@inline function alloc_a_pack(A, M, ::Type{T}) where {T}
    buffer = core_cache_buffer(T, Val(2))
    # _buffer = core_cache_buffer(T, Val(2))
    # buffer = Vector{T}(undef, length(_buffer))
    # Apack = ptrarray0(align(pointer(buffer)), (M, K), (static_sizeof(T), static_sizeof(T) * align(M, T)), DenseDims{(true,false)}())
    Apack = default_zerobased_stridedpointer(align(pointer(buffer)), (static_sizeof(T), static_sizeof(T) * align(M, T)))
    Apack, buffer
end

@inline function packaloopmul!(
    C::AbstractStridedPointer{T},
    A::AbstractStridedPointer,
    B::AbstractStridedPointer,
    α, β, (M, K, N)
) where {T}
    Ãₚ, buffer = alloc_a_pack(A, M, T)
    # @show M,K,N, threadid()
    # @show pointer.((C,A,B)), threadid()
    GC.@preserve buffer begin
        Nᵣ = StaticInt{nᵣ}()
        packamul!(C, Ãₚ, A, B, α, β, CloseOpen(M), CloseOpen(K), CloseOpen(Nᵣ))
        loopmul!(
            gesp(C, (Zero(), Nᵣ)), Ãₚ, gesp(B, (Zero(), Nᵣ)), α, β,
            (CloseOpen(M), CloseOpen(K), CloseOpen(N - Nᵣ))
        )
    end
    nothing
end

# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     ::StaticInt{1}, ::StaticInt{0},
# ) where {Mᵣ,Mᵢ,K,Nᵣ,Nᵢ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}()
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}()
#     @avx for mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aᵣₖᵢ = A[mᵣ,mᵢ,k]
#             Cₘₙ += Aᵣₖᵢ * B[nᵣ,k]
#             Ãₚ[mᵣ,k,mᵢ] = Aᵣₖᵢ
#         end
#         C[mᵣ,mᵢ,nᵣ] = Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,4))
#     @avx for nᵢ ∈ Nᵢrange, mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[mᵣ,k,mᵢ] * B[nᵣ,k,nᵢ]
#         end
#         C[mᵣ,mᵢ,nᵣ,nᵢ] = Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     ::StaticInt{1}, ::StaticInt{1},
# ) where {Mᵣ,Mᵢ,K,Nᵣ,Nᵢ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}()
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}()
#     @avx for mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aᵣₖᵢ = A[mᵣ,mᵢ,k]
#             Cₘₙ += Aᵣₖᵢ * B[nᵣ,k]
#             Ãₚ[mᵣ,k,mᵢ] = Aᵣₖᵢ
#         end
#         C[mᵣ,mᵢ,nᵣ] += Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,4))
#     @avx for nᵢ ∈ Nᵢrange, mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[mᵣ,k,mᵢ] * B[nᵣ,k,nᵢ]
#         end
#         C[mᵣ,mᵢ,nᵣ,nᵢ] += Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     ::StaticInt{1}, β
# ) where {Mᵣ,Mᵢ,K,Nᵣ,Nᵢ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}()
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}()
#     @avx for mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aᵣₖᵢ = A[mᵣ,mᵢ,k]
#             Cₘₙ += Aᵣₖᵢ * B[nᵣ,k]
#             Ãₚ[mᵣ,k,mᵢ] = Aᵣₖᵢ
#         end
#         C[mᵣ,mᵢ,nᵣ] = Cₘₙ + β * C[mᵣ,mᵢ,nᵣ]
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,4))
#     @avx for nᵢ ∈ Nᵢrange, mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[mᵣ,k,mᵢ] * B[nᵣ,k,nᵢ]
#         end
#         C[mᵣ,mᵢ,nᵣ,nᵢ] = Cₘₙ + β * C[mᵣ,mᵢ,nᵣ,nᵢ]
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     α, ::StaticInt{0},
# ) where {Mᵣ,Mᵢ,K,Nᵣ,Nᵢ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}()
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}()
#     @avx for mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aᵣₖᵢ = A[mᵣ,mᵢ,k]
#             Cₘₙ += Aᵣₖᵢ * B[nᵣ,k]
#             Ãₚ[mᵣ,k,mᵢ] = Aᵣₖᵢ
#         end
#         C[mᵣ,mᵢ,nᵣ] = α * Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,4))
#     @avx for nᵢ ∈ Nᵢrange, mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[mᵣ,k,mᵢ] * B[nᵣ,k,nᵢ]
#         end
#         C[mᵣ,mᵢ,nᵣ,nᵢ] = α * Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     α, ::StaticInt{1},
# ) where {Mᵣ,Mᵢ,K,Nᵣ,Nᵢ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}()
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}()
#     @avx for mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aᵣₖᵢ = A[mᵣ,mᵢ,k]
#             Cₘₙ += Aᵣₖᵢ * B[nᵣ,k]
#             Ãₚ[mᵣ,k,mᵢ] = Aᵣₖᵢ
#         end
#         C[mᵣ,mᵢ,nᵣ] += α * Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,4))
#     @avx for nᵢ ∈ Nᵢrange, mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[mᵣ,k,mᵢ] * B[nᵣ,k,nᵢ]
#         end
#         C[mᵣ,mᵢ,nᵣ,nᵢ] += α * Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     α, β
# ) where {Mᵣ,Mᵢ,K,Nᵣ,Nᵢ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}()
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}()
#     @avx for mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aᵣₖᵢ = A[mᵣ,mᵢ,k]
#             Cₘₙ += Aᵣₖᵢ * B[nᵣ,k]
#             Ãₚ[mᵣ,k,mᵢ] = Aᵣₖᵢ
#         end
#         C[mᵣ,mᵢ,nᵣ] = α * Cₘₙ + β * C[mᵣ,mᵢ,nᵣ]
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,4))
#     @avx for nᵢ ∈ Nᵢrange, mᵢ ∈ axes(Ãₚ,3), mᵣ ∈ Mᵣrange, nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[mᵣ,k,mᵢ] * B[nᵣ,k,nᵢ]
#         end
#         C[mᵣ,mᵢ,nᵣ,nᵢ] = α * Cₘₙ + β * C[mᵣ,mᵢ,nᵣ,nᵢ]
#     end
#     nothing
# end

# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{M,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{M,K}},
#     A::AbstractStrideArray{Tuple{M,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     ::StaticInt{1}, ::StaticInt{0}
# ) where {M,K,Nᵣ,Nᵢ}
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}();
#     @avx for m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[nᵣ,k]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,nᵣ] = Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,3))
#     @avx for nᵢ ∈ Nᵢrange, m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[m,k] * B[nᵣ,k,nᵢ]
#         end
#         C[m,nᵣ,nᵢ] = Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{M,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{M,K}},
#     A::AbstractStrideArray{Tuple{M,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     ::StaticInt{1}, ::StaticInt{1}
# ) where {M,K,Nᵣ,Nᵢ}
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}();
#     @avx for m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[nᵣ,k]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,nᵣ] += Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,3))
#     @avx for nᵢ ∈ Nᵢrange, m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[m,k] * B[nᵣ,k,nᵢ]
#         end
#         C[m,nᵣ,nᵢ] += Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{M,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{M,K}},
#     A::AbstractStrideArray{Tuple{M,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     ::StaticInt{1}, β
# ) where {M,K,Nᵣ,Nᵢ}
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}();
#     @avx for m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[nᵣ,k]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,nᵣ] = Cₘₙ + β * C[m,nᵣ]
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,3))
#     @avx for nᵢ ∈ Nᵢrange, m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[m,k] * B[nᵣ,k,nᵢ]
#         end
#         C[m,nᵣ,nᵢ] = Cₘₙ + β * C[m,nᵣ,nᵢ]
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{M,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{M,K}},
#     A::AbstractStrideArray{Tuple{M,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     α, ::StaticInt{0}
# ) where {M,K,Nᵣ,Nᵢ}
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}();
#     @avx for m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[nᵣ,k]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,nᵣ] = α * Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,3))
#     @avx for nᵢ ∈ Nᵢrange, m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[m,k] * B[nᵣ,k,nᵢ]
#         end
#         C[m,nᵣ,nᵢ] = α * Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{M,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{M,K}},
#     A::AbstractStrideArray{Tuple{M,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     α, ::StaticInt{1}
# ) where {M,K,Nᵣ,Nᵢ}
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}();
#     @avx for m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[nᵣ,k]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,nᵣ] += α * Cₘₙ
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,3))
#     @avx for nᵢ ∈ Nᵢrange, m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[m,k] * B[nᵣ,k,nᵢ]
#         end
#         C[m,nᵣ,nᵢ] += α * Cₘₙ
#     end
#     nothing
# end
# function packaloopmul!(
#     C::AbstractStrideArray{Tuple{M,Nᵣ,Nᵢ}},
#     Ãₚ::AbstractStrideArray{Tuple{M,K}},
#     A::AbstractStrideArray{Tuple{M,K}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K,Nᵢ}},
#     α, β
# ) where {M,K,Nᵣ,Nᵢ}
#     Nᵣrange = VectorizationBase.StaticUnitRange{1,Nᵣ}();
#     @avx for m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[nᵣ,k]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,nᵣ] = α * Cₘₙ + β * C[m,nᵣ]
#     end
#     Nᵢrange = VectorizationBase.StaticLowerUnitRange{2}(size(C,3))
#     @avx for nᵢ ∈ Nᵢrange, m ∈ axes(Ãₚ,1), nᵣ ∈ Nᵣrange
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(Ãₚ,2)
#             Cₘₙ += Ãₚ[m,k] * B[nᵣ,k,nᵢ]
#         end
#         C[m,nᵣ,nᵢ] = α * Cₘₙ + β * C[m,nᵣ,nᵢ]
#     end
#     nothing
# end

# function loopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K}},
#     ::StaticInt{1}, ::StaticInt{0}
# ) where {Mᵣ,Mᵢ,K,Nᵣ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}();
#     @avx for mᵢ ∈ axes(A,3), mᵣ ∈ Mᵣrange, nᵣ ∈ axes(C,3)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(A,2)
#             Cₘₙ += A[mᵣ,k,mᵢ] * B[nᵣ,k]
#         end
#         C[mᵣ,mᵢ,nᵣ] = Cₘₙ
#     end
#     nothing
# end
# function loopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K}},
#     ::StaticInt{1}, ::StaticInt{1}
# ) where {Mᵣ,Mᵢ,K,Nᵣ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}();
#     @avx for mᵢ ∈ axes(A,3), mᵣ ∈ Mᵣrange, nᵣ ∈ axes(C,3)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(A,2)
#             Cₘₙ += A[mᵣ,k,mᵢ] * B[nᵣ,k]
#         end
#         C[mᵣ,mᵢ,nᵣ] += Cₘₙ
#     end
#     nothing
# end
# function loopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K}},
#     ::StaticInt{1}, β
# ) where {Mᵣ,Mᵢ,K,Nᵣ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}();
#     @avx for mᵢ ∈ axes(A,3), mᵣ ∈ Mᵣrange, nᵣ ∈ axes(C,3)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(A,2)
#             Cₘₙ += A[mᵣ,k,mᵢ] * B[nᵣ,k]
#         end
#         C[mᵣ,mᵢ,nᵣ] = Cₘₙ + β * C[mᵣ,mᵢ,nᵣ]
#     end
#     nothing
# end
# function loopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K}},
#     α, ::StaticInt{0}
# ) where {Mᵣ,Mᵢ,K,Nᵣ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}();
#     @avx for mᵢ ∈ axes(A,3), mᵣ ∈ Mᵣrange, nᵣ ∈ axes(C,3)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(A,2)
#             Cₘₙ += A[mᵣ,k,mᵢ] * B[nᵣ,k]
#         end
#         C[mᵣ,mᵢ,nᵣ] = α * Cₘₙ
#     end
#     nothing
# end
# function loopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K}},
#     α, ::StaticInt{1}
# ) where {Mᵣ,Mᵢ,K,Nᵣ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}();
#     @avx for mᵢ ∈ axes(A,3), mᵣ ∈ Mᵣrange, nᵣ ∈ axes(C,3)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(A,2)
#             Cₘₙ += A[mᵣ,k,mᵢ] * B[nᵣ,k]
#         end
#         C[mᵣ,mᵢ,nᵣ] += α * Cₘₙ
#     end
#     nothing
# end
# function loopmul!(
#     C::AbstractStrideArray{Tuple{Mᵣ,Mᵢ,Nᵣ}},
#     A::AbstractStrideArray{Tuple{Mᵣ,K,Mᵢ}},
#     B::AbstractStrideArray{Tuple{Nᵣ,K}},
#     α, β
# ) where {Mᵣ,Mᵢ,K,Nᵣ}
#     Mᵣrange = VectorizationBase.StaticUnitRange{1,Mᵣ}();
#     @avx for mᵢ ∈ axes(A,3), mᵣ ∈ Mᵣrange, nᵣ ∈ axes(C,3)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ axes(A,2)
#             Cₘₙ += A[mᵣ,k,mᵢ] * B[nᵣ,k]
#         end
#         C[mᵣ,mᵢ,nᵣ] = α * Cₘₙ + β * C[mᵣ,mᵢ,nᵣ]
#     end
#     nothing
# end

@inline function inlineloopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{0})
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ
    end
    C
end
@inline function inlineloopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{1})
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] += Cₘₙ
    end
    C
end
@inline function inlineloopmul!(C, A, B, ::StaticInt{1}, β)
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ + β * C[m,n]
    end
    C
end
@inline function inlineloopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{0})
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ -= A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ
    end
    C
end
@inline function inlineloopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{1})
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ -= A[m,k] * B[k,n]
        end
        C[m,n] += Cₘₙ
    end
    C
end
@inline function inlineloopmul!(C, A, B, ::StaticInt{-1}, β)
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ -= A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ + β * C[m,n]
    end
    C
end
@inline function inlineloopmul!(C, A, B, α, ::StaticInt{0})
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = α * Cₘₙ
    end
    C
end
@inline function inlineloopmul!(C, A, B, α, ::StaticInt{1})
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] += α * Cₘₙ
    end
    C
end
@inline function inlineloopmul!(C, A, B, α, β)
    M, K, N = matmul_axes(C, A, B)
    # @avx inline=true for n ∈ N, m ∈ M
    @avx inline=true for m ∈ M, n ∈ N
        Cₘₙ = zero(eltype(C))
        for k ∈ K
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n]  = α * Cₘₙ + β * C[m,n]
    end
    C
end

function two_cols_per_vector_add_kn_block!(q::Expr, Kunroll, Nunroll, Koffset, Noffset, Ashuffle, Bshuffle, MindA, MindC, Amask)
    for kk ∈ 0:Kunroll-1
        kt = kk + Koffset
        Aload = Expr(:call, :vload, :ptrA, Expr(:tuple, MindA, kt))
        isnothing(Amask) || push!(Aload.args, Amask)
        push!(q.args, Expr(:(=), Symbol(:A_,kk), Expr(:call, :shufflevector, Expr(:call, :extract_data, Aload), Ashuffle)))
    end
    for kk ∈ 0:Kunroll-1, nn ∈ 0:Nunroll-1
        kt = kk + Koffset
        nt1 = 2(nn + Noffset)
        nt2 = nt1 + 1
        Bsym = Symbol(:B_,kk,:_,nn)
        Csym = Symbol(:C_, nn)
        # push!(q.args, Expr(:(=), Bsym, Expr(:call, :shufflevector, Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt)), Bshuffle)))
        shuffle = Expr(
            :tuple,
            Expr(:call, :(Core.VecElement), Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt1))),
            Expr(:call, :(Core.VecElement), Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt2)))
        )
        push!(q.args, Expr(:(=), Bsym, Expr(:call, :shufflevector, shuffle, Bshuffle)))
        if iszero(Koffset | kk)
            push!(q.args, Expr(:(=), Csym, Expr(:call, :vmul, Symbol(:A_, kk), Bsym)))
        else
            push!(q.args, Expr(:(=), Csym, Expr(:call, :vfmadd, Symbol(:A_, kk), Bsym, Csym)))
        end
    end
    for nn ∈ 0:Nunroll-1
        push!(q.args, Expr(:call, :vnoaliasstore!, :ptrC, Symbol(:C_, nn), Expr(:tuple, MindC, 2(nn + Noffset))))
    end
end

function initmulquote()
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :ptrC, Expr(:call, :stridedpointer, :C)),
        Expr(:(=), :ptrA, Expr(:call, :stridedpointer, :A)),
        Expr(:(=), :ptrB, Expr(:call, :stridedpointer, :B))
    )
end

function two_cols_per_vector_quote(K, N, W, Wshift, Amask = nothing)
    # Calculate 2 columns of C at a time
    q = two_cols_per_vector_quote!(initmulquote(), K, N, W, Wshift, 0, Amask)
    push!(q.args, :C)
    q
end

function two_cols_per_vector_quote!(q, K, N, W, Wshift, Noffbase = 0, Amask = nothing)
    Ndoublereps = N >> 1
    Wh = W >>> 1
    Whs = Wshift - 1

    Kloop = LoopVectorization.Loop(:k, 1, K, Symbol(""), Symbol(""), true, true)
    Nloop = LoopVectorization.Loop(:n, 1, Ndoublereps, Symbol(""), Symbol(""), true, true)
    cost_vec     = Float64[ 64.0, 32.0, 8.0, 0.0 ] # Values chosen to hopefully produce desired behavior... TODO: do this better?
    reg_pressure = Float64[ 1.0, 1.0, 1.0, 0.0, 0.0 ]
    Kunroll, Nunroll, cost = LoopVectorization.solve_unroll(
        :k, :n, cost_vec, reg_pressure, W, :m, Kloop, Nloop, 1, 1
    )
    @assert isfinite(cost)

    Nrep, Nrem = divrem(Ndoublereps, Nunroll)
    Krep, Krem = divrem(K, Kunroll)
    Ashuftup = Expr(:tuple); append!(Ashuftup.args, 0:Wh-1); append!(Ashuftup.args, 0:Wh-1)
    Bshuftup = Expr(:tuple); foreach(_ -> push!(Bshuftup.args, 0), 1:Wh); foreach(_ -> push!(Bshuftup.args, 1), 1:Wh)
    Ashuffle = Expr(:call, Expr(:curly, :Val, Ashuftup))
    Bshuffle = Expr(:call, Expr(:curly, :Val, Bshuftup))
    MindA = Expr(:call, Expr(:curly, :_MM, Wh), 0)
    MindC = Expr(:call, Expr(:curly, :_MM, W), 0)
    for n ∈ 0:Nrep-1
        for k ∈ 0:Krep-1
            two_cols_per_vector_add_kn_block!(q, Kunroll, Nunroll, Kunroll*k, Nunroll*n + Noffbase, Ashuffle, Bshuffle, MindA, MindC, Amask)
        end
        if Krem > 0
            two_cols_per_vector_add_kn_block!(q, Krem, Nunroll, Kunroll*Krep, Nunroll*n + Noffbase, Ashuffle, Bshuffle, MindA, MindC, Amask)
        end
    end
    if Nrem > 0
        for k in 0:Krep-1
            two_cols_per_vector_add_kn_block!(q, Kunroll, Nrem, Kunroll*k, Nunroll*Nrep + Noffbase, Ashuffle, Bshuffle, MindA, MindC, Amask)
        end
        if Krem > 0
            two_cols_per_vector_add_kn_block!(q, Krem, Nrem, Kunroll*Krep, Nunroll*Nrep + Noffbase, Ashuffle, Bshuffle, MindA, MindC, Amask)
        end
    end
    if isodd(N)
        Csym = Symbol(:C_,0)
        for k ∈ 0:K-1
            Asym = Symbol(:A_,k)
            Bsym = Symbol(:B_,k)
            Aload = Expr(:call, :vload, :ptrA, Expr(:tuple, MindA, k))
            isnothing(Amask) || push!(Aload.args, Amask)
            push!(q.args, Expr(:(=), Asym, Expr(:call, :extract_data, Aload)))
            push!(q.args, Expr(:(=), Bsym, Expr(:call, :vload, :ptrB, Expr(:tuple, k, N-1))))
            if iszero(k)
                push!(q.args, Expr(:(=), Csym, Expr(:call, :vmul, Asym, Bsym)))
            else
                push!(q.args, Expr(:(=), Csym, Expr(:call, :vfmadd, Asym, Bsym, Csym)))
            end
        end
        push!(q.args, Expr(:call, :vnoaliasstore!, :ptrC, Csym, Expr(:tuple, MindA, N-1))) #MindA, because we want half-vector
    end
    q
end

# function LinearAlgebra.mul!(
# @generated function LinearAlgebra.mul!(
#     C::AbstractFixedSizeMatrix{M,N,T,1,4,false},
#     A::AbstractFixedSizeMatrix{M,K,T,1,XA},
#     B::AbstractFixedSizeMatrix{K,N,T}
# ) where {M,K,N,XA,T}
# # ) where {M,K,N,T,XA}
#     W, Wshift = VectorizationBase.pick_vector_width_shift(4N, T)
#     nvectors_per_col = W >> 2
#     # if nvectors_per_col ≤ 1 || (K * N > 4VectorizationBase.REGISTER_COUNT) || XA < 4
#     if nvectors_per_col != 2 || (K * N > 4VectorizationBase.REGISTER_COUNT) || XA < 4
#         return Expr(:block, Expr(:meta,:inline), Expr(:call, :jmul!, :C, :A, :B))
#     elseif M == 3#if nvectors_per_col == 2
#         return two_cols_per_vector_quote(K, N, W, Wshift, 0x77)
#     else
#         return two_cols_per_vector_quote(K, N, W, Wshift)
#     # else
#         # return Expr(:block, Expr(:meta,:inline), Expr(:call, :jmul!, :C, :A, :B))
#     end
# end

function four_cols_per_vector_add_kn_block!(q::Expr, Kunroll, Nunroll, Koffset, Noffset, Ashuffle, Bshuffle, MindA, MindC)
    for kk ∈ 0:Kunroll-1
        kt = kk + Koffset
        push!(q.args, Expr(:(=), Symbol(:A_,kk), Expr(:call, :shufflevector, Expr(:call, :extract_data, Expr(:call, :vload, :ptrA, Expr(:tuple, MindA, kt))), Ashuffle)))
    end
    for kk ∈ 0:Kunroll-1, nn ∈ 0:Nunroll-1
        kt = kk + Koffset
        nt1 = 4(nn + Noffset)
        nt2 = nt1 + 1
        Bsym = Symbol(:B_,kk,:_,nn)
        Csym = Symbol(:C_, nn)
        # push!(q.args, Expr(:(=), Bsym, Expr(:call, :shufflevector, Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt)), Bshuffle)))
        shuffle = Expr(
            :tuple,
            Expr(:call, :(Core.VecElement), Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt1))),
            Expr(:call, :(Core.VecElement), Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt2))),
            Expr(:call, :(Core.VecElement), Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt2+1))),
            Expr(:call, :(Core.VecElement), Expr(:call, :vload, :ptrB, Expr(:tuple, kt, nt2+2)))
        )
        push!(q.args, Expr(:(=), Bsym, Expr(:call, :shufflevector, shuffle, Bshuffle)))
        if iszero(Koffset | kk)
            push!(q.args, Expr(:(=), Csym, Expr(:call, :vmul, Symbol(:A_, kk), Bsym)))
        else
            push!(q.args, Expr(:(=), Csym, Expr(:call, :vfmadd, Symbol(:A_, kk), Bsym, Csym)))
        end
    end
    for nn ∈ 0:Nunroll-1
        push!(q.args, Expr(:call, :vnoaliasstore!, :ptrC, Symbol(:C_, nn), Expr(:tuple, MindC, 4(nn + Noffset) )))
    end
end
function four_cols_per_vector_quote(K, N, W, Wshift)
    # Calculate 2 columns of C at a time
    Nquadreps = N >> 2
    Nquadrem = N & 3
    # Nisodd = N & 1
    q = initmulquote()
    Wq = W >>> 2
    Wqs = Wshift - 2

    Kloop = LoopVectorization.Loop(:k, 1, K, Symbol(""), Symbol(""), true, true)
    Nloop = LoopVectorization.Loop(:n, 1, Nquadreps, Symbol(""), Symbol(""), true, true)
    cost_vec     = Float64[ 64.0, 32.0, 8.0, 0.0 ] # Values chosen to hopefully produce desired behavior... TODO: do this better?
    reg_pressure = Float64[ 1.0, 1.0, 1.0, 0.0, 0.0 ]
    Kunroll, Nunroll, cost = LoopVectorization.solve_unroll(
        :k, :n, cost_vec, reg_pressure, W, :m, Kloop, Nloop
    )
    @assert isfinite(cost)

    Nrep, Nrem = divrem(Nquadreps, Nunroll)
    Krep, Krem = divrem(K, Kunroll)
    Ashuftup = Expr(:tuple);
    append!(Ashuftup.args, 0:Wq-1); append!(Ashuftup.args, 0:Wq-1)
    append!(Ashuftup.args, 0:Wq-1); append!(Ashuftup.args, 0:Wq-1)
    Bshuftup = Expr(:tuple);
    foreach(_ -> push!(Bshuftup.args, 0), 1:Wq); foreach(_ -> push!(Bshuftup.args, 1), 1:Wq)
    foreach(_ -> push!(Bshuftup.args, 0), 1:Wq); foreach(_ -> push!(Bshuftup.args, 1), 1:Wq)
    Ashuffle = Expr(:call, Expr(:curly, :Val, Ashuftup))
    Bshuffle = Expr(:call, Expr(:curly, :Val, Bshuftup))
    MindA = Expr(:call, Expr(:curly, :_MM, Wq), 0)
    MindC = Expr(:call, Expr(:curly, :_MM, W), 0)
    for n ∈ 0:Nrep-1
        for k ∈ 0:Krep-1
            four_cols_per_vector_add_kn_block!(q, Kunroll, Nunroll, Kunroll*k, Nunroll*n, Ashuffle, Bshuffle, MindA, MindC)
        end
        if Krem > 0
            four_cols_per_vector_add_kn_block!(q, Krem, Nunroll, Kunroll*Krep, Nunroll*n, Ashuffle, Bshuffle, MindA, MindC)
        end
    end
    if Nrem > 0
        for k in 0:Krep-1
            four_cols_per_vector_add_kn_block!(q, Kunroll, Nrem, Kunroll*k, Nunroll*Nrep, Ashuffle, Bshuffle, MindA, MindC)
        end
        if Krem > 0
            four_cols_per_vector_add_kn_block!(q, Krem, Nrem, Kunroll*Krep, Nunroll*Nrep, Ashuffle, Bshuffle, MindA, MindC)
        end
    end
    if Nquadrem > 0
        two_cols_per_vector_quote!(q, K, N, W >> 1, Wshift - 1, N & -4)
    end
    push!(q.args, :C)
    q
end

# function LinearAlgebra.mul!(
# @generated function LinearAlgebra.mul!(
#     C::AbstractFixedSizeMatrix{M,N,T,1,2,false},
#     A::AbstractFixedSizeMatrix{M,K,T,1,XA},
#     B::AbstractFixedSizeMatrix{K,N,T}
# ) where {M,K,N,XA,T}
# # ) where {M,K,N,T,XA}
#     W, Wshift = VectorizationBase.pick_vector_width_shift(2N, T)
#     nvectors_per_col = W >> 1
#     if nvectors_per_col ≤ 1 || (K * N > 4VectorizationBase.REGISTER_COUNT) || XA < 2
#     # if nvectors_per_col != 2 || (K * N > 4VectorizationBase.REGISTER_COUNT) || XA < 4
#         return Expr(:block, Expr(:meta,:inline), Expr(:call, :jmul!, :C, :A, :B))
#     elseif nvectors_per_col == 2
#         return two_cols_per_vector_quote(K, N, W, Wshift)
#     # elseif nvectors_per_col == 4
#         # return four_cols_per_vector_quote(K, N, W, Wshift)
#     else
#         return Expr(:block, Expr(:meta,:inline), Expr(:call, :jmul!, :C, :A, :B))
#     end
# end


text:
const IT = Interpolations

function pairs(curve::Contour.Curve2)
    collect(zip(Contour.coordinates(curve)...))
end

function LinearInterpolation(xs,ys,A)
    itp = IT.interpolate(A, IT.BSpline(IT.Linear(IT.OnGrid())))
    IT.scale(itp, xs, ys)
end

function CubicInterpolation(xs,ys,A)
    itp = IT.interpolate(A, IT.BSpline(IT.Cubic(IT.Line(IT.OnGrid()))))
    IT.scale(itp, xs, ys)
end

const _DEFAULT_ISOLINE_KWARGS =
    OrderedDict(:colors => ["MidnightBlue",
                            "LightSeaGreen",
                            "Gold",
                            "Tomato"],
                :lift => false,
                :interpolation => :cubic,
		:levels => 10)

"""
    isolines(xs, ys, zs; lift = false, interpolation = :cubic)

Plot the contour lines of the function whose values are represented
by the array (or function) `zs`. If `lift` is true, plot in 3D.

# Examples
```julia-repl
julia> isolines(0:10, 0:10, (x,y) -> (100 - x^2 + y^2)/10, lift = true)
"""
function isolines(args...;kwargs...)
    D = deepcopy(_DEFAULT_ISOLINE_KWARGS)
    merge!(D,Dict(kwargs))
    isolinekwargs, penkwargs = splitkwargs(kwargs,_DEFAULT_ISOLINE_KWARGS)
    grlist = GraphicElement[]
    if length(args) ≥ 3 && args[3] isa Function
        args = (args[1],
                args[2],
                float([args[3](x,y) for x in args[1], y in args[2]]))
    end
    isolinedict = Dict(isolinekwargs)
    C = Contour.contours(args...,isolinedict[:levels])
    m, M = extrema([Contour.level(cl) for cl in Contour.levels(C)])
    if D[:lift]
        if D[:interpolation] == :linear
            interp = LinearInterpolation(args[1:3]...)
        elseif D[:interpolation] == :cubic
            interp = CubicInterpolation(args[1:3]...)
        else
            error("`interpolation` should be either :cubic or :linear")
        end
    end
    for cl in Contour.levels(C)
        L = Contour.level(cl)
        for line in Contour.lines(cl)
            pairs(line)
            pts = (D[:lift] ?
                   [(x,y,interp(x,y)) for (x,y) in pairs(line)] :
                   pairs(line))
            push!(grlist,
                  Path(pts; pen =
                       Pen(color=AsyPlots.cmap(NamedColor.(D[:colors]),
                                               (L-m)/(M-m));penkwargs...)))
        end
    end
    Plot(grlist)
end


text:
using PlanetarySystem
using RecursiveArrayTools
using Test

# The package is tested on random Solar-like systems with a central start initialized at x=(0,0,0) and v=(0,0,0).
# The test verifies the conservation of energy and angular momentum along the evolution, within chosen tolerance.

tol=10^-7 # tolerance

@time @testset "Energy and angular momentum conservation" begin

    for j in 1:5  # loop over different random systems

        # initial conditions:
        N=5 # number of bodies
        M=rand(0.:0.000001:0.0001,N) # mass
        M[1]=1. # central star mass
        pos_x=rand(-30.:0.0001:30.,N) # initial positions
        pos_y=rand(-30.:0.0001:30.,N)
        pos_z=rand(-30.:0.0001:30.,N)
        pos_x[1]=0. # central star
        pos_y[1]=0.
        pos_z[1]=0.
        pos = ArrayPartition(pos_x,pos_y,pos_z)
        vel_x=rand(-0.01:0.00001:0.01,N) # initial velocities
        vel_y=rand(-0.01:0.00001:0.01,N)
        vel_z=rand(-0.01:0.00001:0.01,N)
        vel = ArrayPartition(vel_x,vel_y,vel_z)
        vel_x[1]=0. # central star
        vel_y[1]=0.
        vel_z[1]=0.
        tspan = (0.,200_000.) # time interval

        # solution:
        sol=NBsolution(M, vel, pos, tspan)
        # solution partition:
        vx=sol[1:N,:]
        vy=sol[N+1:2*N,:]
        vz=sol[2*N+1:3*N,:]
        x= sol[3*N+1:4*N,:]
        y= sol[4*N+1:5*N,:]
        z= sol[5*N+1:6*N,:]

        # conservation test:
        if(length(sol.t)>=5000) # test step check

            for i in 1:5000:length(sol.t) # loop over the time steps; test step=500
                @test (H(x, y, z, vx, vy, vz, M)[1] -  H(x, y, z, vx, vy, vz, M)[i])/H(x, y, z, vx, vy, vz, M)[1] ≈ 0 atol=tol #energy
                @test (Lx(x, y, z, vx, vy, vz, M)[1] - Lx(x, y, z, vx, vy, vz, M)[i])/Lx(x, y, z, vx, vy, vz, M)[1] ≈ 0 atol=tol #angular momentum x
                @test (Ly(x, y, z, vx, vy, vz, M)[1] - Ly(x, y, z, vx, vy, vz, M)[i])/Ly(x, y, z, vx, vy, vz, M)[1]  ≈ 0 atol=tol #angular momentum y
                @test (Lz(x, y, z, vx, vy, vz, M)[1] - Lz(x, y, z, vx, vy, vz, M)[i])/Lz(x, y, z, vx, vy, vz, M)[1] ≈ 0 atol=tol #angular momentum z
            end

        else
            println("reduce test step!") # alert
        end

     end

end


text:
@static if VERSION >= v"0.7.0-DEV.2915"

    using Base64
    using Distributed
    import Dates

    compat_search(a...) = (r = search(a...); r === nothing ? 0 : r)
    compat_findfirst(a...) = (r = findfirst(a...); r === nothing ? 0 : r)
    compat_findprev(a...) = (r = findprev(a...); r === nothing ? 0 : r)
    compat_findnext(a...) = (r = findnext(a...); r === nothing ? 0 : r)

else # Julia v0.6

    eval(:(module Base64 end))
    const Dates = Base.Dates
    const Distributed = Base.Distributed

    const compat_search = search
    const compat_findfirst = findfirst
    const compat_findprev = findprev
    const compat_findnext = findnext

    pairs(x) = [k => v for (k,v) in x]

    Base.SubString(s) = SubString(s, 1)

    macro debug(s) DEBUG_LEVEL > 0 ? :(("D- ", $(esc(s)))) : :() end
    macro info(s)  DEBUG_LEVEL >= 0 ? :(println("I- ", $(esc(s)))) : :() end
    macro warn(s)  DEBUG_LEVEL >= 0 ? :(println("W- ", $(esc(s)))) : :() end
    macro error(m, args...)
        args = [:(print("|  "); @show $a) for a in args]
        esc(:(println("E- ", $m); $(args...); nothing))
    end
end

macro uninit(expr)
    if !isdefined(Base, :uninitialized)
        splice!(expr.args, 2)
    end
    return esc(expr)
end

if !isdefined(Base, :Nothing)
    const Nothing = Void
    const Cvoid = Void
end

# https://github.com/JuliaLang/julia/pull/25535
Base.String(x::SubArray{UInt8,1}) = String(Vector{UInt8}(x))


text:
using Turing
using Turing.RandomMeasures
using Test
using Random

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "RandomMeasures.jl" begin
    partitions = [
        [[1, 2, 3, 4]],
        [[1, 2, 3], [4]],
        [[1, 2, 4], [3]],
        [[1, 2], [3, 4]],
        [[1, 2], [3], [4]],
        [[1, 3, 4], [2]],
        [[1, 3], [2, 4]],
        [[1, 3], [2], [4]],
        [[1, 4], [2, 3]],
        [[1], [2, 3, 4]],
        [[1], [2, 3], [4]],
        [[1, 4], [2], [3]],
        [[1], [2, 4], [3]],
        [[1], [2], [3, 4]],
        [[1], [2], [3], [4]]]

    @turing_testset "chinese restaurant processes" begin
        # Data
        data = [-2,2,-1.5,1.5]

        # Base distribution
        mu_0 = mean(data)
        sigma_0 = 4
        sigma_1 = 0.5
        tau0 = 1/sigma_0^2
        tau1 = 1/sigma_1^2

        # DP parameters
        alpha = 0.25

        @model crpimm(y, rpm) = begin
            # Base distribution.
            H = Normal(mu_0, sigma_0)
            # Latent assignments.
            N = length(y)
            z = tzeros(Int, N)

            # Cluster counts.
            cluster_counts = tzeros(Int, N)

            # Cluster locations.
            x = tzeros(Float64, N)
            for i in 1:N
                # Draw assignments using a CRP.
                z[i] ~ ChineseRestaurantProcess(rpm, cluster_counts)
                if cluster_counts[z[i]] == 0
                    # Cluster is new, therefore, draw new location.
                    x[z[i]] ~ H
                end
                cluster_counts[z[i]] += 1
                # Draw observation.
                y[i] ~ Normal(x[z[i]], sigma_1)
            end
        end

        rpm = DirichletProcess(alpha)

        sampler = SMC(1000)
        mf = crpimm(data, rpm)

        # Compute empirical posterior distribution over partitions
        samples = sample(mf, sampler)

        # Check that there is no NaN value associated
        z_samples = Int.(samples[:z].value)
        @test all(!isnan, samples[:x].value[z_samples])
        @test all(!ismissing, samples[:x].value[z_samples])

        empirical_probs = zeros(length(partitions))
        w = map(x -> x.weight, samples.info[:samples])
        sum_weights = sum(w)
        z = z_samples

        for i in 1:size(z,1)
            partition = map(c -> findall(z[i,:,1] .== c), unique(z[i,:,1]))
            partition_idx = findfirst(p -> sort(p) == sort(partition), partitions)
            @test partition_idx != nothing
            empirical_probs[partition_idx] += sum_weights == 0 ? 1 : w[i]
        end

        if sum_weights == 0
            empirical_probs /= length(w)
        end

        l2, discr = correct_posterior(empirical_probs, data, partitions, tau0, tau1, alpha, 1e-7)

        @test l2 < 0.05
        @test discr < 0.2
    end
    @testset "distributions" begin
        @turing_testset "Representations" begin
            d = StickBreakingProcess(DirichletProcess(1.0))
            @test minimum(d) == 0
            @test maximum(d) == 1

            d = SizeBiasedSamplingProcess(DirichletProcess(1.0), 2.0)
            @test minimum(d) == 0
            @test maximum(d) == 2

            d = ChineseRestaurantProcess(DirichletProcess(1.0), [2, 1])
            @test minimum(d) == 1
            @test maximum(d) == 3
        end
        @turing_testset "Dirichlet Process" begin

            α = 0.1
            N = 10_000

            # test SB representation
            d = StickBreakingProcess(DirichletProcess(α))
            Ev = mapreduce(_ -> rand(d), +, 1:N) / N
            @test Ev ≈ mean(Beta(1, α)) atol=0.05

            # test SBS representation
            d = SizeBiasedSamplingProcess(DirichletProcess(α), 2.0)
            Ej = mapreduce(_ -> rand(d), +, 1:N) / N
            @test Ej ≈ mean(Beta(1, α)) * 2 atol=0.05

            # test CRP representation
            d = ChineseRestaurantProcess(DirichletProcess(α), [2, 1])
            ks = map(_ -> rand(d), 1:N)
            c = map(k -> sum(ks .== k), support(d))
            p = c ./ sum(c)

            q = [2, 1, α] ./ (2 + α)
            q ./= sum(q)

            @test p[1] ≈ q[1] atol=0.1
            @test p[2] ≈ q[2] atol=0.1
            @test p[3] ≈ q[3] atol=0.1
        end
        @turing_testset "Pitman-Yor Process" begin

            a = 0.5
            θ = 0.1
            t = 2
            N = 10_000

            # test SB representation
            d = StickBreakingProcess(PitmanYorProcess(a, θ, t))
            Ev = mapreduce(_ -> rand(d), +, 1:N) / N
            @test Ev ≈ mean(Beta(1 - a, θ + t*a)) atol=0.05

            # test SBS representation
            d = SizeBiasedSamplingProcess(PitmanYorProcess(a, θ, t), 2.0)
            Ej = mapreduce(_ -> rand(d), +, 1:N) / N
            @test Ej ≈ mean(Beta(1 - a, θ + t*a)) * 2 atol=0.05

            # test CRP representation
            d = ChineseRestaurantProcess(PitmanYorProcess(a, θ, t), [2, 1])
            ks = map(_ -> rand(d), 1:N)
            c = map(k -> sum(ks .== k), support(d))
            p = c ./ sum(c)

            q = [2 - a, 1 - a, θ + a*t] ./ (3 + θ)
            q ./= sum(q)

            @test p[1] ≈ q[1] atol=0.1
            @test p[2] ≈ q[2] atol=0.1
            @test p[3] ≈ q[3] atol=0.1
        end
    end
    @turing_testset "stick breaking" begin
        # Data
        data = [-2,2,-1.5,1.5]

        # Base distribution
        mu_0 = mean(data)
        sigma_0 = 4
        sigma_1 = 0.5
        tau0 = 1/sigma_0^2
        tau1 = 1/sigma_1^2

        # DP parameters
        alpha = 0.25

        # stick-breaking process based on Papaspiliopoulos and Roberts (2008).
        @model sbimm(y, rpm, trunc) = begin
            # Base distribution.
            H = Normal(mu_0, sigma_0)

            # Latent assignments.
            N = length(y)
            z = tzeros(Int, N)

            # Infinite (truncated) collection of breaking points on unit stick.
            v = tzeros(Float64, trunc)

            # Cluster locations.
            x = tzeros(Float64, trunc)

            # Draw weights and locations.
            for k in 1:trunc
                v[k] ~ StickBreakingProcess(rpm)
                x[k] ~ H
            end

            # Weights.
            w = vcat(v[1], v[2:end] .* cumprod(1 .- v[1:end-1]))

            # Normalize weights to ensure they sum exactly to one.
            # This is required by the Categorical distribution in Distributions.
            w ./= sum(w)

            for i in 1:N
                # Draw location
                z[i] ~ Categorical(w)

                # Draw observation.
                y[i] ~ Normal(x[z[i]], sigma_1)
            end
        end

        rpm = DirichletProcess(alpha)

        sampler = SMC(1000)
        mf = sbimm(data, rpm, 10)

        # Compute empirical posterior distribution over partitions
        samples = sample(mf, sampler)

        # Check that there is no NaN value associated
        z_samples = Int.(samples[:z].value)
        @test all(!isnan, samples[:x].value[z_samples])
        @test all(!ismissing, samples[:x].value[z_samples])

        empirical_probs = zeros(length(partitions))
        w = map(x -> x.weight, samples.info[:samples])
        sum_weights = sum(w)
        z = z_samples

        for i in 1:size(z,1)
            partition = map(c -> findall(z[i,:,1] .== c), unique(z[i,:,1]))
            partition_idx = findfirst(p -> sort(p) == sort(partition), partitions)
            @test partition_idx != nothing
            empirical_probs[partition_idx] += sum_weights == 0 ? 1 : w[i]
        end

        if sum_weights == 0
            empirical_probs /= length(w)
        end

        l2, discr = correct_posterior(empirical_probs, data, partitions, tau0, tau1, alpha, 1e-7)

        # Increased ranges due to truncation error.
        @test l2 < 0.1
        @test discr < 0.3
    end
    @turing_testset "size-based sampling" begin
        # Data
        data = [-2,2,-1.5,1.5]

        # Base distribution
        mu_0 = mean(data)
        sigma_0 = 4
        sigma_1 = 0.5
        tau0 = 1/sigma_0^2
        tau1 = 1/sigma_1^2

        # DP parameters
        alpha = 0.25

        # size-biased sampling process
        @model sbsimm(y, rpm, trunc) = begin
            # Base distribution.
            H = Normal(mu_0, sigma_0)

            # Latent assignments.
            N = length(y)
            z = tzeros(Int, N)

            x = tzeros(Float64, N)
            J = tzeros(Float64, N)
            z = tzeros(Int, N)

            k = 0
            surplus = 1.0

            for i in 1:N
                ps = vcat(J[1:k], surplus)
                z[i] ~ Categorical(ps)
                if z[i] > k
                    k = k + 1
                    J[k] ~ SizeBiasedSamplingProcess(rpm, surplus)
                    x[k] ~ H
                    surplus -= J[k]
                end
                y[i] ~ Normal(x[z[i]], sigma_1)
            end
        end

        rpm = DirichletProcess(alpha)

        sampler = SMC(1000)
        mf = sbsimm(data, rpm, 100)

        # Compute empirical posterior distribution over partitions
        samples = sample(mf, sampler)

        # Check that there is no NaN value associated
        z_samples = Int.(samples[:z].value)
        @test all(!isnan, samples[:x].value[z_samples])
        @test all(!ismissing, samples[:x].value[z_samples])

        empirical_probs = zeros(length(partitions))
        w = map(x -> x.weight, samples.info[:samples])
        sum_weights = sum(w)
        z = z_samples

        for i in 1:size(z,1)
            partition = map(c -> findall(z[i,:,1] .== c), unique(z[i,:,1]))
            partition_idx = findfirst(p -> sort(p) == sort(partition), partitions)
            @test partition_idx != nothing
            empirical_probs[partition_idx] += sum_weights == 0 ? 1 : w[i]
        end

        if sum_weights == 0
            empirical_probs /= length(w)
        end

        l2, discr = correct_posterior(empirical_probs, data, partitions, tau0, tau1, alpha, 1e-7)

        @test l2 < 0.05
        @test discr < 0.2
    end
end


text:
using Statistics
using StatsBase: sample
using Revise
using EvoTrees
using BenchmarkTools

# prepare a dataset
features = rand(Int(1.25e6), 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

params_c = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=1.0, colsample=1.0, nbins=32);

model_c, cache_c = EvoTrees.init_evotree(params_c, X_train, Y_train);

# initialize from cache
params_c = model_c.params
train_nodes = cache_c.train_nodes
splits = cache_c.splits
X_size = size(cache_c.X_bin)

# select random rows and cols
𝑖 = cache_c.𝑖_[sample(params.rng, cache_c.𝑖_, ceil(Int, params_c.rowsample * X_size[1]), replace=false, ordered=true)]
𝑗 = cache_c.𝑗_[sample(params.rng, cache_c.𝑗_, ceil(Int, params_c.colsample * X_size[2]), replace=false, ordered=true)]
# reset gain to -Inf
for feat in cache_c.𝑗_
    splits[feat].gain = -Inf
end

# build a new tree
# 897.800 μs (6 allocations: 736 bytes)
@btime EvoTrees.update_grads!(params_c.loss, params_c.α, cache_c.pred, cache_c.Y, cache_c.δ, cache_c.δ², cache_c.𝑤)
∑δ, ∑δ², ∑𝑤 = sum(cache_c.δ[𝑖]), sum(cache_c.δ²[𝑖]), sum(cache_c.𝑤[𝑖])
gain = EvoTrees.get_gain(params_c.loss, ∑δ, ∑δ², ∑𝑤, params_c.λ)
# assign a root and grow tree
train_nodes[1] = EvoTrees.TrainNode(0, 1, ∑δ, ∑δ², ∑𝑤, gain, 𝑖, 𝑗)
# 69.247 ms (1852 allocations: 38.41 MiB)
@btime tree = grow_tree(cache_c.δ, cache_c.δ², cache_c.𝑤, cache_c.hist_δ, cache_c.hist_δ², cache_c.hist_𝑤, params_c, train_nodes, splits, cache_c.edges, cache_c.X_bin);
push!(model_c.trees, tree)
@btime EvoTrees.predict!(cache_c.pred, tree, cache_c.X)

###########################
# Tree CPU
###########################
δ, δ², 𝑤, hist_δ, hist_δ², hist_𝑤, edges, X_bin = cache_c.δ, cache_c.δ², cache_c.𝑤, cache_c.hist_δ, cache_c.hist_δ², cache_c.hist_𝑤, cache_c.edges, cache_c.X_bin;

T = Float32
L = 1
active_id = ones(Int, 1)
leaf_count = one(Int)
tree_depth = one(Int)
tree = EvoTrees.Tree(Vector{EvoTrees.TreeNode{L,T,Int,Bool}}())

id = 1
node = train_nodes[id]
# 9.613 ms (81 allocations: 13.55 KiB)
@btime EvoTrees.update_hist!(hist_δ[id], hist_δ²[id], hist_𝑤[id], δ, δ², 𝑤, X_bin, node)
# 601.685 ns (6 allocations: 192 bytes) 8 100 feat ~ 60us
@btime EvoTrees.find_split!(view(hist_δ[id], :, j), view(hist_δ²[id], :, j), view(hist_𝑤[id], :, j), params, node, splits[j], edges[j])

###################################################
# GPU
###################################################
params_g = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=1.0, colsample=1.0, nbins=32);

model_g, cache_g = EvoTrees.init_evotree_gpu(params_g, X_train, Y_train);

params_g = model_g.params;
train_nodes = cache_g.train_nodes;
splits = cache_g.splits;
X_size = size(cache_g.X_bin);

# select random rows and cols
𝑖 = cache_g.𝑖_[sample(params_g.rng, cache_g.𝑖_, ceil(Int, params_g.rowsample * X_size[1]), replace=false, ordered=true)]
𝑗 = cache_g.𝑗_[sample(params_g.rng, cache_g.𝑗_, ceil(Int, params_g.colsample * X_size[2]), replace=false, ordered=true)]
# reset gain to -Inf
for feat in cache_g.𝑗_
    splits[feat].gain = -Inf
end

# build a new tree
# 164.200 μs (33 allocations: 1.66 KiB) - 5-6 X time faster on GPU
@btime CUDA.@sync EvoTrees.update_grads_gpu!(params_g.loss, cache_g.δ, cache_g.δ², cache_g.pred, cache_g.Y, cache_g.𝑤)
# sum Gradients of each of the K parameters and bring to CPU
∑δ, ∑δ², ∑𝑤 = Vector(vec(sum(cache_g.δ[𝑖,:], dims=1))), Vector(vec(sum(cache_g.δ²[𝑖,:], dims=1))), sum(cache_g.𝑤[𝑖])
gain = EvoTrees.get_gain(params_g.loss, ∑δ, ∑δ², ∑𝑤, params_g.λ)
# assign a root and grow tree
train_nodes[1] = EvoTrees.TrainNode_gpu(0, 1, ∑δ, ∑δ², ∑𝑤, gain, 𝑖, 𝑗)
# 60.736 ms (108295 allocations: 47.95 MiB) - only 15% faster than CPU
@btime CUDA.@sync tree = EvoTrees.grow_tree_gpu(cache_g.δ, cache_g.δ², cache_g.𝑤, cache_g.hist_δ, cache_g.hist_δ², cache_g.hist_𝑤, params_g, cache_g.K, train_nodes, splits, cache_g.edges, cache_g.X_bin, cache_g.X_bin_cpu);
push!(model_g.trees, tree);
# 2.736 ms (93 allocations: 13.98 KiB)
@btime CUDA.@sync EvoTrees.predict_gpu!(cache_g.pred_cpu, tree, cache_g.X)
# 1.013 ms (37 allocations: 1.19 KiB)
@btime CUDA.@sync cache_g.pred .= CuArray(cache_g.pred_cpu);

###########################
# Tree GPU
###########################
δ, δ², 𝑤, hist_δ, hist_δ², hist_𝑤, K, edges, X_bin, X_bin_cpu = cache_g.δ, cache_g.δ², cache_g.𝑤, cache_g.hist_δ, cache_g.hist_δ², cache_g.hist_𝑤, cache_g.K, cache_g.edges, cache_g.X_bin, cache_g.X_bin_cpu;
T = Float32
active_id = ones(Int, 1)
leaf_count = one(Int)
tree_depth = one(Int)
tree = EvoTrees.Tree_gpu(Vector{EvoTrees.TreeNode_gpu{T,Int,Bool}}())

hist_δ_cpu = zeros(T, size(hist_δ[1]));
hist_δ²_cpu = zeros(T, size(hist_δ²[1]));
hist_𝑤_cpu = zeros(T, size(hist_𝑤[1]));

id = 1
node = train_nodes[id];
# 7.003 ms (106 allocations: 3.09 KiB)
@btime CUDA.@sync EvoTrees.update_hist_gpu!(hist_δ[id], hist_δ²[id], hist_𝑤[id], δ, δ², 𝑤, X_bin, CuVector(node.𝑖), CuVector(node.𝑗), K);
#  32.001 μs (2 allocations: 32 bytes) * 3 adds 100us
@btime hist_δ_cpu .= hist_δ[id];
j = 1
# 2.925 μs (78 allocations: 6.72 KiB) * 100 features ~ 300us
@btime EvoTrees.find_split_gpu!(view(hist_δ_cpu, :, :, j), view(hist_δ²_cpu, :, :, j), view(hist_𝑤_cpu, :, j), params, node, splits[j], edges[j])


text:
# Julia wrapper for header: /Users/dls/jl/openssl-1.0.1e/include/openssl/mdc2.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 MDC2_Init (Ptr{MDC2_CTX},) libcrypto
@c Int32 MDC2_Update (Ptr{MDC2_CTX}, Ptr{Uint8}, Csize_t) libcrypto
@c Int32 MDC2_Final (Ptr{Uint8}, Ptr{MDC2_CTX}) libcrypto
@c Ptr{Uint8} MDC2 (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto


text:
# Configuration for running spatial simulation
#=
Recommended command line to run:
>  julia run.jl examples/example1
or
>  julia -p 4 run.jl examples/example1
=#
export simtype
@everywhere simtype = 2
const num_trials = 1
@everywhere const N_list = [50]        # Meta-population size
N_mut_list = [1.0]
#mutation_stddev_list = [0.005,0.01,0.02,0.04]
const num_subpops = 1                    # Number of subpopulations
const num_attributes_list = [1]        # number attributes for quantitative representation
const ngens = 1       # Generations after burn-in
const burn_in= 3.0    # generations of burn_in as a multiple of N
const ideal = 1.0
const neutral = true
const fit_slope=1.0
const additive_error=false
const wrap_attributes=false


text:
using CryptoUtils
using Test

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