text:
#=
Find nodes for tet elements.
=#
include("tet_quadrature_table.jl");

# Build global nodes for a triangular element with refel and vertices v
function tetrahedron_element_nodes(refel, v)
    return  tetrahedron_refel_to_xyz(refel.r[:,1], refel.r[:,2], refel.r[:,3], v);
end

# Set up refel nodal array
function tetrahedron_refel_nodes!(refel)
    (eqx, eqy, eqz) = tetrahedron_equilateral_nodes(refel.N);
    (r, s, t) = tetrahedron_equilateral_to_rst(eqx,eqy,eqz);

    # Adjust things within 1e-15 of the -1 boundaries to be exactly -1 and close to 0 to be 0
    tet_snap_to_bdry!(r, s, t, 1e-15);

    refel.r = zeros(refel.Np,3);
    refel.wr = zeros(refel.Np);

    refel.r[:,1] = r;
    refel.r[:,2] = s;
    refel.r[:,3] = t;

    # quadrature nodes/weights from a table
    xyzw = tetrahedron_quadrature_nodes_weights(refel.N+1);
    refel.g = xyzw[:,1:3];
    refel.wg = xyzw[:,4];

    # face node maps
    tol = 1e-12;
    tf1(x) = abs(x[1] + 1) < tol;
    tf2(x) = abs(x[2] + 1) < tol;
    tf3(x) = abs(x[3] + 1) < tol;
    tf4(x) = abs(x[1] + x[2] + x[3] + 1) < tol;
    refel.face2local = [get_face2local_map(refel.r, tf1),
                        get_face2local_map(refel.r, tf2),
                        get_face2local_map(refel.r, tf3),
                        get_face2local_map(refel.r, tf4)];

end

# Purpose  : Compute (x,y,z) nodes in equilateral tetrahedron for polynomial of order N
function tetrahedron_equilateral_nodes(N)
    # optimal alpha values for N up to 16 (from book)
    alpopt = [0, 0, 0, 0.1002, 1.1332, 1.5608, 1.3413, 1.2577, 1.1603, 1.10153, 0.6080, 0.4523, 0.8856, 0.8717, 0.9655];

    # Set optimized parameter, alpha, depending on order N
    if (N<16)
        alpha = alpopt[N];
    else
        alpha = 1;
    end

    # total number of nodes
    Np = Int64((N+1)*(N+2)*(N+3)/6);

    # Create equidistributed nodes on equilateral tetrahedron
    L1 = zeros(Np);
    L2 = zeros(Np);
    L3 = zeros(Np);
    L4 = zeros(Np);
    sk = 1;
    for n=1:(N+1)
        for m=1:(N+2-n)
            for q=1:(N+3-n-m)
                L1[sk] = (n-1)/N; # (1 + (-1 + (n-1)*2/N))/2
                L2[sk] = (m-1)/N; # (1 + (-1 + (m-1)*2/N))/2
                L3[sk] = 1 - (q+m+n-3)/N; # -(1 + -1 + (q-1)*2/N + -1 + (m-1)*2/N + -1 + (n-1)*2/N)/2
                L4[sk] = (q-1)/N; # (1 + (-1 + (q-1)*2/N))/2
                sk = sk+1;
            end
        end
    end

    # set vertices of tetrahedron
    invsqrt3 = 0.5773502691896258; # 1/sqrt(3)=0.5773502691896258
    invsqrt6 = 0.4082482904638631; # 1/sqrt(6)=0.4082482904638631
    v1 = [-1, -invsqrt3,   -invsqrt6];
    v2 = [ 1, -invsqrt3,   -invsqrt6];
    v3 = [ 0,  2*invsqrt3, -invsqrt6];
    v4 = [ 0,  0,          3*invsqrt6];

    # orthogonal axis tangents on faces 1-4
    t1 = zeros(4,3);
    t2 = zeros(4,3)
    t1[1,:] = v2.-v1;
    t1[2,:] = v2.-v1;
    t1[3,:] = v3.-v2;
    t1[4,:] = v3.-v1;
    t2[1,:] = v3.-0.5 .* (v1.+v2);
    t2[2,:] = v4.-0.5 .* (v1.+v2);
    t2[3,:] = v4.-0.5 .* (v2.+v3);
    t2[4,:] = v4.-0.5 .* (v1.+v3);

    for n=1:4 # normalize tangents
        t1[n,:] = t1[n,:] ./ norm(t1[n,:]);
        t2[n,:] = t2[n,:] ./ norm(t2[n,:]);
    end

    # Warp and blend for each face (accumulated in shiftXYZ)
    XYZ = zeros(Np,3);
    shift = zeros(Np,3);
    for i=1:Np
        for j=1:3
            XYZ[i,j] = L3[i]*v1[j] + L4[i]*v2[j] + L2[i]*v3[j] + L1[i]*v4[j]; # form undeformed coordinates
        end
    end

    for face=1:4
        if face==1     La = L1; Lb = L2; Lc = L3; Ld = L4;
        elseif face==2 La = L2; Lb = L1; Lc = L3; Ld = L4;
        elseif face==3 La = L3; Lb = L1; Lc = L4; Ld = L2;
        elseif face==4 La = L4; Lb = L1; Lc = L3; Ld = L2;
        end

        # compute warp tangential to face
        (warp1, warp2) = tetrahedron_warpfactor(N, Np, alpha, Lb, Lc, Ld);

        blend = Lb.*Lc.*Ld;   # compute volume blending

        denom = (Lb .+ .5 .* La).*(Lc .+ .5 .* La).*(Ld .+ .5 .* La);   # modify linear blend
        tol = 1e-10
        for ids=1:Np
            if denom[ids] > tol
                blend[ids] = (1 .+ (alpha.*La[ids]).^2).*blend[ids]./denom[ids];
            end
        end

        # compute warp & blend
        shift = shift .+ (blend.*warp1) * t1[face,:]' .+ (blend.*warp2) * t2[face,:]';

        # fix face warp
        for ids=1:Np
            numtrue = 0;
            if Lb[ids]>tol
                numtrue += 1;
            end
            if Lc[ids]>tol
                numtrue += 1;
            end
            if Ld[ids]>tol
                numtrue += 1;
            end
            if La[ids]<tol && (numtrue < 3)
                shift[ids,:] = warp1[ids]*t1[face,:] + warp2[ids]*t2[face,:];
            end
        end
    end

    XYZ = XYZ .+ shift;

    return (XYZ[:,1], XYZ[:,2], XYZ[:,3]);
end

# Compute scaled warp function at order N based on rout interpolation nodes
function tetrahedron_warpfactor(N, Np, alpha, L1, L2, L3)
    # Compute LGL and equidistant node distribution
    (LGLr, w) = jacobi_LGL_quad(N);

    # 2) compute blending function at each node for each edge
    blend1 = L2.*L3;
    blend2 = L1.*L3;
    blend3 = L1.*L2;

    # 3) amount of warp for each node, for each edge
    warp1 = zeros(Np);
    warp2 = zeros(Np);
    warp3 = zeros(Np);

    dx = zeros(Np);
    dy = zeros(Np);

    xeq = zeros(N+1);
    for i=1:N+1
        xeq[i] = -1 + 2*(N+1-i)/N;
    end

    for wi=1:Np
        for i=1:N+1
            d1 = -LGLr[i] - xeq[i];
            d2 = -LGLr[i] - xeq[i];
            d3 = -LGLr[i] - xeq[i];
            for j=2:N
                if i!=j
                    d1 = d1.*(L3[wi]-L2[wi]-xeq[j])/(xeq[i]-xeq[j]);
                    d2 = d2.*(L1[wi]-L3[wi]-xeq[j])/(xeq[i]-xeq[j]);
                    d3 = d3.*(L2[wi]-L1[wi]-xeq[j])/(xeq[i]-xeq[j]);
                end
            end

            if i!=1
                d1 = -d1/(xeq[i]-xeq[1]);
                d2 = -d2/(xeq[i]-xeq[1]);
                d3 = -d3/(xeq[i]-xeq[1]);
            end

            if i!=(N+1)
                d1 = d1/(xeq[i]-xeq[N+1]);
                d2 = d2/(xeq[i]-xeq[N+1]);
                d3 = d3/(xeq[i]-xeq[N+1]);
            end

            warp1[wi] = warp1[wi]+d1;
            warp2[wi] = warp2[wi]+d3;
            warp3[wi] = warp3[wi]+d3;
        end

        # 4) combine blend & warp
        warp1[wi] = 4 * blend1[wi] * warp1[wi] * (1 + (alpha * L1[wi])^2);
        warp2[wi] = 4 * blend2[wi] * warp2[wi] * (1 + (alpha * L2[wi])^2);
        warp3[wi] = 4 * blend3[wi] * warp3[wi] * (1 + (alpha * L3[wi])^2);

        # 5) evaluate shift in equilateral triangle
        # dx[wi] = 1*warp1[wi] + cos(2*pi/3)*warp2[wi] + cos(4*pi/3)*warp3[wi];
        # dy[wi] = 0*warp1[wi] + sin(2*pi/3)*warp2[wi] + sin(4*pi/3)*warp3[wi];
        dx[wi] = warp1[wi] - 0.5*warp2[wi] - 0.5*warp3[wi];
        dy[wi] = 0.8660254037844387*warp2[wi] - 0.8660254037844387*warp3[wi];
    end

    return (dx, dy);
end

# From (x,y,z) in equilateral tetrahedron to (r,s,t) coordinates in standard tetrahedron
function tetrahedron_equilateral_to_rst(x,y,z)
    invsqrt3 = 0.5773502691896258; # 1/sqrt(3)=0.5773502691896258
    invsqrt6 = 0.4082482904638631; # 1/sqrt(6)=0.4082482904638631
    v1 = [-1, -invsqrt3,   -invsqrt6];
    v2 = [ 1, -invsqrt3,   -invsqrt6];
    v3 = [ 0,  2*invsqrt3, -invsqrt6];
    v4 = [ 0,  0,          3*invsqrt6];

    # back out right tet nodes
    rhs = zeros(3,length(x));
    rhs[1,:] = x .- (0.5 * (v2[1]+v3[1]+v4[1]-v1[1]));
    rhs[2,:] = y .- (0.5 * (v2[2]+v3[2]+v4[2]-v1[2]));
    rhs[3,:] = z .- (0.5 * (v2[3]+v3[3]+v4[3]-v1[3]));

    A = [0.5 .* (v2 .- v1) 0.5 .* (v3 .- v1) 0.5 .* (v4 .- v1)];
    RST = A \ rhs;

    return (RST[1,:], RST[2,:], RST[3,:]);
end

# # From (r,s,t) coordinates in reference tetrahedron to (x,y,z) in tetrahedron with vertices v
# # v is a 3x4 array [x1 x2 x3 x4; y1 y2 y3 y4; z1 z2 z3 z4]
# function tetrahedron_refel_to_xyz(r, s, t, v)
#     d = v[:,1];
#     A = [v[:,2].-d v[:,3].-d v[:,4].-d];

#     np = length(r);
#     mv = zeros(3,np);
#     for i=1:np
#         tmp = [(r[i]+1)/2, (s[i]+1)/2, (t[i]+1)/2];
#         mv[:,i] = A*tmp + d;
#     end

#     x = mv[1,:]
#     y = mv[2,:]
#     z = mv[3,:]

#     return (x, y, z);
# end

# function get_face2local_map(r, compare)
#     n = size(r,1);
#     map = zeros(Int,n);
#     nf = 0;
#     for i=1:n
#         if compare(r[i,:])
#             nf = nf+1;
#             map[nf] = i;
#         end
#     end

#     return map[1:nf];
# end

function tet_snap_to_bdry!(r, s, t, tol)
    n = length(r);
    for i=1:n
        if abs(r[i]+1) < tol
            r[i] = -1;
        end
        if abs(s[i]+1) < tol
            s[i] = -1;
        end
        if abs(t[i]+1) < tol
            t[i] = -1;
        end

        if abs(r[i]) < tol
            r[i] = 0;
        end
        if abs(s[i]) < tol
            s[i] = 0;
        end
        if abs(t[i]) < tol
            t[i] = 0;
        end
    end
end


text:
using Test

using CUDAdrv
include(joinpath(@__DIR__, "..", "test", "array.jl"))   # real applications: use CuArrays.jl

dev = CuDevice(0)
ctx = CuContext(dev)

md = CuModuleFile(joinpath(@__DIR__, "vadd.ptx"))
vadd = CuFunction(md, "kernel_vadd")

dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)
c = similar(a)

d_a = CuTestArray(a)
d_b = CuTestArray(b)
d_c = CuTestArray(c)

len = prod(dims)
cudacall(vadd, Tuple{CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}}, d_a, d_b, d_c; threads=len)

@test a+b ≈ Array(d_c)

CUDAdrv.unsafe_destroy!(ctx)

ext: jl


text:
@ghdef mutable struct Installation
    id::Union{Int, Nothing}
end

namefield(inst::Installation) = inst.id

"""
    GitHub.create_access_token([api,] inst::Installation;
        auth::JWTAuth=jwt, options...)

Create an `OAuth2` access token for the `App` authenticated by `jwt` to make authenticated API
requests to the installation `inst` on an organization or individual account.

- https://developer.github.com/v3/apps/#create-a-new-installation-token
"""
@api_default function create_access_token(api::GitHubAPI, inst::Installation; auth::JWTAuth, headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.machine-man-preview+json"
    payload = gh_post_json(api, "/app/installations/$(inst.id)/access_tokens", auth = auth,
        headers=headers, options...)
    OAuth2(payload["token"])
end
@api_default function create_access_token(api::GitHubAPI, inst::Installation, auth::JWTAuth; options...)
    create_access_token(api, inst; auth=auth, options...)
end

"""
    GitHub.installation([api,] repo::Repo;
        auth::JWTAuth=jwt, options...)

Get the `Installation` on `repo` of the `App` authenticated by `jwt`.

- https://developer.github.com/v3/apps/#get-a-repository-installation
"""
@api_default function installation(api::GitHubAPI, repo::Repo; auth::JWTAuth, headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.machine-man-preview+json"
    result = gh_get_json(api, "/repos/$(name(repo))/installation", auth = auth,
        headers=headers, options...)
    Installation(result)
end
@api_default function installation(api::GitHubAPI, repo::Repo, auth::JWTAuth; options...)
    # TODO: deprecate
    installation(api, repo; auth=auth, options...)
end

"""
    insts, pagedata = GitHub.installations([api];
        auth::JWTAuth=jwt, options...)

List all `Installation`s corresponding to the app authenticated by `jwt`.

- https://developer.github.com/v3/apps/#list-installations
"""
@api_default function installations(api::GitHubAPI; auth::JWTAuth, headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.machine-man-preview+json"
    results, page_data = gh_get_paged_json(api, "/app/installations";
        auth=auth, headers=headers, options...)
    map(Installation, results), page_data
end
@api_default function installations(api::GitHubAPI, auth::JWTAuth; options...)
    # TODO: deprecate
    installations(api; auth=auth, options...)
end

"""
    repos, pagedata = GitHub.repos([api];
        auth::OAuth=token, options...)

List `Repo`s that the installation corresponding to `token` can access.

- https://developer.github.com/v3/apps/installations/#list-repositories
"""
@api_default function repos(api::GitHubAPI; auth::OAuth2, headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.machine-man-preview+json"
    results, page_data = github_paged_get(api, "/installation/repositories";
        auth=auth, headers=headers, options...)
    mapreduce(x->map(Repo, JSON.parse(HTTP.payload(x, String))["repositories"]), vcat, results; init=Repo[]), page_data
end
@api_default function repos(api::GitHubAPI, inst::Installation; options...)
    # TODO: deprecate
    repos(api; options...)
end


text:
function nearest_unoccupied_planet(ship::Ship, planets)
    d = Inf
    p = nothing
    for planet in planets
        if isowned(planet)
            if planet.owner_id == ship.owner_id && isfull(planet)
                continue
            end
        end
        dist = distance_between(ship, planet)
        if dist < d
            p = planet
            d = dist
        end
    end
    return p
end

function closest_point_to(e1::Entity, e2::Entity, min_distance::Float64 = 3.0)
    angle = angle_between(e1, e2)
    r = radius(e2) + min_distance
    x = e2.x - r * cos(deg2rad(angle))
    y = e2.y - r * sin(deg2rad(angle))
    return Position(x, y)
end

"""
    obstacles_between(game_map, entity1, entity2, ignore_ships, ignore_planets)
Check whether there is a straight-line path between two entities, without planetary or ships obstacles in between. Returns list of obstacles, empty if there is no obstacles were found.
"""
function obstacles_between(game_map, self, other, ignore_ships, ignore_planets)
    obstacles = Vector{Entity}()
    !ignore_ships   && check_collisions!(obstacles, all_ships(game_map), self, other)
    !ignore_planets && check_collisions!(obstacles, all_planets(game_map), self, other)
    return obstacles
end

function check_collisions!(obstacles, entitylist, self, other)
    for entity in entitylist
        any([self, other] .== entity) && continue
        if intersect_segment_circle(self, other, entity, radius(self) + 0.1)
            push!(obstacles, entity)
        end
    end
end

"""
intersect_segment_circle(start, end, circle, fudge=0.5)

Test whether a line segment and circle intersect.

# Arguments
- `start::Entity`: The start of the line segment. (Needs x, y attributes)
- `end::Entity`: The end of the line segment. (Needs x, y attributes)
- `circle::Entity`: The circle to test against. (Needs x, y, r attributes)
- `fudge::Entity`: A fudge factor; additional distance to leave between the segment and circle. (Probably set this to the ship radius, 0.5.)
"""
function intersect_segment_circle(startl::Entity, endl::Entity, circle::Entity, fudge=0.5)
    # Derived with SymPy
    # Parameterize the segment as start + t * (end - start),
    # and substitute into the equation of a circle
    # Solve for t
    dx = endl.x - startl.x
    dy = endl.y - startl.y

    a = dx^2 + dy^2
    # Start and end are the same point
    a < 1e-6 && return distance_between(startl, circle) <= radius(circle) + fudge

    b = -2.0 * (startl.x^2 - startl.x*endl.x - startl.x*circle.x + endl.x*circle.x +
                startl.y^2 - startl.y*endl.y - startl.y*circle.y + endl.y*circle.y)

    # Time along segment when closest to the circle (vertex of the quadratic)
    t = min(-b / (2 * a), 1.0)
    t < 0 && return false

    closest_x = startl.x + dx * t
    closest_y = startl.y + dy * t
    closest_distance = distance_between(Position(closest_x, closest_y), circle)

    return closest_distance <= radius(circle) + fudge
end

max_stars_repo_name: Cerbes/Halite-II

max_issues_repo_name: Cerbes/Halite-II

max_forks_repo_name: Cerbes/Halite-II


text:
using JLD, PyPlot

for m in 25
    for sgn in 0
        for numChanges in [1, 3, 5]
            file = jldopen("res/power_$(m)_$(sgn)_$(numChanges).jld", "r")
            power = read(file, "power")
            close(file)

            fig = figure(figsize=(3, 3), dpi=300)

            scatter(0.0:0.05:0.4, power[1:9, 1], s=5, marker="o")
            scatter(0.0:0.05:0.4, power[1:9, 2], s=5, marker="v")

            minorticks_on()
            grid()
            grid(which="minor", ls="dotted", lw=".5")

            xlabel("δ", size="xx-small")
            ylabel("proportion of rejections", size="xx-small")

            ylim(0, 1.05)

            tick_params("both", labelsize="xx-small")

            tight_layout()

            savefig("res/power_$(m)_$(sgn)_$(numChanges).png")

            close(fig)
        end
    end
end

for m in 25
    for sgn in 0
        for numChanges in [1, 3, 5]
            file = jldopen("res/power_$(m)_$(sgn)_$(numChanges).jld", "r")
            power = read(file, "power")
            close(file)

            fig = figure(figsize=(3, 3), dpi=300)

            scatter(0.0:0.05:0.4, power[1:9, 3], s=5, marker="o")
            scatter(0.0:0.05:0.4, power[1:9, 4], s=5, marker="v")

            minorticks_on()
            grid()
            grid(which="minor", ls="dotted", lw=".5")

            xlabel("δ", size="xx-small")
            ylabel("proportion of rejections", size="xx-small")

            ylim(0, 1.05)

            tick_params("both", labelsize="xx-small")

            tight_layout()

            savefig("res/power_studentized_$(m)_$(sgn)_$(numChanges).png")

            close(fig)
        end
    end
end

for m in 25
    for sgn in 0
        fig = figure(figsize=(27, 9), dpi=300)

        k = 0
        for numChanges in [1, 3, 5]
            for lbInd = 1:9
                k += 1

                ax = subplot(3, 9, k)

                file = jldopen("res/coverage_$(m)_$(sgn)_$(numChanges)_$(lbInd).jld", "r")
                coverage = read(file, "coverage")
                close(file)

                plot(0.0:0.05:1.0, 0.0:0.05:1.0, color="grey", linestyle=":", linewidth=.25)

                plot(0.05:0.05:0.95, coverage[1:19, 1], linewidth=.25, marker="o", markersize=2)
                plot(0.05:0.05:0.95, coverage[1:19, 2], linewidth=.25, marker="v", markersize=2)

                ax.set_xlim(0.0, 1.0)
                ax.set_ylim(0.0, 1.0)

                ax[:tick_params]("both", labelsize=2)
            end
        end

        tight_layout()

        savefig("res/coverage_$(m)_$(sgn).png")

        close(fig)
    end
end

for m in 25
    for sgn in 0
        fig = figure(figsize=(27, 9), dpi=300)

        k = 0
        for numChanges in [1, 3, 5]
            for lbInd = 1:9
                k += 1

                ax = subplot(3, 9, k)

                file = jldopen("res/coverage_$(m)_$(sgn)_$(numChanges)_$(lbInd).jld", "r")
                coverage = read(file, "coverage")
                close(file)

                plot(0.0:0.05:1.0, 0.0:0.05:1.0, color="grey", linestyle=":", linewidth=.25)

                plot(0.05:0.05:0.95, coverage[1:19, 3], linewidth=.25, marker="o", markersize=2)
                plot(0.05:0.05:0.95, coverage[1:19, 4], linewidth=.25, marker="v", markersize=2)

                ax.set_xlim(0.0, 1.0)
                ax.set_ylim(0.0, 1.0)

                ax[:tick_params]("both", labelsize=2)
            end
        end

        tight_layout()

        savefig("res/coverage_studentized_$(m)_$(sgn).png")

        close(fig)
    end
end


text:
# TODO add docs!!!

abstract type AbstractBitPrimitive{J} <: ArrowVector{J} end
export BitPrimitive

struct BitPrimitive <: AbstractBitPrimitive{Bool}
    length::Int
    values::Primitive{UInt8}  # in principle this can be anything, but don't see the need
end
export BitPrimitive

function BitPrimitive(data::Vector{UInt8}, i::Integer, len::Integer)
    BitPrimitive(len, Primitive{UInt8}(data, i, bytesforbits(len)))
end
function BitPrimitive(data::Vector{UInt8}, i::Integer, x::AbstractVector{Bool})
    BitPrimitive(len, Primitive{UInt8}(data, i, bitpack(x)))
end
function BitPrimitive(data::Vector{UInt8}, i::Integer, x::AbstractVector{UInt8},
                      len::Integer=sizeof(UInt8)*length(x))
    BitPrimitive(len, Primitive{UInt8}(data, i, x))
end

function BitPrimitive(v::AbstractVector{UInt8}, len::Integer=sizeof(UInt8)*length(v))
    BitPrimitive(len, Primitive{UInt8}(v))
end
function BitPrimitive(v::AbstractVector{Bool})
    BitPrimitive(length(v), Primitive{UInt8}(bitpack(v)))
end

function BitPrimitive(::Type{<:Array}, v::AbstractVector{UInt8}, len::Integer=sizeof(UInt8)*length(v))
    BitPrimitive(length(v), Primitive{UInt8}(Array, bitpack(v)))
end

BitPrimitive(b::BitPrimitive) = BitPrimitive(b.length, b.values)

struct NullableBitPrimitive <: AbstractBitPrimitive{Union{Bool,Missing}}
    length::Int
    bitmask::Primitive{UInt8}
    values::Primitive{UInt8}
end
export NullableBitPrimitive

function NullableBitPrimitive(data::Vector{UInt8}, bitmask_idx::Integer, values_idx::Integer,
                              len::Integer)
    bmask = Primitive{UInt8}(data, bitmask_idx, bitmaskbytes(len))
    vals = Primitive{UInt8}(data, values_idx, bytesforbits(len))
    NullableBitPrimitive(len, bmask, vals)
end

function NullableBitPrimitive(data::Vector{UInt8}, bitmask_idx::Integer, values_idx::Integer,
                              x::AbstractVector{T}) where T<:Union{Bool,Union{Bool,Missing}}
    bmask = Primitive{UInt8}(data, bitmask_idx, bitmaskpadded(x))
    vals = Primitive{UInt8}(data, values_idx, bytesforbits(x))
    setnonmissing!(vals, bitpack(x))
    NullableBitPrimitive(length(x), bmask, vals)
end
function NullableBitPrimitive(data::Vector{UInt8}, i::Integer, x::AbstractVector{T}
                             ) where T<:Union{Bool,Union{Bool,Missing}}
    NullableBitPrimitive(data, i, i+bitmaskbytes(x), x)
end

function NullableBitPrimitive(v::AbstractVector{T}) where T<:Union{Bool,Union{Bool,Missing}}
    bmask = Primitive{UInt8}(bitmaskpadded(v))
    vals = Primitive{UInt8}(bitpack(v))
    NullableBitPrimitive(length(v), bmask, vals)
end

function NullableBitPrimitive(::Type{<:Array}, v::AbstractVector{T}
                             ) where T<:Union{Bool,Union{Bool,Missing}}
    b = Vector{UInt8}(bitmaskbytes(v)+bytesforbits(v))
    NullableBitPrimitive(b, 1, v)
end

NullableBitPrimitive(b::NullableBitPrimitive) = NullableBitPrimitive(b.length, b.bitmask, b.values)

unsafe_getvalue(A::AbstractBitPrimitive, idx) = unsafe_getbit(values(A), idx)

getvalue(A::AbstractBitPrimitive, idx) = getbit(values(A), idx)
setvalue!(A::AbstractBitPrimitive, x::Bool, i::Integer) = setbit!(values(A), x, i)
function setvalue!(A::AbstractBitPrimitive, x::AbstractVector{Bool}, idx::AbstractVector{<:Integer})
    setbit!(values(A), x, idx)
end


text:
using Test, Serialization, Documenter, Pickle, PyCall, SparseArrays

DocMeta.setdocmeta!(Pickle, :DocTestSetup, :(using Pickle); recursive=true)

include("./pyscript.jl")
include("./torch/thscript.jl")

builtin_type_samples = Dict(
  "str" => "Julia!",
  "int" => 42,
  "bool" => Dict(
    true => false,
    false => true,
  ),
  "float" => 3.1415926,
  "bytes" => b"1234",
  "tuple" => (1, 2.0, "3", b"4"),
  "set" => Set((1,2,3,12,21)),
  "bigint" => 1234567890987654321012345678909876543210,
  "list" => ["February", 14, 2012],
)

const tests = [
  "load",
  "store",
  "np",
  "sparse",
]

Pickle.BATCHSIZE[] = 3

@testset "Pickle" begin
  @info "BATCHSIZE is set to: $(Pickle.BATCHSIZE[])"
  @info "Test doctest"
  doctest(Pickle)

  for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    @info "Test $(uppercase(t))"
    include(fp)
  end

  @info "Test Pickle.Torch"
  include("./torch/torch.jl")

end


text:
#adapted from osmnx: https://github.com/gboeing/osmnx/blob/main/osmnx/simplification.py

"""
Predicate wether v is an edge endpoint in the simplified version of g
"""
function is_endpoint(g::AbstractGraph, v)
    neighbors = all_neighbors(g, v)
    if v in neighbors # has self loop
        return true
    elseif outdegree(g, v) == 0 || indegree(g, v) == 0 # sink or source
        return true
    elseif length(neighbors) != 2 || indegree(g, v) != outdegree(g, v) # change to/from one way
        return true
    end
    return false
end

"""
iterator over all endpoints in g
"""
endpoints(g::AbstractGraph) = (v for v in vertices(g) if is_endpoint(g, v))

"""
iterator over all paths in g which can be contracted
"""
function paths_to_reduce(g::AbstractGraph)
    (path_to_endpoint(g, (u, v)) for u in endpoints(g) for v in outneighbors(g, u))
end

"""
path to the next endpoint starting in edge (ep, ep_succ)
"""
function path_to_endpoint(g::AbstractGraph, (ep, ep_succ)::Tuple{T,T}) where {T<:Integer}
    path = [ep, ep_succ]
    head = ep_succ
    # ep_succ not in endpoints -> has 2 neighbors and degree 2 or 4
    while !is_endpoint(g, head)
        neighbors = [n for n in outneighbors(g, head) if n != path[end-1]]
        @assert length(neighbors) == 1 "found unmarked endpoint!"
        head, = neighbors
        push!(path, head)
        (head == ep) && return path # self loop
    end
    return path
end

"""
Return the total weight of a path given as a Vector of Ids.
"""
function total_weight(g::OSMGraph, path::Vector{<:Integer})
    sum((g.weights[path[[i, i+1]]...] for i in 1:length(path)-1))
end

function ways_in_path(g::OSMGraph, path::Vector{<:Integer})
    ways = Set{Int}()
    for i in 1:(length(path)-1)
        edge = [g.index_to_node[path[i]], g.index_to_node[path[i+1]]]
        push!(ways, g.edge_to_way[edge])
    end
    return collect(ways)
end

"""
Build a new graph which simplifies the topology of osmg.graph.
The resulting graph only contains intersections and dead ends from the original graph.
The geometry of the contracted nodes is kept in the edge_gdf DataFrame
"""
function simplify_graph(osmg::OSMGraph{U, T, W}) where {U, T, W}
    g = osmg.graph
    relevant_nodes = collect(endpoints(g))
    n_relevant = length(relevant_nodes)
    graph = DiGraph(n_relevant)
    weights = similar(osmg.weights, (n_relevant, n_relevant))
    node_coordinates = Vector{Vector{W}}(undef, n_relevant)
    node_to_index = OrderedDict{T,U}()
    index_to_node = OrderedDict{U,T}()

    index_mapping = Dict{U,U}()
    for (new_i, old_i) in enumerate(relevant_nodes)
        index_mapping[old_i] = new_i
        node_coordinates[new_i] = osmg.node_coordinates[old_i]
        node = osmg.index_to_node[old_i]
        index_to_node[new_i] = node
        node_to_index[node] = new_i
    end

    edges = Dict{NTuple{3,U}, Vector{U}}()
    edge_count = Dict{Tuple{U,U}, Int}()
    for path in paths_to_reduce(g)
        u = index_mapping[first(path)]
        v = index_mapping[last(path)]
        path_weight = total_weight(osmg, path)
        if add_edge!(graph, (u, v))
            key = 0
            weights[u, v] = path_weight
            edge_count[u,v] = 1
        else # parallel edge
            key = edge_count[u,v]
            edge_count[u,v] += 1
            weights[u, v] = min(path_weight, weights[u, v])
        end
        edges[u,v,key] = path
    end

    edge_to_way = Dict{NTuple{3,U}, Vector{T}}()
    for (edge, path) in edges
        edge_to_way[edge] = ways_in_path(osmg, path)
    end

    return SimplifiedOSMGraph(
                osmg,
                node_coordinates,
                node_to_index,
                index_to_node,
                edge_to_way,
                graph,
                edges,
                weights,
                nothing
            )
end


text:
using CSV
using DataFrames
using DelimitedFiles
using Glob
using Plots; pgfplotsx()

# Create plots
function plot_gmres_ilu0()

	# Plot labels
	labels = ["Parallel ILU0", "Serial ILU0", "No precond."]

   	# Plot linestyles
   	linestyles = [:solid, :dashdot, :dash, :dot]

	# Benchmark folders
	folders = ["1e6_1000-samples", "1e7_100-samples", "atmosdd_100-samples"]

	# Solver

	# Iterate over matrix size number of samples folders
	for (ifolder, folder) in enumerate(folders)

		# Plot 1
		# Left: Solver runtimes none, ilu0, parilu0,
		# Right: Speedup parilu0 over ilu0

		# Get folders with csv files
		csvfolder_ser = joinpath(folder, "ser",
			readdir(joinpath(folder, "ser"))[1], "csv")
		csvfolder_par = joinpath(folder, "par",
			readdir(joinpath(folder, "par"))[1], "csv")

		# Create plots
		plt1 = plot(ylabel="Runtime in seconds", xlabel="Number of threads",
					legend=:topright, framestyle=:box)
		plt2 = plot(ylabel="Parallel speedup", xlabel="Number of threads",
					legend=:none, framestyle=:box)

		# Runtimes solver csv files
		csvfile_ser = joinpath(csvfolder_ser,
			"runtimes-solver_par_csr_dgks_ilu0.csv")
		csvfile_par = joinpath(csvfolder_par,
			"runtimes-solver_par_csr_dgks_parilu0.csv")
		csvfiles = [csvfile_par, csvfile_ser]
		if folder == "atmosdd_100-samples"
			csvfile_none = joinpath(csvfolder_ser,
				"runtimes-solver_par_csr_dgks_none.csv")
			csvfiles = [csvfile_par, csvfile_ser, csvfile_none]
		end

		# Add runtimes from all files to the first plot
		for (ifile, file) in enumerate(csvfiles)
		    df = DataFrame!(CSV.File(file, header=true))
		    x = convert(Vector, df[!, 1])
		    y = convert(Vector, df[!, 5])
			plot!(plt1, x, y, linestyles=linestyles[ifile], label=labels[ifile])
		end

		# Compute speedups and add to second plot
		df_par = DataFrame!(CSV.File(csvfile_par, header=true))
		df_ser = DataFrame!(CSV.File(csvfile_ser, header=true))
		y_par = convert(Vector, df_par[!, 5])
		y_ser = convert(Vector, df_ser[!, 5])
		x = convert(Vector, df_par[!, 1])
		y = y_ser ./ y_par
		plot!(plt2, x, y, linestyles=linestyles[1], label=labels[1])

		# Plot runtimes and speedups side by side
		plot(plt1, plt2, layout=(1,2), size=(530,530/2))
		savefig(string("gmres_ilu0_", folder, "_solver.png"))
		savefig(string("gmres_ilu0_", folder, "_solver.tex"))

	end # for

	# Precon

	# Iterate over matrix size number of samples folders
	for (ifolder, folder) in enumerate(folders)

		# Plot 1
		# Left: Solver runtimes none, ilu0, parilu0,
		# Right: Speedup parilu0 over ilu0

		# Get folders with csv files
		csvfolder_ser = joinpath(folder, "ser",
			readdir(joinpath(folder, "ser"))[1], "csv")
		csvfolder_par = joinpath(folder, "par",
			readdir(joinpath(folder, "par"))[1], "csv")

		# Create plots
		plt1 = plot(ylabel="Runtime in seconds", xlabel="Number of threads",
					legend=:topright, framestyle=:box)
		plt2 = plot(ylabel="Parallel speedup", xlabel="Number of threads",
					legend=:none, framestyle=:box)

		# Runtimes solver csv files
		csvfile_ser = joinpath(csvfolder_ser,
			"runtimes-precon_par_csr_dgks_ilu0.csv")
		csvfile_par = joinpath(csvfolder_par,
			"runtimes-precon_par_csr_dgks_parilu0.csv")
		csvfiles = [csvfile_par, csvfile_ser]

		# Add runtimes from all files to the first plot
		for (ifile, file) in enumerate(csvfiles)
		    df = DataFrame!(CSV.File(file, header=true))
		    x = convert(Vector, df[!, 1])
		    y = convert(Vector, df[!, 5])
			plot!(plt1, x, y, linestyles=linestyles[ifile], label=labels[ifile])
		end

		# Compute speedups and add to second plot
		# Compute speedups and add to second plot
		df_par = DataFrame!(CSV.File(csvfile_par, header=true))
		df_ser = DataFrame!(CSV.File(csvfile_ser, header=true))
		y_par = convert(Vector, df_par[!, 5])
		y_ser = convert(Vector, df_ser[!, 5])
		x = convert(Vector, df_par[!, 1])
		y = y_ser ./ y_par
		plot!(plt2, x, y, linestyles=linestyles[1], label=labels[1])

		# Plot runtimes and speedups side by side
		plot(plt1, plt2, layout=(1,2), size=(530,530/2))
		savefig(string("gmres_ilu0_", folder, "_precon.png"))
		savefig(string("gmres_ilu0_", folder, "_precon.tex"))

	end # for

end # function

max_stars_repo_name: jkrch/JuliaParGMRES-Plots

max_issues_repo_name: jkrch/JuliaParGMRES-Plots

max_forks_repo_name: jkrch/JuliaParGMRES-Plots


text:
using Random: BigFloat
using SemiseparableMatrices, BandedMatrices, Test, LinearAlgebra, MatrixFactorizations, LazyArrays, ArrayLayouts, Random
import BandedMatrices: _BandedMatrix, _banded_qr!
import SemiseparableMatrices: bandpart, fillpart, AlmostBandedLayout, VcatAlmostBandedLayout, resizedata!, _almostbanded_qr!
import MatrixFactorizations: QRPackedQ
import LazyArrays: CachedArray

Random.seed!(0)

@testset "AlmostBandedMatrix" begin
    @testset "Constructors" begin
        A = AlmostBandedMatrix{Float64}(undef, (10,11), (2,1), 2)
        A[1,1] = 2
        @test A[1,1] == 2
        A[4,1] = 0
        @test A[4,1] == 0.0
        @test_throws BandError A[4,1] = 2
        @test_throws ErrorException A[1,3] = 5
        @test almostbandwidths(A) == (2,1)
        @test almostbandedrank(A) == 2
    end

    @testset "copyto!" begin
        n = 10
        A = Vcat(Ones(1,n), BandedMatrix((0 => -Ones(n-1), 1 => 1:(n-1)), (n-1,n)))
        @test MemoryLayout(A) == VcatAlmostBandedLayout()
        @test almostbandwidths(A) == (1,0)
        @test almostbandedrank(A) == 1

        dest = AlmostBandedMatrix{Float64}(undef, size(A), almostbandwidths(A), almostbandedrank(A))
        copyto!(dest, A)
        @test dest == A

        dest = AlmostBandedMatrix{Float64}(undef, (1,2), (1,1), 1)
        copyto!(dest, view(A,1:1,1:2))
        @test dest == A[1:1,1:2]

        @test AlmostBandedMatrix(A) == Matrix(A)

        V = view(A,1:5,2:5)
        @test MemoryLayout(V) == VcatAlmostBandedLayout()
        @test almostbandwidths(V) == (2,-1)
        @test almostbandedrank(V) == 1
        @test A[1:5,2:5] isa AlmostBandedMatrix

        @test AlmostBandedMatrix(V) == Matrix(V) == V == A[1:5,2:5]

        V = view(A,Base.OneTo(5),Base.OneTo(5))
        @test AlmostBandedMatrix(V) == V == Matrix(V)

        V = view(AlmostBandedMatrix(A),1:5,2:5)
        @test MemoryLayout(V) == AlmostBandedLayout()
        @test AlmostBandedMatrix(V) == V == Matrix(V)

        A = Vcat(randn(2,n), BandedMatrix((0 => -Ones(n-1), 1 => 1:(n-1), 2 => Ones(n-2)), (n-2,n)))
        @test MemoryLayout(A) isa VcatAlmostBandedLayout
        @test almostbandedrank(A) == 2
        @test almostbandwidths(A) == (2,0)
        @test AlmostBandedMatrix(A) == Matrix(A)

        A = Vcat(randn(1,n), randn(1,n), BandedMatrix((0 => -Ones(n-1), 1 => 1:(n-1), 2 => Ones(n-2)), (n-2,n)))
        @test MemoryLayout(A) isa VcatAlmostBandedLayout
        @test almostbandedrank(A) == 2
        @test almostbandwidths(A) == (2,0)
        @test AlmostBandedMatrix(A) == Matrix(A)
    end

    @testset "cache" begin
        n = 10
        A = Vcat(Ones(1,n), BandedMatrix((0 => -Ones(n-1), 1 => 1:(n-1)), (n-1,n)))
        C = cache(A); resizedata!(C,2,2);
        @test C.data[1:2,1:2] == [1.0 1.0; -1.0 1.0]
        C = cache(A);
        resizedata!(C,10,10);
        @test C.data[Base.OneTo.(C.datasize)...] == A[1:10,1:10]
        @test C isa LazyArrays.CachedMatrix{Float64,<:AlmostBandedMatrix}

        l,u = almostbandwidths(A)
        r = almostbandedrank(A)
        data = AlmostBandedMatrix{Float64}(undef,(2l+u+1,0),(l,l+u),r) # pad super
        C = CachedArray(data,A);  resizedata!(C,1,2);
        @test C.data[1,2] == 1.0
        resizedata!(C,5,3);
        @test C.data[1:5,1:3] == A[1:5,1:3]

        A = Vcat(randn(2,n), BandedMatrix((0 => -Ones(n-1), 1 => 1:(n-1), 2 => Ones(n-2)), (n-2,n)))
        C = cache(A);
        resizedata!(C,10,10);
        @test C.data[Base.OneTo.(C.datasize)...] == A[1:10,1:10]
        @test C isa LazyArrays.CachedMatrix{Float64,<:AlmostBandedMatrix}

        @testset "BigFloat" begin
            A = Vcat(randn(2,n), BandedMatrix{BigFloat}((0 => -Ones(n-1), 1 => 1:(n-1), 2 => Ones(n-2)), (n-2,n)))
            C = cache(A);
            resizedata!(C,10,10);
            @test C.data[Base.OneTo.(C.datasize)...] == A[1:10,1:10]
            @test C isa LazyArrays.CachedMatrix{BigFloat,<:AlmostBandedMatrix}
        end
    end

    @testset "Slices" begin
        n = 80
        A = AlmostBandedMatrix(BandedMatrix(fill(2.0,n,n),(1,1)), LowRankMatrix(fill(3.0,n), ones(1,n)))
        V = view(A,1:3,1:3)
        @test MemoryLayout(A) == MemoryLayout(V) == AlmostBandedLayout()
        @test AlmostBandedMatrix(V) == A[1:3,1:3] == V
        @test A[1:3,1:3] isa AlmostBandedMatrix
    end

    @testset "colsupport" begin
        A = AlmostBandedMatrix{Float64}(undef, (10,11), (2,1), 2)
        @test colsupport(A,3) ≡ Base.OneTo(5)
        @test rowsupport(A,5) ≡ 3:11
    end

    @testset "Triangular" begin
        n = 80
        A = AlmostBandedMatrix(BandedMatrix(fill(2.0,n,n),(1,1)), LowRankMatrix(fill(3.0,n), ones(1,n)))
        b = randn(n)
        @test MemoryLayout(UpperTriangular(A)) == TriangularLayout{'U','N',AlmostBandedLayout}()
        @test UpperTriangular(Matrix(A)) \ b ≈ UpperTriangular(A) \ b
        @test UnitUpperTriangular(Matrix(A)) \ b ≈ UnitUpperTriangular(A) \ b
        @test LowerTriangular(Matrix(A)) \ b ≈ LowerTriangular(A) \ b
        @test UnitLowerTriangular(Matrix(A)) \ b ≈ UnitLowerTriangular(A) \ b
    end

    @testset "QR" begin
        n = 80
        A = AlmostBandedMatrix(BandedMatrix(fill(2.0,n,n),(1,1)), LowRankMatrix(fill(3.0,n), ones(1,n)))
        A[band(0)] .+= 1:n
        Ã = deepcopy(A)
        B,L = bandpart(A),fillpart(A)
        Ā = AlmostBandedMatrix(A,(1,2))
        @test A == Ā == B + triu(Matrix(L),2)
        F = qr(A)
        @test F.Q isa LinearAlgebra.QRPackedQ{Float64,<:BandedMatrix}
        @test F.R isa UpperTriangular{Float64,<:SubArray{Float64,2,<:AlmostBandedMatrix}}
        @test F.Q' * A ≈ F.R
        @test A == Ã

        Ā = AlmostBandedMatrix(A,(1,2))
        τ = Vector{Float64}(undef,n)
        @inferred(SemiseparableMatrices._almostbanded_qr!(Ā,τ))

        b = randn(n)
        @test A \ b ≈ Matrix(A) \ b
        @test all(A \ b .=== F \ b .=== F.R \ (F.Q'*b))
        Q̃ = QRPackedQ(F.factors,F.τ)
        @test Q̃ ≈ F.Q
        @test lmul!(Q̃,copy(b)) ≈ lmul!(F.Q,copy(b)) ≈ Matrix(F.Q)*b
        @test lmul!(Q̃',copy(b)) ≈ lmul!(F.Q',copy(b)) ≈ Matrix(F.Q)'*b

        A = Vcat(randn(2,n), BandedMatrix((0 => -Ones(n-1), 1 => 1:(n-1), 2 => Ones(n-2)), (n-2,n)))
        @test MemoryLayout(A) isa VcatAlmostBandedLayout
        @test qr(A) isa MatrixFactorizations.QR{Float64,<:AlmostBandedMatrix}
        @test qr(A) \ b ≈ Matrix(A) \ b

        A = Vcat(randn(1,n), randn(1,n), BandedMatrix((0 => -Ones(n-1), 1 => 1:(n-1), 2 => Ones(n-2)), (n-2,n)))
        @test MemoryLayout(A) isa VcatAlmostBandedLayout
    end

    @testset "one-col qr" begin
        A = AlmostBandedMatrix{Float64}(undef,(2,2),(1,1),1)
        A.bands.data .= randn.()
        A.fill.args[1] .= randn.()
        A.fill.args[2] .= randn.()
        Ã = deepcopy(A)
        _almostbanded_qr!(A,[1.0],1)
        @test qr(Ã[:,1]).Q' * Ã ≈ UpperTriangular(A)
    end

    @testset "almost banded degenerate band" begin
        A = Vcat(randn(2,6), BandedMatrix(2 => 1:5)[1:4,:])
        @test almostbandwidths(A) == (1,0)
        @test colsupport(A,1) == Base.OneTo(2)
        b = randn(6)
        @test A\b ≈ Matrix(A)\b
    end
end


text:
using DataFrames

"""Print in the format we used in the project group"""
function write_peaks(path::Union{AbstractString,IO}, df::DataFrames.DataFrame)
    d = copy(df)
    if !(:name in names(d))
        d[:name] = [Symbol(@sprintf("P%d", i)) for i=1:size(d, 1)]
        d = d[circshift(names(d), 1)]
    end
    # d[:t] .*= 100
    sort!(d, cols=order(:r))
    if isa(path, IO)
        printtable(path, d, separator='\t', header=true)
    else
        writetable(path, d, separator='\t', header=true)
    end
end

"""Read a peak list and return a DataFrame"""
function read_peaks(path::AbstractString)
    readtable(path, separator='\t')
end

export read_peaks, write_peaks


text:
struct LLVMNativeCode    # thin wrapper
    p::Ptr{Cvoid}
end

function xlinfo(f, tt)
    # get the method instance
    world = typemax(UInt)
    g = (args...) -> Cassette.overdub(ctx, f, args...)
    meth = which(g, tt)
    sig_tt = Tuple{typeof(g), tt.parameters...}
    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig_tt, meth.sig)::Core.SimpleVector

    if VERSION >= v"1.2.0-DEV.320"
        meth = Base.func_for_method_checked(meth, ti, env)
    else
        meth = Base.func_for_method_checked(meth, ti)
    end

    return ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                 (Any, Any, Any, UInt), meth, ti, env, world)
end

"""
Returns an LLVMNativeCode object for the function call `f` with TupleTypes `tt`.
"""
function raise_exception(insblock::BasicBlock, ex::Value)
end

# const jlctx = Ref{LLVM.Context}()

# function __init__()
#     jlctx[] = LLVM.Context(convert(LLVM.API.LLVMContextRef,
#                                    cglobal(:jl_LLVMContext, Nothing)))
# end

"""
    irgen(func, tt;
          optimize = true,
          optimize_llvm = true,
          fix_globals = true,
          overdub = true,
          module_setup = (m) -> nothing)

Generates Julia IR targeted for static compilation.
`ccall` and `cglobal` uses have pointer references changed to symbols
meant to be linked with libjulia and other libraries.

`optimize` controls Julia-side optimization. `optimize_llvm` controls
optimization on the LLVM side.

If `overdub == true` (the default), Cassette is used to swap out
`ccall`s with a tuple of library and symbol.

`module_setup` is an optional function to control setup of modules. It takes an LLVM
module as input.
"""
function irgen(@nospecialize(func), @nospecialize(tt);
               optimize = true,
               optimize_llvm = true,
               fix_globals = true,
               overdub = true,
               module_setup = (m) -> nothing)
    # get the method instance
    isa(func, Core.Builtin) && error("function is not a generic function")
    world = typemax(UInt)
    gfunc = overdub ? (args...) -> Cassette.overdub(ctx, func, args...) : func
    meth = which(gfunc, tt)
    sig_tt = Tuple{typeof(gfunc), tt.parameters...}
    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig_tt, meth.sig)::Core.SimpleVector

    if VERSION >= v"1.2.0-DEV.320"
        meth = Base.func_for_method_checked(meth, ti, env)
    else
        meth = Base.func_for_method_checked(meth, ti)
    end

    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                  (Any, Any, Any, UInt), meth, ti, env, world)

    current_method = nothing
    last_method_instance = nothing
    call_stack = Vector{Core.MethodInstance}()
    global method_map = Dict{String,Core.MethodInstance}()
    global dependencies = MultiDict{Core.MethodInstance,LLVM.Function}()
    # set-up the compiler interface
    function hook_module_setup(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        module_setup(LLVM.Module(ref))
    end
    function hook_raise_exception(insblock::Ptr{Cvoid}, ex::Ptr{Cvoid})
        insblock = convert(LLVM.API.LLVMValueRef, insblock)
        ex = convert(LLVM.API.LLVMValueRef, ex)
        raise_exception(BasicBlock(insblock), Value(ex))
    end
    function postprocess(ir)
        # get rid of jfptr wrappers
        for llvmf in functions(ir)
            startswith(LLVM.name(llvmf), "jfptr_") && unsafe_delete!(ir, llvmf)
        end

        return
    end
    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        global ir = LLVM.Module(ref)
        postprocess(ir)

        # find the function that this module defines
        llvmfs = filter(llvmf -> !isdeclaration(llvmf) &&
                                 linkage(llvmf) == LLVM.API.LLVMExternalLinkage,
                        collect(functions(ir)))
        llvmf = nothing
        if length(llvmfs) == 1
            llvmf = first(llvmfs)
        elseif length(llvmfs) > 1
            llvmfs = filter!(llvmf -> startswith(LLVM.name(llvmf), "julia_"), llvmfs)
            if length(llvmfs) == 1
                llvmf = first(llvmfs)
            end
        end
        insert!(dependencies, last_method_instance, llvmf)
        method_map[name(llvmf)] = current_method
    end
    function hook_emit_function(method_instance, code, world)
        push!(call_stack, method_instance)
    end
    function hook_emitted_function(method, code, world)
        current_method = method
        last_method_instance = pop!(call_stack)
        # @show code
        # dump(method, maxdepth=2)
        # global mymeth = method
    end

    params = Base.CodegenParams(cached=false,
                                track_allocations=false,
                                code_coverage=false,
                                static_alloc=false,
                                prefer_specsig=true,
                                module_setup=hook_module_setup,
                                module_activation=hook_module_activation,
                                raise_exception=hook_raise_exception,
                                emit_function=hook_emit_function,
                                emitted_function=hook_emitted_function,
                                )

    # get the code
    mod = let
        ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                    (Any, UInt, Bool, Bool, Base.CodegenParams),
                    linfo, world, #=wrapper=#false, #=optimize=#false, params)
        if ref == C_NULL
            #  error(jlctx[], "the Julia compiler could not generate LLVM IR")
        end

        llvmf = LLVM.Function(ref)
        LLVM.parent(llvmf)
    end

    # the main module should contain a single jfptr_ function definition,
    # e.g. jlcall_kernel_vadd_62977

    # definitions = filter(f->!isdeclaration(f), functions(mod))
    definitions = Iterators.filter(f->!isdeclaration(f), collect(functions(mod)))
    # definitions = collect(functions(mod))
    wrapper = let
        fs = collect(Iterators.filter(f->startswith(LLVM.name(f), "jfptr_"), definitions))
        @assert length(fs) == 1
        fs[1]
    end

    # the jlcall wrapper function should point us to the actual entry-point,
    # e.g. julia_kernel_vadd_62984
    entry_tag = let
        m = match(r"jfptr_(.+)_\d+", LLVM.name(wrapper))
        @assert m != nothing
        m.captures[1]
    end
    unsafe_delete!(mod, wrapper)
    entry = let
        re = Regex("julia_$(entry_tag)_\\d+")
        llvmcall_re = Regex("julia_$(entry_tag)_\\d+u\\d+")
        fs = collect(Iterators.filter(f->occursin(re, LLVM.name(f)) &&
                               !occursin(llvmcall_re, LLVM.name(f)), definitions))
        if length(fs) != 1
            compiler_error(func, tt, cap, "could not find single entry-point";
                           entry=>entry_tag, available=>[LLVM.name.(definitions)])
        end
        fs[1]
    end

    LLVM.name!(entry, string(nameof(func)))

    # link in dependent modules
    cache = Dict{String,String}()
    for called_method_instance in keys(dependencies)
        llvmfs = dependencies[called_method_instance]

        # link the first module
        llvmf = popfirst!(llvmfs)
        llvmfn = LLVM.name(llvmf)
        link!(mod, LLVM.parent(llvmf))
        # process subsequent duplicate modules
        for dup_llvmf in llvmfs
            if Base.JLOptions().debug_level >= 2
                # link them too, to ensure accurate backtrace reconstruction
                link!(mod, LLVM.parent(dup_llvmf))
            else
                # don't link them, but note the called function name in a cache
                dup_llvmfn = LLVM.name(dup_llvmf)
                cache[dup_llvmfn] = llvmfn
            end
        end
    end
    # resolve function declarations with cached entries
    for llvmf in filter(isdeclaration, collect(functions(mod)))
        llvmfn = LLVM.name(llvmf)
        if haskey(cache, llvmfn)
            def_llvmfn = cache[llvmfn]
            replace_uses!(llvmf, functions(mod)[def_llvmfn])
            unsafe_delete!(LLVM.parent(llvmf), llvmf)
        end
    end
    # rename functions to something easier to decipher
    # especially helps with overdubbed functions
    for (fname, mi) in method_map
        id = split(fname, "_")[end]
        basename = mi.def.name
        args = join(collect(mi.specTypes.parameters)[2:end], "_")
        if basename == :overdub  # special handling for Cassette
            basename = string(mi.specTypes.parameters[3])
            basename = replace(basename, r"^typeof\(" => "")
            basename = replace(basename, r"\)$" => "")
            args = join(collect(mi.specTypes.parameters)[4:end], "_")
        end
        newname = join([basename, args, id], "_")
        if haskey(functions(mod), fname)
            name!(functions(mod)[fname], newname)
        end
    end

    d = find_ccalls(gfunc, tt)
    fix_ccalls!(mod, d)
    if fix_globals
        fix_globals!(mod)
    end
    if optimize_llvm
        optimize!(mod)
    end
    return mod
end

"""
    optimize!(mod::LLVM.Module)

Optimize the LLVM module `mod`. Crude for now.
Returns nothing.
"""
function optimize!(mod::LLVM.Module)
    for llvmf in functions(mod)
        startswith(LLVM.name(llvmf), "jfptr_") && unsafe_delete!(mod, llvmf)
        startswith(LLVM.name(llvmf), "julia_") && LLVM.linkage!(llvmf, LLVM.API.LLVMExternalLinkage)
    end
    # triple = "wasm32-unknown-unknown-wasm"
    # triple!(mod, triple)
    # datalayout!(mod, "e-m:e-p:32:32-i64:64-n32:64-S128")
    # LLVM.API.@apicall(:LLVMInitializeWebAssemblyTarget, Cvoid, ())
    # LLVM.API.@apicall(:LLVMInitializeWebAssemblyTargetMC, Cvoid, ())
    # LLVM.API.@apicall(:LLVMInitializeWebAssemblyTargetInfo, Cvoid, ())
    triple = "i686-pc-linux-gnu"
    tm = TargetMachine(Target(triple), triple)

    ModulePassManager() do pm
        # add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        ccall(:jl_add_optimization_passes, Cvoid,
              (LLVM.API.LLVMPassManagerRef, Cint, Cint),
              LLVM.ref(pm), Base.JLOptions().opt_level, 1)

        dead_arg_elimination!(pm)
        global_optimizer!(pm)
        global_dce!(pm)
        strip_dead_prototypes!(pm)

        run!(pm, mod)
    end
    mod
end

function write_object(mod::LLVM.Module, path)
    host_triple = triple()
    host_t = Target(host_triple)
    TargetMachine(host_t, host_triple, "", "", LLVM.API.LLVMCodeGenLevelDefault, LLVM.API.LLVMRelocPIC) do tm
        emit(tm, mod, LLVM.API.LLVMObjectFile, path)
    end
end


text:
using JLBoost
using Test

@testset "JLBoost.jl" begin

end


text:
import LinearAlgebra, ProjectiveVectors
function setup_prediction_test()
    F = SPSystem(equations(katsura(6)))
    x = ProjectiveVectors.PVector(rand(Complex{Float64}, 7))
    xnext = copy(x)
    t = rand()
    H = HomotopyWithCache(StraightLineHomotopy(F, F), x, t)
    J = jacobian(H, x, t)
    Jac = HC.Jacobian(J)
    ẋ = Vector{ComplexF64}(undef, 7)
    H, x, xnext, t, ẋ, Jac
end

function setup_overdetermined_prediction_test()
    @polyvar x y
    F = SPSystem([x^2+y, y^2-3x*y, y+x+3, y+2x-5])
    x = ProjectiveVectors.PVector(rand(Complex{Float64}, 2))
    xnext = copy(x)
    t = rand()
    H = HomotopyWithCache(StraightLineHomotopy(F, F), x, t)
    Jac = HC.Jacobian(HC.jacobian(H, x, t))
    ẋ = Vector{ComplexF64}(undef, 2)
    H, x, xnext, t, ẋ, Jac
end

function test_predictor(PREDICTOR, CACHE)
    H, x, xnext, t, ẋ, Jac = setup_prediction_test()

    predictor = PREDICTOR()
    @test predictor isa PREDICTOR
    predictor_cache = cache(predictor, H, x, ẋ, t)
    setup!(predictor_cache, H, x, ẋ, t, Jac)
    @test predictor_cache isa CACHE
    # check that this doesn't throw
    update!(predictor_cache, H, x, ẋ, t, Jac)
    @test_nowarn predict!(xnext, predictor, predictor_cache, H, x, t, 0.05, ẋ, Jac)

    H, x, xnext, t, ẋ, Jac = setup_overdetermined_prediction_test()
    predictor = PREDICTOR()
    @test predictor isa PREDICTOR
    predictor_cache = cache(predictor, H, x, ẋ, t)
    setup!(predictor_cache, H, x, ẋ, t, Jac)
    @test predictor_cache isa CACHE
    # check that this doesn't throw
    update!(predictor_cache, H, x, ẋ, t, Jac)
    @test_nowarn predict!(xnext, predictor, predictor_cache, H, x, t, 0.05, ẋ, Jac)
end

@testset "Predictors" begin

    @testset "NullPredictor" begin
        test_predictor(NullPredictor, HomotopyContinuation.NullPredictorCache)
    end

    @testset "Euler" begin
        test_predictor(Euler, HomotopyContinuation.EulerCache)
    end

    @testset "Heun" begin
        test_predictor(Heun, HomotopyContinuation.HeunCache)
    end

    @testset "Midpoint" begin
        test_predictor(Midpoint, HomotopyContinuation.MidpointCache)
    end

    @testset "Ralston" begin
        test_predictor(Ralston, HomotopyContinuation.RalstonCache)
    end

    @testset "RK3" begin
        test_predictor(RK3, HomotopyContinuation.RK3Cache)
    end

    @testset "Pade21" begin
        test_predictor(Pade21, HomotopyContinuation.Pade21Cache)
    end

    @testset "RK4" begin
        test_predictor(RK4, HomotopyContinuation.RK4Cache)
    end
end


text:
typedef struct Pixel_t{
  int r, g, b;
} * Pixel;

typedef struct Screen_t{
  Pixel[][] pixels;
  boolean active;
} * Screen;

void initPixel(Pixel[][] pixels){
  int i = 0;
  while(i < pixels.length){
    int j = 0;
    while(j < pixels[i].length){
      Pixel p = new Pixel_t;
      p->r = i;
      p->g = j;
      p->b = i*j;
      pixels[i][j] = p;
      j++;
    }
    i++;
  }
}

Screen createScreen(){
  Screen screen = new Screen_t;
  screen->pixels = new Pixel[640][480];
  initPixel(screen->pixels);
  screen->active = true;
  return screen;
}

int main(){
  Screen screen = createScreen();
  printInt(screen->pixels[4][2]->r);
  printInt(screen->pixels[4][2]->g);
  printInt(screen->pixels[4][2]->b);
  // prints 4
  //        2
  //        8
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
# This is an implementation of second order physics constrained learning

export pcl_sparse_solve, pcl_square_sum, pcl_hessian, pcl_linear_op, pcl_compress

@doc raw"""
    pcl_sparse_solve(indices::Array{Int64, 2},
        vals::Array{Float64, 1},
        u::Array{Float64, 1},
        hessian_u::Array{Float64, 2}, grad_u::Array{Float64, 1})

Hessian backpropagation for

$$Au = f$$

`A` is given by a sparse matrix (the nonzero entries are a vector $a$). This function
back-propagates the Hessian w.r.t. $u$ to the Hessian w.r.t. $A$.
"""
function pcl_sparse_solve(indices::Array{Int64, 2},
    vals::Array{Float64, 1},
    u::Array{Float64, 1},
    hessian_u::Array{Float64, 2}, grad_u::Array{Float64, 1})
    A = sparse(indices[:,1], indices[:,2], vals)
    indof = length(vals)
    outdof = length(u)
    invA = inv(Array(A))
    J = zeros(indof, outdof)
    for k = 1:indof
        i = indices[k,1]
        j = indices[k,2]
        J[k,:] = -invA[:,i] * u[j]
    end
    H = zeros(indof, indof)
    x = A'\grad_u
    for l = 1:indof
        for r = 1:indof
            li, lj = indices[l, :]
            ri, rj = indices[r, :]
            H[l, r] = -(x[li] * J[r, lj] + x[ri] * J[l, rj])
        end
    end
    return J*hessian_u*J' + H
end

@doc raw"""
    pcl_square_sum(scaling::Float64 = 1.0)

Returns the Hessian matrix for
$$\text{scaling} * \|u - u_o\|_2^2$$

`pcl_square_sum` is a reduction op, which can provides the highest level Hessian matrix.
"""
pcl_square_sum(n::Int64; scaling::Float64 = 1.0) = 2Array(I(n)) * scaling

@doc raw"""
    pcl_linear_op(J::Array{Float64, 2}, W::Array{Float64,2})

For a linear operator

$$y = J^Tx + b$$

The PCL backpropagation is given by
$$H = JWJ^T$$
"""
function pcl_linear_op(J::Array{Float64, 2}, W::Array{Float64,2})
    J * W * J'
end

@doc raw"""
    pcl_hessian(y::PyObject, x::PyObject)

Returns the Hessian tensor for the operator

$$y = F(x)$$

$x\in \mathbb{R}^m$ and $y\in \mathbb{R}^n$ are 1D vectors. The function returns a tuple

- `H`: $m\times m$ Hessian matrix
- `W`: $n\times n$ Hessian matrix from downstream of the computational graph
- `dy`: length $n$ tensor; gradient from downstream of the computational graph
"""
function pcl_hessian(y::PyObject, x::PyObject, loss::PyObject)
    @assert length(size(x))==length(size(y))==1
    J = jacobian(y, x)
    if isnothing(J)
        error("`jacobian(y, x) = nothing. Please check the dependency of `y` on `x`")
    end
    W = placeholder(Float64, shape = (length(y), length(y)))
    dy = independent(gradients(loss, y))
    H = J' * W * J + hessian(dot(dy, y), x)
    return H, W
end

@doc raw"""
    pcl_compress(indices::Array{Int64, 2})

Computes the Jacobian matrix for `compress`. Assume that `compress` does the following transformation:
```
indices, values (v) --> new_indices, new_values (u)
```
This function computes
$$J_{ij} = \frac{\partial u_j}{\partial v_i}$$
"""
function pcl_compress(indices::Array{Int64, 2})
    N = size(indices,1)
    v = zeros(N)
    nout = zeros(Int32,1)
    indices = indices'[:] .- 1
    out = @eval ccall((:pcl_SparseCompressor, $(ADCME.LIBADCME)), Ptr{Cdouble},
        (Ptr{Clonglong}, Ptr{Cdouble}, Cint, Ptr{Cint}),
         $indices, $v, Int32($N), $nout)
    J = unsafe_wrap(Array, out, (N,Int64(nout[1])); own = true)
    return J
end


text:
f = open("hightemp.txt")
for line in readlines(f)
  print(replace(line,"\t", " "))
end


text:
using DataGraphs
using Dictionaries
using Graphs

g = grid((4,))
dg = DataGraph{String,Symbol}(g)
@show !isassigned(dg, Edge(1, 2))
@show !isassigned(dg, 1 => 2)
@show !isassigned(dg, Edge(1 => 2))
@show !isassigned(dg, 1 => 3)
@show !isassigned(dg, 1)
@show !isassigned(dg, 2)
@show !isassigned(dg, 3)
@show !isassigned(dg, 4)

@show has_edge(dg, 1, 2)
@show has_edge(dg, 1 => 2)
@show !has_edge(dg, 1, 3)
@show !has_edge(dg, 1 => 3)
@show has_vertex(dg, 1)
@show has_vertex(dg, 4)
@show !has_vertex(dg, 0)
@show !has_vertex(dg, 5)

dg[1] = "V1"
dg[2] = "V2"
dg[3] = "V3"
dg[4] = "V4"
@show isassigned(dg, 1)
@show dg[1] == "V1"
@show dg[2] == "V2"
@show dg[3] == "V3"
@show dg[4] == "V4"

dg[1 => 2] = :E12
dg[2 => 3] = :E23
dg[Edge(3, 4)] = :E34
#@show isassigned(dg, (1, 2))
@show isassigned(dg, Edge(2, 3))
@show isassigned(dg, 3 => 4)
@show dg[Edge(1, 2)] == :E12
@show dg[2 => 3] == :E23
@show dg[3 => 4] == :E34


text:
__precompile__()

module mbInterpolate
export interpolateLanczos, interpolateTrigonometric

using mbUtil

"""
    interpolateLanczos(V, a::Integer)

simple periodic (!) Lanczos interpolation
"""
function interpolateLanczos(V,a::Integer)
	return mbUtil.vectorize() do y
		sum = zero(V[1])
		N = length(V)
		for i in floor(Integer,y)+(-a+1:a)
			local tmp = y-i
			sum += sinc(tmp)*sinc(tmp/a)*V[mod(i-1,N) + 1]
		end
		return sum
	end
end
interpolateLanczos(V, a::Integer, x) = interpolateLanczos(V,a)(x)

"""
    interpolateTrigonometric(a₀, a, b)
Returns trigonometric polynomial.
Use with 2a,-2b and divide by 2m+1 to use with rfft coefficients.
"""
function interpolateTrigonometric(a₀, a, b)
	return mbUtil.vectorize() do x
		a₀ + reduce(0, 1:length(a)) do I,i
			tmp = i*x
			I + a[i]*cos(tmp) + b[i]*sin(tmp)
		end
	end
end
interpolateTrigonometric(a₀, a, b, x) = interpolateTrigonometric(a₀, a, b)(x)

end #module


text:
using Test
using Butterfly

max_stars_repo_name: zclaytor/butterpy

max_issues_repo_name: zclaytor/butterpy

max_forks_repo_name: zclaytor/butterpy


text:
using Compat
using CurveFit
using Compat.Test
using Compat.LinearAlgebra
using Polynomials
# write your own tests here
@test 1 == 1

if VERSION < v"0.7-"
    lspf(start, stop, len) = linspace(start, stop, len)
else
    lspf(start, stop, len) = range(start, stop=stop, length=len)
end

# Testing linear fit
x = lspf(1,10,10);
fun0(x) = 1.0 + 2.0*x
y = fun0.(x)
f = linear_fit(x, y)

@test f[1] ≈ 1.0 atol=1.0e-7
@test f[2] ≈ 2.0 atol=1.0e-7

f = curve_fit(LinearFit, x, y)
@test f(1.5) ≈ fun0(1.5) atol=1.0e-7

# power
fun1(x) = 1.0 + 2.0*log(x)
y = fun1.(x)
f = log_fit(x, y)
@test f[1] ≈ 1.0 atol=1.0e-7
@test f[2] ≈ 2.0 atol=1.0e-7
f = curve_fit(LogFit, x, y)
@test f(1.5) ≈ fun1(1.5) atol=1.0e-7

# log
fun2(x) = 2.0*x.^0.8
y = fun2(x)
f = power_fit(x, y)
@test f[1] ≈ 2.0 atol=1.0e-7
@test f[2] ≈ 0.8 atol=1.0e-7
f = curve_fit(PowerFit, x, y)
@test f(1.5) ≈ fun2(1.5) atol=1.0e-7

# Exp
fun3(x) = 2.0*exp.(0.8*x)
y = fun3(x)
f = exp_fit(x, y)
@test f[1] ≈ 2.0 atol=1.0e-7
@test f[2] ≈ 0.8 atol=1.0e-7
f = curve_fit(ExpFit, x, y)
@test f(1.5) ≈ fun3(1.5) atol=1.0e-7

# Poly
fun4(x) = 1.0 + 2.0*x + 3.0*x^2 + 0.5*x^3
y = fun4.(x)
f = poly_fit(x, y, 4)
@test f[1] ≈ 1.0 atol=1.0e-7
@test f[2] ≈ 2.0 atol=1.0e-7
@test f[3] ≈ 3.0 atol=1.0e-7
@test f[4] ≈ 0.5 atol=1.0e-7
@test f[5] ≈ 0.0 atol=1.0e-7

f = curve_fit(Poly, x, y, 4)
@test polyval(f, 1.5) ≈ fun4(1.5) atol=1.0e-7

# King's law
U = [lspf(1, 20, 20);]
A = 5.0
B = 1.5
n = 0.5
E = sqrt.(A .+ B*U.^n)
fun5(E) = ((E.^2 - A)/B).^(1 ./ n)

f = linear_king_fit(E, U)
@test f[1] ≈ A atol=1.0e-7
@test f[2] ≈ B atol=1.0e-7
f= curve_fit(LinearKingFit, E, U)
@test f(3.0) ≈ fun5(3.0) atol=1.0e-7

# Modified King's law
n = 0.42

E = sqrt.(A .+ B.*U.^n)
fun6(E) = ((E^2 - A)/B)^(1 / n)

f = king_fit(E, U)
@test f[1] ≈ A atol=1.0e-7
@test f[2] ≈ B atol=1.0e-7
@test f[3] ≈ n atol=1.0e-7
f= curve_fit(KingFit, E, U)
@test f(3.0) ≈ fun6(3.0) atol=1.0e-5

# Linear Rational fit
r =  RationalPoly([1.0, 0.0, -2.0], [1.0, 2.0, 3.0])
y = r(x)
f = linear_rational_fit(x, y, 2, 3)
@test f[1] ≈ 1.0 atol=1.0e-8
@test f[2] ≈ 0.0 atol=1.0e-8
@test f[3] ≈ -2.0 atol=1.0e-8
@test f[4] ≈ 2.0 atol=1.0e-8
@test f[5] ≈ 3.0 atol=1.0e-8
@test f[6] ≈ 0.0 atol=1.0e-8

# Nonlinear Rational fit
f = rational_fit(x, y, 2, 3)

@test f[1] ≈ 1.0 atol=1.0e-7
@test f[2] ≈ 0.0 atol=1.0e-7
@test f[3] ≈ -2.0 atol=1.0e-7
@test f[4] ≈ 2.0 atol=1.0e-7
@test f[5] ≈ 3.0 atol=1.0e-7
@test f[6] ≈ 0.0 atol=1.0e-7

f = curve_fit(RationalPoly, x, y, 2, 3)
@test f(1.5) ≈ r(1.5) atol=1.0e-8
@test f(4.5) ≈ r(4.5) atol=1.0e-8


text:
function __init__()

    # println("calling function __init__") # raised error when starting Julia REPL with "$sysimg_path.dll": "calling function __init__fatal: error thrown and no exception handler available."

    # remove MySql dependency
    #
    ## MySQL.MySQLHandle objects involve raw Ptr objects and need to be defined in `__init__`
    ## (runtime instead of compile time), since memory layout does not remain the same across
    ## process restarts (e.g. between compile time and runtime).
    #
    global const DB_INFO = Dict{String,String}(
        "database" => "chaipcr",
        "host"     => "127.0.0.1",
        "username" => "root",
        "password" => ""
    )
    global const DB_CONN_DICT = OrderedDict(map([
        ("default", DB_INFO["database"]),
        ("t1", "test_1ch"),
        ("t2", "test_2ch")
    ]) do db_tuple
        db_short, db_real = db_tuple
        db_short => try
            MySQL.mysql_connect(DB_INFO["host"], DB_INFO["username"], DB_INFO["password"], db_real)
        catch err
            if isa(err, MySQL.MySQLInternalError)
                warn("Database \"$db_real\" does not exist, returning MySQL.MySQLInternalError")
            else
                warn("Unknown error (other than database absence) occurred when attempt connection to \"$db_real\". ease report this to code owner.")
            end
            err
        end
    end) # do db_name

    # test
    if all(map(conn_outcome -> isa(conn_outcome, MySQL.MySQLHandle), values(DB_CONN_DICT)))
        println("start of test")
        @time global const TEST_RESULTS = test(;
            comprehensive=false,
            debug=true,
            verbose=true # `false` on PC, true on BBB to debug
       )
        println("end of test")
    else
        warn("Not all default or test databases were connected successfully. Test is not performed during initialization.")
    end

    return nothing
end

#

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
abstract type AbstractStorage end

struct BSONStorage <: AbstractStorage end
struct JLD2Storage <: AbstractStorage end
struct JLSOStorage <: AbstractStorage end
struct JSONStorage <: AbstractStorage end
struct YAMLStorage <: AbstractStorage end

const EXTENSIONS = Dict(
    "bson" => BSONStorage,
    "jld2" => JLD2Storage,
    "jlso" => JLSOStorage,
    "json" => JSONStorage,
    "yml" => YAMLStorage,
    "yaml" => YAMLStorage,
)

function get_storage(path, default)
    _, ext = splitext(path)
    key = lowercase(isempty(ext) ? default : ext[2:end])
    haskey(EXTENSIONS, key) || throw(ArgumentError("Unknown extension $ext"))
    storage = EXTENSIONS[key]
    isempty(ext) && (path *= ".$default")
    return storage, path
end

function resp_to_dict(resp)
    req = resp.request
    return Dict(
        "status" => resp.status,
        "headers" => resp.headers,
        "body" => repr_body(resp.body),
        "version" => string(resp.version),
        "request" => Dict(
            "method" => req.method,
            "headers" => req.headers,
            "target" => req.target,
            "body" => repr_body(req.body),
            "txcount" => req.txcount,
            "version" => string(req.version),
        ),
    )
end

function dict_to_resp(dict)
    resp = Response()
    resp.status = dict["status"]
    resp.headers = dicts_to_pairs(dict["headers"])
    resp.body = Vector{UInt8}(dict["body"])
    resp.version = VersionNumber(dict["version"])
    resp.request.method = dict["request"]["method"]
    resp.request.body = Vector{UInt8}(dict["request"]["body"])
    resp.request.headers = dicts_to_pairs(dict["request"]["headers"])
    resp.request.target = dict["request"]["target"]
    resp.request.txcount = dict["request"]["txcount"]
    resp.request.version = VersionNumber(dict["request"]["version"])
    return resp
end

repr_body(body) = isvalid(String, body) ? String(body) : resp.body
dicts_to_pairs(dicts) = map(d -> first(pairs(d)), dicts)
write_json(path, data) = open(io -> JSON.print(io, data, 2), path, "w")

load(::Type{BSONStorage}, path) = BSON.load(path)
store(::Type{BSONStorage}, path; responses, format) =
    BSON.bson(path; responses=responses, format=format)

function load(::Type{JLD2Storage}, path)
    JLD2.@load path responses format
    return Dict(:responses => responses, :format => format)
end

store(::Type{JLD2Storage}, path; responses, format) = JLD2.@save path responses format

load(::Type{JLSOStorage}, path) = JLSO.load(path)
store(::Type{JLSOStorage}, path; responses, format) =
    JLSO.save(path, :responses => responses, :format => format)

for (Storage, load_fun, store_fun) in (
    (JSONStorage, JSON.parsefile, write_json),
    (YAMLStorage, YAML.load_file, YAML.write_file),
)
    @eval begin
        function load(::Type{$Storage}, path)
            data = $load_fun(path)
            return Dict(
                :responses => map(dict_to_resp, data["responses"]),
                :format => VersionNumber(data["format"]),
            )
        end

        function store(::Type{$Storage}, path; responses, format)
            data = Dict(
                :responses => map(resp_to_dict, responses),
                :format => string(format),
            )
            $store_fun(path, data)
        end
    end
end


text:
"""
    brightness(c) -> Float64

对给定的颜色三通道`像素值`，求出其亮度并以浮点数的形式输出

# Notice
由于实验中记录的全息图都是灰度图像，其三个通道分量颜色都是一致的，故取某一颜色分量作为其亮度即可

# Arguments
- c: RGB类型的像素值
"""
brightness(c::AbstractRGB) = c.r |> Float64

"""
    brightness(::AbstractMatrix{T}) -> AbstractMatrix{Float64}

对给定的颜色三通道`图像`，求出其亮度并以`浮点数`的形式输出

# Notice
由于实验中记录的全息图都是灰度图像，其三个通道分量颜色都是一致的，故取某一颜色分量作为其亮度即可
"""
brightness(img::AbstractMatrix{<:AbstractRGB}) = brightness.(img)

"""
    load_holo(PATH::String, name::String; convert=true) -> AbstractMatrix

给定`图像`路径以及`图像`的名称(包括后缀)，实现对全息图的载入

# Arguments
- `PATH`: `图像`所在路径
- `name`: `图像`的名称，比如`tail.bmp`
- `convert`: 可选参数，默认值是`true`。当为`true`时，该函数返回`图像`的灰度值浮点数矩阵；当为`false`时，该函数直接返回图像
"""
function load_holo(PATH::String, name::String; convert = true)
    img = load(joinpath(PATH, name), view = true)

    return (convert == true) ? brightness(img) : img
end

function _normalize(x::AbstractMatrix)
    min_value, max_value = minimum(x), maximum(x)
    min_max = max_value - min_value
    return (x .- min_value) ./ min_max
end

function _normalize!(x::AbstractMatrix)
    min_value, max_value = minimum(x), maximum(x)
    min_max = max_value - min_value

    for i in eachindex(x)
        x[i] = (x[i] - min_value) / min_max
    end

    return x
end

function _normalize_threads(x::AbstractMatrix)
    min_value, max_value = minimum(x), maximum(x)
    min_max = max_value - min_value
    src = similar(x, size(x))

    @inbounds Threads.@threads for i in eachindex(x)
        src[i] = (x[i] - min_value) / min_max
    end

    return src
end

function _normalize_threads!(x::AbstractMatrix)
    min_value, max_value = minimum(x), maximum(x)
    min_max = max_value - min_value

    @inbounds Threads.@threads for i in eachindex(x)
        x[i] = (x[i] - min_value) / min_max
    end

    return x
end

"""
    normalize(x::AbstractMatrix; nthreads=false) -> AbstractMatrix

对矩阵`x`进行归一化，并返回归一化后的矩阵

# Arguments
- `x`: 矩阵
- `nthreads`: 决定是否开多线程来进行归一化，默认为否

# Notice
对于不大的矩阵（低于10_000 × 10_000），不建议开多线程，这是因为多线程会有额外的内存和时间开销
"""
normalize(x::AbstractMatrix; nthreads = false) =
    nthreads ? _normalize_threads(x) : _normalize(x)

"""
    normalize!(x::AbstractMatrix; nthreads=false) -> AbstractMatrix

对矩阵`x`进行`原地`归一化，即原地修改矩阵并返回修改后的`x`

# Arguments
- `x`: 矩阵
- `x`: 矩阵
- `nthreads`: 决定是否开多线程来进行归一化，默认为否

# Notice
对于不大的矩阵（低于10_000 × 10_000），不建议开多线程，这是因为多线程会有额外的内存和时间开销
"""
normalize!(x::AbstractMatrix; nthreads = false) =
    nthreads ? _normalize_threads!(x) : _normalize!(x)

"""
    color(x::AbstractMatrix{<:AbstractFloat}) -> AbstractMatrix
将浮点数类型矩阵转换为灰度图像

# Arguments
- x: 浮点数矩阵
"""
color(x::AbstractMatrix{<:AbstractFloat}) = Gray.(x)

"""
    reconst(holo::AbstractMatrix, P::FFTW.FFTWPlan, scale::Integer; shift=false) -> AbstractMatrix
对无透镜傅立叶全息图进行数值重建

# Arguments
- `holo`: 无透镜傅立叶全息图
- `P`: 其类型为FFTW.FFTWPlan，表示重建采用`fft`傅立叶变换
- `scale`: 由于零级像强度过大，故需要用`scale`拉升±1级像的强度
- `shift`: 是否旋转图像，默认为否（除非是想展示完整的重建像，否则不建议将`shift`设为`true`，因为其内存和时间开销比较大）
- `nthreads`: 是否开启多线程，默认为是（值得注意的是，该线程与傅立叶变换多线程的开启无关）
"""
function reconst(
    holo::AbstractMatrix,
    P::FFTW.FFTWPlan,
    scale::Integer;
    shift = false,
    nthreads = true,
)
    if shift
        fft_img =
            @pipe fftshift(P * fftshift(holo)) .|> abs |> normalize!(_, nthreads = nthreads)
    else
        fft_img = @pipe P * holo .|> abs |> normalize!(_, nthreads = nthreads)
    end
    fft_img .*= scale
    return clamp!(fft_img, 0, 1)
end


text:
function strictly_lessthan_singlevariable_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, CP.Domain{Int})
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        CP.Strictly{MOI.LessThan{Int}, Int},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, CP.Domain(Set([1, 2])))
    c2 = MOI.add_constraint(model, x2, MOI.EqualTo(2))

    c3 = MOI.add_constraint(
        model,
        1 * x1 - 1 * x2,
        CP.Strictly(MOI.LessThan(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 2
    end
end

function strictly_greaterthan_singlevariable_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, CP.Domain{Int})
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        CP.Strictly{MOI.GreaterThan{Int}, Int},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, CP.Domain(Set([1, 2])))
    c2 = MOI.add_constraint(model, x2, MOI.EqualTo(2))

    c3 = MOI.add_constraint(
        model,
        1 * x1 - 1 * x2,
        CP.Strictly(MOI.GreaterThan(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 2
    end
end

const strictlytests = Dict(
    "strictly_lessthan_singlevariable" =>
        strictly_lessthan_singlevariable_test,
    "strictly_greaterthan_singlevariable" =>
        strictly_greaterthan_singlevariable_test,
)

MOIT.@moitestset strictly


text:
"""
```
abstract AbstractScenario
```

Abstract supertype for all alternative scenarios.

### Required Fields

- `key::Symbol`
- `description::String`
- `vintage::String`
"""
abstract type AbstractScenario end

abstract type SingleScenario <: AbstractScenario end

type Scenario <: SingleScenario
    key::Symbol
    description::String
    target_names::Vector{Symbol}
    instrument_names::Vector{Symbol}
    targets::DataFrame
    instruments::DataFrame
    shock_scaling::Float64
    draw_states::Bool
    altpolicy::AltPolicy
    vintage::String
end

function Base.show(io::IO, scen::Scenario)
    @printf io "%-14s %s\n" "Key:" scen.key
    @printf io "%-14s %s\n" "Description:" scen.description
    @printf io "%-14s %s\n" "Targets:" scen.target_names
    @printf io "%-14s %s\n" "Instruments:" scen.instrument_names
    if scen.shock_scaling != 1.0
        @printf io "%-14s %s\n" "Shock Scaling:" scen.shock_scaling
    end
    if scen.draw_states
        @printf io "%-14s %s\n" "Draw States:" scen.draw_states
    end
    if scen.altpolicy.key != :historical
        @printf io "%-14s %s\n" "Alt Policy:" scen.altpolicy
    end
    @printf io "%-14s %s"   "Vintage:" scen.vintage
end

"""
```
Scenario(key, description, target_names, instrument_names, vintage;
    shock_scaling = 1.0, draw_states = false,
    altpolicy = AltPolicy(:historical, eqcond, solve))
```

Scenario constructor. If `instrument_names` is empty, then all model shocks will
be used.
"""
function Scenario(key::Symbol, description::String,
                  target_names::Vector{Symbol},
                  instrument_names::Vector{Symbol},
                  vintage::String;
                  shock_scaling::Float64 = 1.0,
                  draw_states::Bool = false,
                  altpolicy::AltPolicy = AltPolicy(:historical, eqcond, solve))
    Scenario(key, description, target_names, instrument_names,
             DataFrame(), DataFrame(), shock_scaling, draw_states, altpolicy, vintage)
end

n_targets(scen::Scenario) = length(scen.target_names)
n_instruments(scen::Scenario) = length(scen.instrument_names)
n_target_horizons(scen::Scenario) = size(scen.targets, 1)

function targets_to_data(m::AbstractModel, scen::Scenario)
    df = DataFrame()

    # Assign dates
    horizons   = n_target_horizons(scen)
    start_date = date_forecast_start(m)
    end_date   = DSGE.iterate_quarters(start_date, horizons - 1)
    df[:date]  = quarter_range(start_date, end_date)

    for var in keys(m.observables)
        df[var] = if var in scen.target_names
            scen.targets[var]
        else
            fill(NaN, horizons)
        end
    end
    return df
end

type SwitchingScenario <: SingleScenario
    key::Symbol
    description::String
    vintage::String
    original::SingleScenario
    default::SingleScenario
    probs_enter::Vector{Float64}
    probs_exit::Vector{Float64}

    function SwitchingScenario(key, description, vintage, original, default,
                               probs_enter, probs_exit)
        if length(probs_enter) != length(probs_exit)
            error("Lengths of probs_enter and probs_exit must be the same")
        end
        if !all(p -> 0.0 <= p <= 1.0, probs_enter)
            error("Elements of probs_enter must be between 0 and 1")
        end
        if !all(p -> 0.0 <= p <= 1.0, probs_exit)
            error("Elements of probs_exit must be between 0 and 1")
        end
        return new(key, description, vintage, original, default, probs_enter, probs_exit)
    end
end

function Base.show(io::IO, scen::SwitchingScenario)
    @printf io "%-12s %s\n" "Key:" scen.key
    @printf io "%-12s %s\n" "Description:" scen.description
    @printf io "%-12s %s\n" "Original:" scen.original.key
    @printf io "%-12s %s\n" "Default:" scen.default.key
    @printf io "%-12s %s"   "Vintage:" scen.vintage
end

"""
```
SwitchingScenario(key, original, default, probs_enter, probs_exit;
                  description = original.description,
                  vintage = original.vintage)
```

Constructs an instance of `SwitchingScenario` from two scenarios and
two vectors of entry/exit probabilities.
"""
function SwitchingScenario(key::Symbol, original::Scenario, default::Scenario,
                           probs_enter::Vector{Float64},
                           probs_exit::Vector{Float64};
                           description::String = original.description,
                           vintage::String = original.vintage)

    SwitchingScenario(key, description, vintage, original, default,
                      probs_enter, probs_exit)
end

"""
```
ScenarioAggregate(key, description, scenarios, sample, proportions, total_draws,
    replace, vintage)

ScenarioAggregate(key, description, scenarios, vintage)

ScenarioAggregate(key, description, scenarios, proportions, total_draws,
    replace, vintage)
```

Composite type for aggregated `AbstractScenario`s. Scenarios in
`scenarios::Vector{AbstractScenario}` are sampled into the aggregate
with the probability of the corresponding entry in `proportions` (whose values
must sum to 1). The field `replace` indicates whether to sample with
replacement.
"""
type ScenarioAggregate <: AbstractScenario
    key::Symbol
    description::String
    scenarios::Vector{AbstractScenario}
    sample::Bool # if false, scenario draws are kept in original proportions and `proportions` and `total_draws` are filled in after reading in draws
    proportions::Vector{Float64}
    total_draws::Int
    replace::Bool # whether to sample with replacement
    vintage::String

    function ScenarioAggregate(key, description, scenarios, sample,
                               proportions, total_draws, replace, vintage)
        if sample
            if length(scenarios) != length(proportions)
                error("Lengths of scenarios and proportions must be the same")
            end
            if !all(p -> 0.0 <= p <= 1.0, proportions)
                error("Elements of proportions must be between 0 and 1")
            end
            if sum(proportions) != 1.0
                error("Elements of proportions must sum to 1")
            end
        end
        return new(key, description, scenarios, sample,
                   proportions, total_draws, replace, vintage)
    end
end

# `sample = true` constructor
function ScenarioAggregate(key::Symbol, description::String, scenarios::Vector{AbstractScenario},
                           proportions::Vector{Float64}, total_draws::Int, replace::Bool,
                           vintage::String)
    sample = true
    return ScenarioAggregate(key, description, scenarios, sample, proportions, total_draws,
                             replace, vintage)
end

# `sample = false` constructor
function ScenarioAggregate(key::Symbol, description::String, scenarios::Vector{AbstractScenario},
                           vintage::String)
    sample = false
    proportions = Float64[]
    total_draws = -1
    replace = false
    return ScenarioAggregate(key, description, scenarios, sample, proportions, total_draws,
                             replace, vintage)
end

function Base.show(io::IO, agg::ScenarioAggregate)
    @printf io "%-24s %s\n" "Key:" agg.key
    @printf io "%-24s %s\n" "Description:" agg.description
    @printf io "%-24s %s\n" "Scenarios:" map(scen -> scen.key, agg.scenarios)
    @printf io "%-24s %s\n" "Sample:" agg.sample
    @printf io "%-24s %s\n" "Proportions:" agg.proportions
    @printf io "%-24s %s\n" "Total Draws:" agg.total_draws
    @printf io "%-24s %s\n" "Sample with Replacement:" agg.replace
    @printf io "%-24s %s"   "Vintage:" agg.vintage
end


text:
module FooTests

import Foo

end


text:
"""
$(TYPEDEF)

A nonlinear system of equations.

# Fields
$(FIELDS)

# Examples

```
@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β])
```
"""
struct NonlinearSystem <: AbstractSystem
    """Vector of equations defining the system."""
    eqs::Vector{Equation}
    """Unknown variables."""
    states::Vector{Variable}
    """Parameters."""
    ps::Vector{Variable}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems
    """
    systems::Vector{NonlinearSystem}
end

function NonlinearSystem(eqs, states, ps;
                         name = gensym(:NonlinearSystem),
                         systems = NonlinearSystem[])
    NonlinearSystem(eqs, convert.(Variable,states), convert.(Variable,ps), name, systems)
end

function calculate_jacobian(sys::NonlinearSystem)
    rhs = [eq.rhs for eq in sys.eqs]
    jac = expand_derivatives.(calculate_jacobian(rhs, [dv() for dv in states(sys)]))
    return jac
end

function generate_jacobian(sys::NonlinearSystem, vs = states(sys), ps = parameters(sys);
                           sparse = false, kwargs...)
    jac = calculate_jacobian(sys)
    if sparse
        jac = SparseArrays.sparse(jac)
    end
    return build_function(jac, convert.(Variable,vs), convert.(Variable,ps),
                          conv = AbstractSysToExpr(sys))
end

function generate_function(sys::NonlinearSystem, vs = states(sys), ps = parameters(sys); kwargs...)
    rhss = [eq.rhs for eq ∈ sys.eqs]
    vs′ = convert.(Variable,vs)
    ps′ = convert.(Variable,ps)
    return build_function(rhss, vs′, ps′;
                          conv = AbstractSysToExpr(sys), kwargs...)
end

"""
```julia
function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem,u0map,tspan,
                                          parammap=DiffEqBase.NullParameters();
                                          jac = false, sparse=false,
                                          checkbounds = false,
                                          linenumbers = true, multithread=false,
                                          kwargs...) where iip
```

Generates an NonlinearProblem from a NonlinearSystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem,u0map,tspan,
                                          parammap=DiffEqBase.NullParameters();
                                          jac = false, sparse=false,
                                          checkbounds = false,
                                          linenumbers = true, multithread=false,
                                          kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)

    f = generate_function(sys;checkbounds=checkbounds,linenumbers=linenumbers,
                              multithread=multithread,sparse=sparse,expression=Val{false})
    u0 = varmap_to_vars(u0map,dvs)
    p = varmap_to_vars(parammap,ps)
    NonlinearProblem(f,u0,tspan,p;kwargs...)
end


text:
primitive type Single32 <: AbstractFloat 64 end

Single32(x::Single32) = x

Single32(x::Float64)  = reinterpret(Single32, x)
Single32(x::Float32)  = Core.Intrinsics.fpext(Single32, x)
Single32(x::Float16)  = reinterpret(Single32, Float64(x))
Single32(x::BigFloat) = Single32(Float64(x))
Single32(x::BigInt)   = Single32(Float64(x))

Base.Float64(x::Single32) = reinterpret(Float64, x)
Base.Float32(x::Single32) = Core.Intrinsics.fptrunc(Float32, x)
Base.Float16(x::Single32) = Float16(Float32(x))

for st in (Int8, Int16, Int32, Int64)
    @eval begin
        (::Type{Single32})(x::($st)) = Core.Intrinsics.sitofp(Single32, x)
        promote_rule(::Type{Single32}, ::Type{$st}) = Single32
    end
end
for ut in (Bool, UInt8, UInt16, UInt32, UInt64)
    @eval begin
        (::Type{Single32})(x::($ut)) = Core.Intrinsics.uitofp(Single32, x)
        promote_rule(::Type{Single32}, ::Type{$ut}) = Single32
    end
end

for T in (Int128, UInt128)
  @eval begin
     promote_rule(::Type{Single32}, ::Type{$T}) = Single32
     Single32(x::$T) = Single32(Float64(x))
   end
end

Bool(x::Single32) = x==0 ? false : x==1 ? true : throw(InexactError(:Bool, Bool, x))

Single32(x::Irrational{S}) where S = Single32(Float64(x))

Single32(x::Complex{Float64}) = Single32(Float64(x))
Single32(x::Complex{Float32}) = Single32(Float32(x))

show(io::IO, x::Single32) = show(io, Float32(x))
string(x::Single32) = string(Float32(x))
repr(x::Single32) = string("Single32(",repr(Float64(x)),")")

widen(::Type{Single32}) = Float64
widen(x::Single32) = Float64(x)

hash(x::Single32) = hash(Float64(x))
hash(x::Single32, h::UInt64) = hash(Float64(x), h)

maxintfloat(::Type{Single32}) = maxintfloat(Float32)
maxintfloat(::Type{Single32}, ::Type{I}) where {I<:Integer} = maxintfloat(Float32, I)

signbit(x::Single32) = signbit(Float64(x))

zero(::Type{Single32}) = Single32(zero(Float64))
one(::Type{Single32}) = Single32(one(Float64))
iszero(x::Single32) = iszero(Float64(x))
isone(x::Single32) = isone(Float64(x))
isinteger(x::Single32) = isinteger(Float32(x))

prevfloat(x::Single32, n::Int) = Single32(prevfloat(Float32(x), n))
nextfloat(x::Single32, n::Int) = Single32(nextfloat(Float32(x), n))

typemin(::Type{Single32})  = typemin(Float32)
typemax(::Type{Single32})  = typemax(Float32)
floatmin(::Type{Single32}) = floatmin(Float32)
floatmax(::Type{Single32}) = floatmax(Float32)
typemin(x::Single32)  = typemin(Single32)
typemax(x::Single32)  = typemax(Single32)
floatmin(x::Single32) = floatmin(Single32)
floatmax(x::Single32) = floatmax(Single32)

eps(::Type{Single32}) = eps(Float32)
eps(x::Single32) = eps(Float32(x))
exponent(x::Single32) = exponent(Float64(x))
significand(x::Single32) = significand(Float64(x))
precision(::Type{Single32}) = precision(Float64)

trunc(x::Single32) = trunc(Float64(x))
floor(x::Single32) = floor(Float64(x))
ceil(x::Single32)  = ceil(Float64(x))

round(x::Single32, r::RoundingMode{:ToZero})  = trunc(Float64(x))
round(x::Single32, r::RoundingMode{:Down})    = floor(Float64(x))
round(x::Single32, r::RoundingMode{:Up})      = ceil(Float64(x))
round(x::Single32, r::RoundingMode{:Nearest}) = round(Float64(x))

trunc(::Type{I}, x::Single32) where {I<:Integer} = trunc(I, Float64(x))
floor(::Type{I}, x::Single32) where {I<:Integer} = floor(I, Float64(x))
ceil(::Type{I}, x::Single32)  where {I<:Integer} = ceil(I, Float64(x))

round(::Type{I}, x::Single32, r::RoundingMode{:ToZero})  where {I<:Integer} = trunc(I, Float64(x))
round(::Type{I}, x::Single32, r::RoundingMode{:Down})    where {I<:Integer} = floor(I, Float64(x))
round(::Type{I}, x::Single32, r::RoundingMode{:Up})      where {I<:Integer} = ceil(I, Float64(x))
round(::Type{I}, x::Single32, r::RoundingMode{:Nearest}) where {I<:Integer} = round(I, Float64(x))
round(::Type{I}, x::Single32) where {I<:Integer} = round(I, Float64(x))

Base.BigFloat(x::Single32) = BigFloat(Float64(x))
Base.Int128(x::Single32) = Int128(Float64(x))
Base.Int64(x::Single32) = Int64(Float64(x))
Base.Int32(x::Single32) = Int32(Float64(x))
Base.Int16(x::Single32) = Int16(Float64(x))

promote_rule(::Type{Single32}, ::Type{Float32}) = Single32
promote_rule(::Type{Single32}, ::Type{Float64}) = Single32

for T in (BigFloat, Float64, Float32, Float16,
          Int128, Int64, Int32, Int16, Int8,
          UInt128, UInt64, UInt32, UInt16, UInt8)
  @eval begin
    convert(::Type{Single32}, x::$T) = Single32(x)
    convert(::Type{$T}, x::Single32) = ($T)(x)
  end
end

Single32(x::Rational{T}) where {T} = Single32(Float64(x))
Rational{T}(x::Single32; tol=eps(x)) where {T} = Rational{T}(rationalize(Float64(x), tol=tol))

promote_rule(::Type{Single32}, ::Type{Rational}) = Single32
convert(::Type{Single32}, x::Rational) = Single32(x)

# comparison

==(x::Single32, y::Single32) = Float32(x) === Float32(y)
==(x::Single32, y::Float64)  = Float64(x) === y
==(x::Float64, y::Single32)  = x === Float64(y)
!=(x::Single32, y::Single32) = Float32(x) !== Float32(y)
!=(x::Single32, y::Float64)  = Float64(x) !== y
!=(x::Float64, y::Single32)  = x !== Float64(y)
<( x::Single32, y::Single32) = Float32(x) < Float32(y)
<( x::Single32, y::Float64)  = Float64(x) < y
<( x::Float64, y::Single32)  = x < Float64(y)
<=(x::Single32, y::Single32) = Float32(x) <= Float32(y)
<=(x::Single32, y::Float64)  = Float64(x) <= y
<=(x::Float64, y::Single32)  = x <= Float64(y)

isequal(x::Single32, y::Single32) = isequal(Float32(x), Float32(y))
isequal(x::Single32, y::Float64)  = isequal(Float64(x), y)
isequal(x::Float64, y::Single32)  = isequal(x, Float64(y))
isless( x::Single32, y::Single32) = isless(Float32(x), Float32(y))
isless( x::Single32, y::Float64)  = isless(Float64(x), y)
isless( x::Float64, y::Single32)  = isless(x, Float64(y))

==(x::Single32, y::Float32)  = Float32(x) === y
==(x::Float32, y::Single32)  = x === Float32(y)
!=(x::Single32, y::Float32)  = Float32(x) !== y
!=(x::Float32, y::Single32)  = x !== Float32(y)
<( x::Single32, y::Float32)  = Float32(x) < y
<( x::Float32, y::Single32)  = x < Float32(y)
<=(x::Single32, y::Float32)  = Float32(x) <= y
<=(x::Float32, y::Single32)  = x <= Float32(y)

isequal(x::Single32, y::Float32)  = isequal(Float32(x), y)
isequal(x::Float32, y::Single32)  = isequal(x, Float32(y))
isless( x::Single32, y::Float32)  = isless(Float32(x), y)
isless( x::Float32, y::Single32)  = isless(x, Float32(y))

for Op in (:cmp, :(==), :(!=), :(>=), :(<=), :(>), :(<), :isless, :isequal)
    @eval begin
        $Op(x::Single32, y::Float16) = $Op(Float64(x), Float64(y))
        $Op(x::Float16, y::Single32) = $Op(Float64(x), Float64(y))

        $Op(x::Single32, y::BigFloat) = $Op(Float64(x), y)
        $Op(x::BigFloat, y::Single32) = $Op(x, Float64(y))

        $Op(x::Single32, y::Int128) = $Op(Float64(x), y)
        $Op(x::Int128, y::Single32) = $Op(x, Float64(y))
        $Op(x::Single32, y::Int64)  = $Op(Float64(x), y)
        $Op(x::Int64, y::Single32)  = $Op(x, Float64(y))
        $Op(x::Single32, y::Int32)  = $Op(Float64(x), y)
        $Op(x::Int32, y::Single32)  = $Op(x, Float64(y))
        $Op(x::Single32, y::Int16)  = $Op(Float64(x), y)
        $Op(x::Int16, y::Single32)  = $Op(x, Float64(y))

        $Op(x::Single32, y::UInt128) = $Op(Float64(x), y)
        $Op(x::UInt128, y::Single32) = $Op(x, Float64(y))
        $Op(x::Single32, y::UInt64)  = $Op(Float64(x), y)
        $Op(x::UInt64, y::Single32)  = $Op(x, Float64(y))
        $Op(x::Single32, y::UInt32)  = $Op(Float64(x), y)
        $Op(x::UInt32, y::Single32)  = $Op(x, Float64(y))
        $Op(x::Single32, y::UInt16)  = $Op(Float64(x), y)
        $Op(x::UInt16, y::Single32)  = $Op(x, Float64(y))
    end
end

decompose(x::Single32) = decompose(Float32(x))


text:
ENV["GKSwstype"] = "nul"

using Plots

using Printf
using Statistics
using CUDA

using Oceananigans
using Oceananigans.Advection
using Oceananigans.Units

using Oceananigans.Advection: EnergyConservingScheme
using Oceananigans.OutputReaders: FieldTimeSeries

using Oceananigans.Advection: ZWENO, WENOVectorInvariantVel, WENOVectorInvariantVort, VectorInvariant, VelocityStencil, VorticityStencil
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.Operators: Δx, Δy
using Oceananigans.TurbulenceClosures

CUDA.device!(1)

include("bickley_utils.jl")

"""
    run_bickley_jet(output_time_interval = 2, stop_time = 200, arch = CPU(), Nh = 64, ν = 0,
                    momentum_advection = VectorInvariant())

Run the Bickley jet validation experiment until `stop_time` using `momentum_advection`
scheme or formulation, with horizontal resolution `Nh`, viscosity `ν`, on `arch`itecture.
"""
function run_bickley_jet(; output_time_interval = 2, stop_time = 200, arch = CPU(), Nh = 64,
                           momentum_advection = VectorInvariant())

    grid = bickley_grid(; arch, Nh, halo=(4, 4, 4))

    @inline toplft(x, y) = (((x > π/2) & (x < 3π/2)) & (((y > π/3) & (y < 2π/3)) | ((y > 4π/3) & (y < 5π/3))))
    @inline botlft(x, y) = (((x > π/2) & (x < 3π/2)) & (((y < -π/3) & (y > -2π/3)) | ((y < -4π/3) & (y > -5π/3))))
    @inline toprgt(x, y) = (((x < -π/2) & (x > -3π/2)) & (((y > π/3) & (y < 2π/3)) | ((y > 4π/3) & (y < 5π/3))))
    @inline botrgt(x, y) = (((x < -π/2) & (x > -3π/2)) & (((y < -π/3) & (y > -2π/3)) | ((y < -4π/3) & (y > -5π/3))))
    @inline bottom(x, y) = Int(toplft(x, y) | toprgt(x, y) | botlft(x, y) | botrgt(x, y))

    grid = ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> -1))

    c = sqrt(10.0)
    Δt = 0.1 * grid.Δxᶜᵃᵃ / c

    timescale = (5days / (6minutes) * Δt)
    @show prettytime(timescale)

    @inline νhb(i, j, k, grid, lx, ly, lz) = (1 / (1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2 ))^2 / timescale
    biharmonic_viscosity = HorizontalScalarBiharmonicDiffusivity(ν=νhb, discrete_form=true)

    model = HydrostaticFreeSurfaceModel(momentum_advection = momentum_advection,
                                        tracer_advection = WENO5(),
                                        grid = grid,
                                        tracers = :c,
                                        closure = nothing,
                                        free_surface = ExplicitFreeSurface(gravitational_acceleration=10.0),
                                        coriolis = nothing,
                                        buoyancy = nothing)

    # ** Initial conditions **
    #
    # u, v: Large-scale jet + vortical perturbations
    #    c: Sinusoid

    set_bickley_jet!(model)

    wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=10.0)

    simulation = Simulation(model, Δt=Δt, stop_time=stop_time)

    progress(sim) = @printf("Iter: %d, time: %s, Δt: %s, max|u|: %.3f, max|η|: %.3f \n",
                            iteration(sim), prettytime(sim), prettytime(sim.Δt),
                            maximum(abs, model.velocities.u), maximum(abs, model.free_surface.η))

    simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
    wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=10.0)

    simulation.callbacks[:wizard]   = Callback(wizard, IterationInterval(10))

    # Output: primitive fields + computations
    u, v, w, c = merge(model.velocities, model.tracers)

    ζ = Field(∂x(v) - ∂y(u))

    outputs = merge(model.velocities, model.tracers, (ζ=ζ, η=model.free_surface.η))

    name = typeof(model.advection.momentum).name.wrapper
    if model.advection.momentum isa WENOVectorInvariantVel
        name = "WENOVectorInvariantVel"
    end
    if model.advection.momentum isa WENOVectorInvariantVort
        name = "WENOVectorInvariantVort"
    end

    @show experiment_name = "bickley_jet_Nh_$(Nh)_Upwind"

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, outputs,
                                schedule = TimeInterval(output_time_interval),
                                prefix = experiment_name,
                                force = true)

    @info "Running a simulation of an unstable Bickley jet with $(Nh)² degrees of freedom..."

    start_time = time_ns()

    run!(simulation)

    return experiment_name
end

"""
    visualize_bickley_jet(experiment_name)

Visualize the Bickley jet data associated with `experiment_name`.
"""
function visualize_bickley_jet(experiment_name)

    @info "Making a fun movie about an unstable Bickley jet..."

    filepath = experiment_name * ".jld2"

    ζ_timeseries = FieldTimeSeries(filepath, "ζ", boundary_conditions=nothing, location=(Face, Face, Center))
    c_timeseries = FieldTimeSeries(filepath, "c", boundary_conditions=nothing, location=(Face, Center, Center))

    grid = c_timeseries.grid

    xζ, yζ, zζ = nodes(ζ_timeseries)
    xc, yc, zc = nodes(c_timeseries)

    anim = @animate for (i, iteration) in enumerate(c_timeseries.times)

        @info "    Plotting frame $i from iteration $iteration..."

        ζ = ζ_timeseries[i]
        c = c_timeseries[i]
        t = ζ_timeseries.times[i]

        ζi = interior(ζ)[:, :, 1]
        ci = interior(c)[:, :, 1]

        kwargs = Dict(
                      :aspectratio => 1,
                      :linewidth => 0,
                      :colorbar => :none,
                      :ticks => nothing,
                      :clims => (-1, 1),
                      :xlims => (-grid.Lx/2, grid.Lx/2),
                      :ylims => (-grid.Ly/2, grid.Ly/2)
                     )

        ζ_plot = heatmap(xζ, yζ, clamp.(ζi, -1, 1)'; color = :balance, kwargs...)
        c_plot = heatmap(xc, yc, clamp.(ci, -1, 1)'; color = :thermal, kwargs...)

        ζ_title = @sprintf("ζ at t = %.1f", t)
        c_title = @sprintf("u at t = %.1f", t)

        plot(ζ_plot, c_plot, title = [ζ_title c_title], size = (4000, 2000))
    end

    mp4(anim, experiment_name * ".mp4", fps = 8)
end

advection_schemes = [WENO5(vector_invariant=VelocityStencil()),
                     WENO5(vector_invariant=VorticityStencil()),
                     WENO5(),
                     VectorInvariant()]

for Nx in [128]
    for advection in advection_schemes
        experiment_name = run_bickley_jet(arch=GPU(), momentum_advection=advection, Nh=Nx)
        # visualize_bickley_jet(experiment_name)
    end
end