text:
    using StatGeochem
    using Plots
    using LsqFit: curve_fit

## --- Load partition coefficients from GERM into a dict

    # path = joinpath(moduleresourcepath,"PartitionCoefficients")
    path = joinpath("src/resources/","PartitionCoefficients")

    # Obtain the names of downloaded mineral files
    filenames = readdir(path)

    # Load every `.tsv` file into a struct
    pr = Dict{String,Any}()
    pr["minerals"] = String[]
    for n in filenames
        m = match(r"(.*).tsv", n)
        if ~ isnothing(m)
            name = m[1]
            push!(pr["minerals"], name)
            pr[name] = importdataset(joinpath(path, name*".tsv"), '\t', mindefinedcolumns=2)
        end
    end

    # Approximate content of various rock types, for when actual melt SiO2 not reported
    str = ["basalt","andesit","dacit","rhyolit","hawaiite","alkali basalt","basanit","benmoreite","camptonite","eclogit","garnet pyroxenit","kimberlit","komatiit","lampro","latite","leucitit","lherzolit","mugearite","pantellerite","trachyt","peridotit","phonolit","picrit","syen","morb","tholeiit","tonalit","granodiorit","granit","shoshonit","carbonatit","aplite","monzonit","synthetic"]
    val = [49,57.5,67.5,73,48.5,47,44,57,42,49,49,40,45,46.5,62.5,47,42,52,73,61,42,59,45.5,61,49,50,65,67.5,72,55,2,73,67.4,60]
    err = [5,6.5,5.5,5,3.5,4,3,4,3,5,5,4,5,11.5,11.5,5,5,5,5,9,5,5,5.5,9,5,4,8,7.5,5,5,5,5,5,20]

    # Determine the unique reference / rock type pairs
    samples = Dict{String,Vector{Tuple{String,String,Float64,Float64}}}()
    allsamples = Tuple{String,String,Float64,Float64}[]
    allelements = String[]
    for m in pr["minerals"]
        if ~haskey(pr[m],"SiO2")
            pr[m]["SiO2"] = fill(NaN, size(pr[m]["Reference"]))
        end
        if ~haskey(pr[m],"SiO2_sigma")
            pr[m]["SiO2_sigma"] = fill(NaN, size(pr[m]["Reference"]))
        end
        # Estimate SiO2 from rock type
        needssi = isnan.(pr[m]["SiO2"])
        for i=1:length(str)
            t = containsi.(pr[m]["Rock Type"], str[i]) .& needssi
            # Average together, such that "basaltic andesite" is halfway between "basalt" and "andesite"
            pr[m]["SiO2"][t] = nanmean([pr[m]["SiO2"][t] val[i]*ones(count(t))], dim=2)
            pr[m]["SiO2_sigma"][t] = nanmean([pr[m]["SiO2_sigma"][t] err[i]*ones(count(t))], dim=2)
        end

        # Save unique samples
        println(m)
        flush(stdout)
        samples[m] = collect(zip(string.(pr[m]["Reference"]), string.(pr[m]["Rock Type"]), pr[m]["SiO2"], pr[m]["SiO2_sigma"]))
        append!(allsamples, unique(samples[m]))

        # Save unique elements
        append!(allelements, unique(pr[m]["Elem"]))
    end
    allsamples = sort(unique(allsamples))
    allelements = sort(unique(allelements))

    # Second dict: sort the data by sample ID
    pd = Dict{String, Union{Vector, Dict{String, Union{Vector{Float64}, Vector{String}}}}}()
    pd["samples"] = allsamples
    pd["Reference"] = allsamples .|> x -> x[1]
    pd["Rock Type"] = allsamples .|> x -> x[2]
    pd["SiO2"] = allsamples .|> x -> x[3]
    pd["SiO2_sigma"] = allsamples .|> x -> x[4]
    pd["minerals"] = pr["minerals"]
    for m in pd["minerals"]
        pd[m] = Dict{String, Union{Vector{Float64}, Vector{String}}}()
        pd[m]["elements"] = allelements
        for e in allelements
            pd[m][e] = fill(NaN, length(allsamples))
            pd[m][e*"_sigma"] = fill(NaN, length(allsamples))

            t = pr[m]["Elem"] .== e
            kd = map(x->(x>0 ? log10(x) : NaN), pr[m]["Value"][t])
            kd_sigma = log10.(pr[m]["Value"][t] .+ pr[m]["SD"][t]) .- kd
            kdl = map(x->(x>0 ? log10(x) : NaN), pr[m]["Low"][t])
            kdh = map(x->(x>0 ? log10(x) : NaN), pr[m]["High"][t])

            if any(t)
                rows = findmatches(samples[m][t], allsamples)
                pd[m][e][rows] = nanmean([kd nanmean([kdl kdh], dim=2)], dim=2)
                pd[m][e*"_sigma"][rows] = nanmean([kd_sigma (kdh .- kdl)/4], dim=2)
            end
        end
    end

## --- Fit REEs as a function of activation energy and bulk modulus

    ree3 = ["La","Pr","Nd","Sm","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu",] # 3+ rare earth elements
    r = [0.106,0.101,0.100,0.096,0.094,0.092,0.091,0.089,0.088,0.087,0.086,0.085,] # Atomic Radii
    rp = 0.085:0.001:0.107 # Radius range to plot over

    # Equation we're fitting to (from Blundy and Wood, 1994):
    # lnD0 + a * (r0/2*(x-r0)^2 + 1/3*(x-r0)^3)
    # where a = 4π E Na / RT
    @. blundy_wood(x,param) = param[1] + param[2] * (param[3]/2*(x-param[3])^2 + 1/3*(x-param[3])^3)

    # for mineral={'Apatite','Amphibole','Clinopyroxene','Orthopyroxene','Garnet','Sphene','Allanite','Baddeleyite'}
    for m in pd["minerals"]
        h = plot(title=m, legend=:bottomleft)
        for i = 1:length(pd["samples"])
            kD = fill(NaN,length(ree3))
            for j = 1:length(ree3)
                kD[j] = pd[m][ree3[j]][i]
            end

            if count(.~isnan.(kD)) > 1
                # plot!(h, r, kD, label="$i", seriestype=:scatter, color=lines[mod(i,length(lines))+1], msalpha=0)
                plot!(h, r, kD, label="", seriestype=:scatter, color=lines[mod(i,length(lines))+1], msalpha=0)
            end

            if (count(.~isnan.(kD)) > 3) && (nanrange(r[.~isnan.(kD)]) > 0.013)
                # Fit to Blundy and Wood curve
                t = .~( isnan.(kD) .| isinf.(kD) )
                param = [maximum(kD[t]), -10000, 0.095] # initial guess
                f = curve_fit(blundy_wood, r[t], kD[t], param) # Fit
                plot!(h, r, blundy_wood(r,f.param), label="", color=lines[mod(i,length(lines))+1]) # Plot

                # Replace stored partition coefficients with new fits
                for j = 1:length(ree3)
                    pd[m][ree3[j]][i] = blundy_wood(r[j],f.param)
                end
            end
        end
        savefig(h, "Calc_$(m)_REE.pdf")
    end

    # Interpolate Eu partition coefficients where missing, assuming
    # 60% Eu as Eu2+ (c.f. Ba, Sr, Ca) and 40% as Eu3+ (c.f. Sm, Gd)
    # for m in = ["Albite","Anorthite","Orthoclase","Apatite"]
    for m in pd["minerals"]
        for i = 1:length(pd["samples"])
            if isnan(pd[m]["Eu"][i])
                pd[m]["Eu"][i] = log10(0.6*10^nanmean([pd[m]["Ba"][i], pd[m]["Sr"][i]]) + 0.4*10^nanmean([pd[m]["Sm"][i], pd[m]["Gd"][i]]))
            end
        end
    end
    for m in ("Monazite", "Xenotime", "Allanite")
        for i = 1:length(pd["samples"])
            if isnan(pd[m]["Eu"][i])
                pd[m]["Eu"][i] = log10(0.6*0 + 0.4*10^nanmean([pd[m]["Sm"][i], pd[m]["Gd"][i]]))
            end
        end
    end

## --- Convert data to average as a function of SiO2

    # Convert from row-based to Si-based Dict
    kd = Dict{String, Union{Vector{String}, Vector{Float64}, Dict{String,Union{Float64, Vector{String}, Vector{Float64}}}}}()
    kd["minerals"] = pd["minerals"]
    for m in kd["minerals"]
        kd[m] = Dict{String,Union{Float64, Vector{String}, Vector{Float64}}}()
        kd[m]["elements"] = allelements
        for e in allelements
            t = .!isnan.(pd[m][e])
            if (count(t) > 2) && (nanrange(pd["SiO2"][t]) > 8)
                kd[m][e] = mcfit(pd["SiO2"], pd["SiO2_sigma"], pd[m][e], pd[m][e*"_sigma"], 40, 80, 41, binwidth=5)[2]
            else
                kd[m][e] = ones(41) * nanmean(pd[m][e])
            end
            kd[m][e*"_sigma"] = nanstd(pd[m][e])
        end
    end
    kd["SiO2"] = collect(40:80.)

    # Set Albite partiton coefficients
    for e in kd["Albite"]["elements"]
        kd["Albite"][e] = nanmean([kd["Albite"][e] kd["Orthoclase"][e] kd["Anorthite"][e]], dim=2)
    end
    kd["note"] = ["kd for Albite is nanmean of AlkaliFeldspar, Orthoclase, and Anorthite",]

## --- Save results

    # using MAT
    # matwrite(joinpath(path,"partitioncoeffs.mat"),p)
    f = open(joinpath(path,"PartitionCoefficients.jl"), "a")
    print(f, "\ngerm_kd = $kd\nexport germ_kd\n")
    close(f)

## --- End of File


text:
using Test, MathJaxRenderer, VisualRegressionTests

# Test cases from Latexify.jl.

keywords = Dict(
    "small" => (
        x_zoom = 4,
        y_zoom = 4,
        background_color = :orange,
    ),
    "large" => (
        x_zoom = 2,
        y_zoom = 2,
        background_color = :lightblue,
    ),
)

@testset "MathJaxRenderer" begin
    @testset "Integration" begin
        tex = joinpath(@__DIR__, "data", "tex")
        for file in readdir(tex)
            name = first(splitext(file))
            file = joinpath(tex, file)
            math = read(file, Math)
            kws = keywords[name]
            @test !isempty(sprint(io -> show(io, MIME("image/svg+xml"), math; kws...)))
            @test !isempty(sprint(io -> show(io, MIME("image/png"), math; kws...)))
            @test !isempty(sprint(io -> show(io, MIME("application/pdf"), math; kws...)))
            @test !isempty(sprint(io -> show(io, MIME("application/postscript"), math; kws...)))
        end
    end
    @testset "Visuals" begin
        testfun = (filename) -> write(filename, Math("\\frac{1}{1 + x}"); x_zoom = 5, y_zoom = 5)
        @visualtest testfun joinpath(@__DIR__, "data/references/fraction.png") popup = false tol = 0.2
        testfun = (filename) -> write(filename, Math("\\begin{align}\n3.29e+07 =& 1.23e+00 =& P_{1} \\\\\n\\frac{x}{y} =& 1.00e+10 =& 1.29e+03\n\\end{align}"); x_zoom = 5, y_zoom = 5)
        @visualtest testfun joinpath(@__DIR__, "data/references/align.png") popup = false tol = 0.2
    end
end

ext: jl


text:
module MyExamplePkg

greet() = println("Hello")

include("extra_file.jl")
export my_f

end


text:
push!(LOAD_PATH, "../src")
using Test
using ImplicitGlobalGrid; GG = ImplicitGlobalGrid
import MPI
import ImplicitGlobalGrid: @require
macro coords(i) :(GG.global_grid().coords[$i]) end

## Test setup
MPI.Init();
nprocs = MPI.Comm_size(MPI.COMM_WORLD);
@require nprocs == 1  # This test requires nprocs to be 1.

@testset "$(basename(@__FILE__))" begin
    @testset "1. *_g functions" begin
        lx = 8;
        ly = 8;
        lz = 8;
        nx = 5;
        ny = 5;
        nz = 5;
        P   = zeros(nx,  ny,  nz  );
        Vx  = zeros(nx+1,ny,  nz  );
        Vz  = zeros(nx,  ny,  nz+1);
        A   = zeros(nx,  ny,  nz+2);
        Sxz = zeros(nx-2,ny-1,nz-2);
        init_global_grid(nx, ny, nz, dimx=1, dimy=1, dimz=1, periodz=1, quiet=true, init_MPI=false);
        @testset "nx_g / ny_g / nz_g" begin
            @test nx_g() == nx
            @test ny_g() == ny
            @test nz_g() == nz-2
        end;
        @testset "x_g / y_g / z_g" begin
            dx  = lx/(nx_g()-1);
            dy  = ly/(ny_g()-1);
            dz  = lz/(nz_g()-1);
            # (for P)
            @test [x_g(ix,dx,P) for ix = 1:size(P,1)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [y_g(iy,dy,P) for iy = 1:size(P,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [z_g(iz,dz,P) for iz = 1:size(P,3)] == [8.0, 0.0, 4.0, 8.0, 0.0]
            # (for Vx)
            @test [x_g(ix,dx,Vx) for ix = 1:size(Vx,1)] == [-1.0, 1.0, 3.0, 5.0, 7.0, 9.0]
            @test [y_g(iy,dy,Vx) for iy = 1:size(Vx,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [z_g(iz,dz,Vx) for iz = 1:size(Vx,3)] == [8.0, 0.0, 4.0, 8.0, 0.0]
            # (for Vz)
            @test [x_g(ix,dx,Vz) for ix = 1:size(Vz,1)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [y_g(iy,dy,Vz) for iy = 1:size(Vz,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [z_g(iz,dz,Vz) for iz = 1:size(Vz,3)] == [ 6.0, 10.0,  2.0,  6.0, 10.0,  2.0]
            #                       base grid (z dim):        [ 8.0,  0.0,  4.0,  8.0,  0.0]
            #                       possible alternative:  [ 6.0, -2.0,  2.0,  6.0, -2.0,  2.0]  # This would be a possible alternative way to define {x,y,z}_g; however, we decided that the grid should start at 0.0 in this case and the overlap be at the end (we avoid completely any negative z_g).
            #                       wrong:                 [ 6.0, -2.0,  2.0,  6.0, 10.0,  2.0]  # The 2nd and the 2nd-last cell must be the same due to the overlap of 3.
            # (for A)
            @test [x_g(ix,dx,A) for ix = 1:size(A,1)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [y_g(iy,dy,A) for iy = 1:size(A,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [z_g(iz,dz,A) for iz = 1:size(A,3)] == [4.0, 8.0, 0.0, 4.0, 8.0, 0.0, 4.0]
            #                       base grid (z dim):       [ 8.0, 0.0, 4.0, 8.0, 0.0]
            # (for Sxz)
            @test [x_g(ix,dx,Sxz) for ix = 1:size(Sxz,1)] ==   [2.0, 4.0, 6.0]
            #                       base grid (x dim):    [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [y_g(iy,dy,Sxz) for iy = 1:size(Sxz,2)] == [1.0, 3.0, 5.0, 7.0]
            #                       base grid (y dim):    [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [z_g(iz,dz,Sxz) for iz = 1:size(Sxz,3)] ==   [0.0, 4.0, 8.0]
            #                       base grid (z dim):   [ 8.0, 0.0, 4.0, 8.0, 0.0]
        end;
        finalize_global_grid(finalize_MPI=false);
    end;

    @testset "2. *_g functions with non-default overlap" begin
        lx = 8;
        ly = 8;
        lz = 8;
        nx = 5;
        ny = 5;
        nz = 8;
        P   = zeros(nx,  ny,  nz  );
        Vx  = zeros(nx+1,ny,  nz  );
        Vz  = zeros(nx,  ny,  nz+1);
        A   = zeros(nx,  ny,  nz+2);
        Sxz = zeros(nx-2,ny-1,nz-2);
        init_global_grid(nx, ny, nz, dimx=1, dimy=1, dimz=1, periodz=1, overlapx=3, overlapz=3, quiet=true, init_MPI=false);
        @testset "nx_g / ny_g / nz_g" begin
            @test nx_g() == nx
            @test ny_g() == ny
            @test nz_g() == nz-3
        end;
        @testset "x_g / y_g / z_g" begin
            dx  = lx/(nx_g()-1);
            dy  = ly/(ny_g()-1);
            dz  = lz/(nz_g()-1);
            # (for P)
            @test [x_g(ix,dx,P) for ix = 1:size(P,1)] == [0.0, 2.0, 4.0, 6.0, 8.0]          # (same as in the first test)
            @test [y_g(iy,dy,P) for iy = 1:size(P,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]          # (same as in the first test)
            @test [z_g(iz,dz,P) for iz = 1:size(P,3)] == [8.0, 0.0, 2.0, 4.0, 6.0, 8.0, 0.0, 2.0]
            # (for Vz)
            @test [x_g(ix,dx,Vz) for ix = 1:size(Vz,1)] == [0.0, 2.0, 4.0, 6.0, 8.0]        # (same as in the first test)
            @test [y_g(iy,dy,Vz) for iy = 1:size(Vz,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]        # (same as in the first test)
            @test [z_g(iz,dz,Vz) for iz = 1:size(Vz,3)] == [7.0, 9.0, 1.0, 3.0, 5.0, 7.0, 9.0, 1.0, 3.0]
            #                       base grid (z dim):       [8.0, 0.0, 2.0, 4.0, 6.0. 8.0, 0.0, 2.0]
            #                       possible alternative:  [7.0,-1.0, 1.0, 3.0, 5.0, 7.0,-1.0, 1.0, 3.0]
            # (for A)
            @test [x_g(ix,dx,A) for ix = 1:size(A,1)] == [0.0, 2.0, 4.0, 6.0, 8.0]          # (same as in the first test)
            @test [y_g(iy,dy,A) for iy = 1:size(A,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]          # (same as in the first test)
            @test [z_g(iz,dz,A) for iz = 1:size(A,3)] == [6.0, 8.0, 0.0, 2.0, 4.0, 6.0, 8.0, 0.0, 2.0, 4.0]
            #                       base grid (z dim):        [8.0, 0.0, 2.0, 4.0, 6.0, 8.0, 0.0, 2.0]
            # (for Sxz)
            @test [x_g(ix,dx,Sxz) for ix = 1:size(Sxz,1)] ==   [2.0, 4.0, 6.0]         # (same as in the first test)
            #                       base grid (x dim):    [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [y_g(iy,dy,Sxz) for iy = 1:size(Sxz,2)] == [1.0, 3.0, 5.0, 7.0]       # (same as in the first test)
            #                       base grid (y dim):    [0.0, 2.0, 4.0, 6.0, 8.0]
            @test [z_g(iz,dz,Sxz) for iz = 1:size(Sxz,3)] ==   [0.0, 2.0, 4.0, 6.0, 8.0, 0.0]
            #                       base grid (z dim):    [8.0, 0.0, 2.0, 4.0, 6.0, 8.0, 0.0, 2.0]
        end;
        finalize_global_grid(finalize_MPI=false);
    end;

    @testset "3. *_g functions (simulated 3x3x3 processes)" begin
        lx = 20;
        ly = 20;
        lz = 16;
        nx = 5;
        ny = 5;
        nz = 5;
        P  = zeros(nx,  ny,  nz  );
        A  = zeros(nx+1,ny-2,nz+2);
        init_global_grid(nx, ny, nz, dimx=1, dimy=1, dimz=1, periodz=1, quiet=true, init_MPI=false);
        # (Set dims, nprocs and nxyz_g in GG.global_grid().)
        dims      = [3,3,3];
        nxyz      = GG.global_grid().nxyz;
        periods   = GG.global_grid().periods;
        overlaps  = GG.global_grid().overlaps;
        nprocs    = prod(dims);
        nxyz_g    = dims.*(nxyz.-overlaps) .+ overlaps.*(periods.==0);
        GG.global_grid().dims   .= dims;
        GG.global_grid().nxyz_g .= nxyz_g;
        @testset "nx_g / ny_g / nz_g" begin
            @test nx_g() == nxyz_g[1]
            @test ny_g() == nxyz_g[2]
            @test nz_g() == nxyz_g[3]
        end;
        @testset "x_g / y_g / z_g" begin
            dx  = lx/(nx_g()-1);
            dy  = ly/(ny_g()-1);
            dz  = lz/(nz_g()-1);
            # (for P)
            @coords(1)=0;  @test [x_g(ix,dx,P) for ix = 1:size(P,1)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @coords(1)=1;  @test [x_g(ix,dx,P) for ix = 1:size(P,1)] == [6.0, 8.0, 10.0, 12.0, 14.0]
            @coords(1)=2;  @test [x_g(ix,dx,P) for ix = 1:size(P,1)] == [12.0, 14.0, 16.0, 18.0, 20.0]
            @coords(2)=0;  @test [y_g(iy,dy,P) for iy = 1:size(P,2)] == [0.0, 2.0, 4.0, 6.0, 8.0]
            @coords(2)=1;  @test [y_g(iy,dy,P) for iy = 1:size(P,2)] == [6.0, 8.0, 10.0, 12.0, 14.0]
            @coords(2)=2;  @test [y_g(iy,dy,P) for iy = 1:size(P,2)] == [12.0, 14.0, 16.0, 18.0, 20.0]
            @coords(3)=0;  @test [z_g(iz,dz,P) for iz = 1:size(P,3)] == [16.0, 0.0, 2.0, 4.0, 6.0]
            @coords(3)=1;  @test [z_g(iz,dz,P) for iz = 1:size(P,3)] == [4.0, 6.0, 8.0, 10.0, 12.0]
            @coords(3)=2;  @test [z_g(iz,dz,P) for iz = 1:size(P,3)] == [10.0, 12.0, 14.0, 16.0, 0.0]
            # (for A)
            @coords(1)=0;  @test [x_g(ix,dx,A) for ix = 1:size(A,1)] == [-1.0, 1.0, 3.0, 5.0, 7.0, 9.0]
            @coords(1)=1;  @test [x_g(ix,dx,A) for ix = 1:size(A,1)] == [5.0, 7.0, 9.0, 11.0, 13.0, 15.0]
            @coords(1)=2;  @test [x_g(ix,dx,A) for ix = 1:size(A,1)] == [11.0, 13.0, 15.0, 17.0, 19.0, 21.0]
            @coords(2)=0;  @test [y_g(iy,dy,A) for iy = 1:size(A,2)] == [2.0, 4.0, 6.0]
            @coords(2)=1;  @test [y_g(iy,dy,A) for iy = 1:size(A,2)] == [8.0, 10.0, 12.0]
            @coords(2)=2;  @test [y_g(iy,dy,A) for iy = 1:size(A,2)] == [14.0, 16.0, 18.0]
            @coords(3)=0;  @test [z_g(iz,dz,A) for iz = 1:size(A,3)] == [14.0, 16.0, 0.0, 2.0, 4.0, 6.0, 8.0]
            @coords(3)=1;  @test [z_g(iz,dz,A) for iz = 1:size(A,3)] == [2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0]
            @coords(3)=2;  @test [z_g(iz,dz,A) for iz = 1:size(A,3)] == [8.0, 10.0, 12.0, 14.0, 16.0, 0.0, 2.0]
        end;
        finalize_global_grid(finalize_MPI=false);
    end;
end;

## Test tear down
MPI.Finalize()


text:
module Example

using MemoryAccesses

function schedule_timings_dac(indices::AbstractUnitRange)
    if length(indices) == 0
    elseif length(indices) == 1
        MemoryAccesses.record(indices[1])
    else
        f = first(indices)
        l = last(indices)
        m = (l - f + 1) ÷ 2 + f - 1
        task = Threads.@spawn schedule_timings_dac(m+1:l)
        schedule_timings_dac(f:m)
        wait(task)
    end
    return nothing
end

function schedule_timings(nitems::Integer)
    MemoryAccesses.init(max(2^10, nitems))
    schedule_timings_dac(1:nitems)
    @assert !MemoryAccesses.isfull()
    return MemoryAccesses.copy()
end

function two_arrays_dac(indices::AbstractUnitRange, xs, ys)
    if length(indices) == 0
    elseif length(indices) == 1
        i = indices[1]
        MemoryAccesses.record(@inbounds view(xs, i:i))
        MemoryAccesses.record(@inbounds view(ys, i:i))
    else
        f = first(indices)
        l = last(indices)
        m = (l - f + 1) ÷ 2 + f - 1
        task = Threads.@spawn two_arrays_dac(m+1:l, xs, ys)
        two_arrays_dac(f:m, xs, ys)
        wait(task)
    end
    return nothing
end

function two_arrays(nitems::Integer)
    xs = zeros(nitems)
    ys = zeros(nitems)
    MemoryAccesses.init(max(2^10, nitems))
    MemoryAccesses.reference(xs, "xs")
    MemoryAccesses.reference(ys, "ys")
    two_arrays_dac(1:nitems, xs, ys)
    @assert !MemoryAccesses.isfull()
    return MemoryAccesses.copy()
end

end  # module


text:
# Autogenerated wrapper script for libfreenect_jll for x86_64-unknown-freebsd
export libfreenect, libfreenect_sync

using libusb_jll
JLLWrappers.@generate_wrapper_header("libfreenect")
JLLWrappers.@declare_library_product(libfreenect, "libfreenect.so.0")
JLLWrappers.@declare_library_product(libfreenect_sync, "libfreenect_sync.so.0")
function __init__()
    JLLWrappers.@generate_init_header(libusb_jll)
    JLLWrappers.@init_library_product(
        libfreenect,
        "lib/libfreenect.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfreenect_sync,
        "lib/libfreenect_sync.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

function randlines(n=300, t=0.4)
	#n = 100 #1000 #1000 #20000
	#t = 0.4 #0.15 #0.4 #0.15
	V = zeros(Float64,2,2*n)
	EV = [zeros(Int64,2) for k=1:n]
	for k=1:n
		v1 = rand(Float64,2)
		v2 = rand(Float64,2)
		vm = (v1+v2)/2
		transl = rand(Float64,2)
		V[:,k] = (v1-vm)*t + transl
		V[:,n+k] = (v2-vm)*t + transl
		EV[k] = [k,n+k]
	end

	V = GL.normalize2(V)
	model = (V,EV)
	Sigma = Lar.spaceindex(model)

	model = V,EV;
	W,EW = Lar.fragmentlines(model);
	U,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells));
	EV = convert(Lar.Cells, union(EVs...))
	V,FVs,EVs = Lar.arrange2D(U,EV)
end

# ////////////////////////////////////////////////////////////
# generation of 2D arrangement
V,FVs,EVs = randlines()

# native OpenGL visualization
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,3,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,1));
GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,1,1));


text:
using Unitful
export inradians, inseconds, inframes, inHz

# having these in a submodule allows the user to decide if
# they want to import the shorthand unit names (or just use e.g. u"ms")
module Units
    using Unitful
    using Unitful: Hz, s, kHz, ms, °, rad, dB

    @dimension 𝐒 "𝐒" Frame

    @refunit frames "frames" Frames 𝐒 true
    const SampDim = Unitful.Dimensions{(Unitful.Dimension{:Frame}(1//1),)}
    const FrameQuant{N} = Quantity{N,𝐒}

    const localunits = Unitful.basefactors
    const localpromotion = Unitful.promotion
    function __init__()
        merge!(Unitful.basefactors,localunits)
        merge!(Unitful.promotion, localpromotion)
    end

    export kframes, frames, Hz, s, kHz, ms, dB, °, rad
end
using .Units
using .Units: FrameQuant

"""
    inradians([Type],x)

Given an angle value, convert to a value of Type (defaults to Float64) in
radians. Unitless numbers are assumed to be in radians and are silently
passed through.

## Examples

```jldoctest
julia> inradians(180°)
3.141592653589793

julia> inradians(2π)
6.283185307179586

julia> inradians(0.5π*rad)
1.5707963267948966
```

"""
inradians(x::Number, rate=missing) = x
inradians(x::Quantity,rate=missing) = ustrip(uconvert(rad, x))
inradians(x::Unitful.Time, rate=missing) =
    ustrip(uconvert(Unitful.NoUnits, x * inHz(rate)*Hz))*2π * rad
inradians(::Type{T},x::Number, rate=missing) where T <: Integer =
    trunc(T,inradians(x,rate))
inradians(::Type{T},x::Number, rate=missing) where T =
    convert(T,inradians(x,rate))

"""
    inframes([Type,]quantity[, rate])

Translate the given quantity to a (unitless) number of time frames,
given a particular framerate. Note that this isn't quantized to integer numbers
of frames. If given a `Type`, the result will first be coerced to the given type.

If the given quantity is Unitful, we use the given units. If it is not we assume
it is already a value in frames.

For some units (e.g. frames) you will need to specify a frame rate. If not
specified the rate is `missing`.

# Example

```jldoctest
julia> inframes(0.5s, 44100Hz)
22050.0

julia> inframes(Int,10.5frames)
10
```

"""
inframes(frame::FrameQuant, rate=missing) = ustrip(uconvert(frames, frame))
inframes(time::Unitful.Time, rate=missing) = inseconds(time)*inHz(rate)
inframes(frame::Number, rate=missing) = frame

function inframes(::Type{T}, frame::Number, rate=missing) where T
    n = inframes(frame,rate)
    ismissing(n) && return missing
    T <: Integer ? trunc(T,n) : convert(T,n)
end
inframes(::Missing,fs=missing) = missing
inframes(::Type,::Missing,fs=missing) = missing

"""
    inHz(quantity)

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

If the given quantity is Unitful, we use the given units. If it is not we assume
it is already a value in Hz.

## Examples

```jldoctest
julia> inHz(1.0kHz)
1000.0
```

"""
inHz(x::Quantity) = ustrip(uconvert(Hz, x))
inHz(x::Number) = x
inHz(::Missing) = missing
inHz(::Type{T},x) where T <: Integer = trunc(T,inHz(x))
inHz(::Type{T},x) where T = convert(T,inHz(x))
inHz(::Type,x::Missing) = missing

"""
    inseconds(quantity[, rate])

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

If the given quantity is Unitful, we use the given units. If it is not we assume
it is already a value in seconds.

For some units (e.g. frames) you will need to specify a frame rate. If not
specified the rate is `missing`.

## Examples

```jldoctest
julia> inseconds(50.0ms)
0.05

julia> inseconds(441frames, 44100Hz)
0.01
```

"""
inseconds(x::Unitful.Time, rate=missing) = ustrip(uconvert(s,x))
inseconds(x::FrameQuant, rate=missing) = inframes(x,rate) / inHz(rate)
function inseconds(x::Quantity, rate=missing)
    error("Don't know how to convert $x to seconds.")
end
inseconds(x::Number, rate=missing) = x
function inseconds(::Type{T},x::Number,rate=missing) where T
    n = inseconds(x,rate)
    ismissing(n) && return missing
    T <: Integer ? trunc(T,n) : convert(T,n)
    convert(T,inseconds(x,rate))
end
inseconds(::Missing,r=missing) = missing
inseconds(::Type,::Missing,r=missing) = missing

maybeseconds(n::Number) = n*s
maybeseconds(n::Quantity) = n


text:
using Javis#, Animations

nframes = 200

function ground(args...)
    background("white") # canvas background
    sethue("black") # pen color
end

myvideo = Video(300, 300)
Background(1:200, ground)

function draw_line(p1 = O, p2 = O, color = "black", action = :stroke, edge = "solid")
    sethue(color)
    setdash(edge)
    line(p1, p2, action)
end

# Draw the coordinate lines with solid in the positive side and dashed in the negative one
vert_line = Object(
    (args...) -> draw_line(Point(0, 170), Point(0, 0), "black", :stroke, "longdashed"),
)
horiz_line = Object(
    (args...) -> draw_line(Point(-170, 0), Point(0, 0), "black", :stroke, "longdashed"),
)

vert_line_solid =
    Object((args...) -> draw_line(Point(0, 0), Point(0, -170), "black", :stroke, "solid"))
horiz_line_solid =
    Object((args...) -> draw_line(Point(0, 0), Point(170, 0), "black", :stroke, "solid"))

"""
Take the position of an object as Point, and draw
a squircle that is its projection in the positive orthant
"""
function draw_proj(pos, color = "blue")
    pos_proj = Point(max(pos.x, 0), min(pos.y, 0))
    sethue(color)
    squircle(pos_proj, 10, 10, :fill)
    return pos_proj
end

function object(p = O, color = "black", radius = 12)
    sethue(color)
    circle(p, radius, :fill)
    return p
end

red_ball = Object(1:nframes, (args...) -> object(O, "red"), Point(100, 0))
act!(red_ball, Action(anim_rotate_around(2π, O)))

Object(1:nframes, (args...) -> draw_proj(pos(red_ball)))

render(myvideo; pathname = joinpath(@__DIR__, "gifs/circle_projection.gif"), framerate = 60)


text:
module EnclosingWrapper

using JavaCall

const Class = @jimport java.lang.Class
const Encloser = @jimport org.hipparchus.geometry.enclosing.Encloser
const EnclosingBall = @jimport org.hipparchus.geometry.enclosing.EnclosingBall
const Iterable = @jimport java.lang.Iterable
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Object = @jimport java.lang.Object
const Point = @jimport org.hipparchus.geometry.Point
const SupportBallGenerator = @jimport org.hipparchus.geometry.enclosing.SupportBallGenerator
const WelzlEncloser = @jimport org.hipparchus.geometry.enclosing.WelzlEncloser

export Encloser
export EnclosingBall
export SupportBallGenerator
export WelzlEncloser
export ball_on_support
export contains
export enclose
export equals
export get_center
export get_class
export get_radius
export get_support
export get_support_size
export hash_code
export notify
export notify_all
export select_farthest
export to_string
export wait

include("encloser.jl")
include("enclosing_ball.jl")
include("support_ball_generator.jl")
include("welzl_encloser.jl")

end


text:
safeint(::Type{Bool}) = Bool
baseint(::Type{Bool}) = Bool
safeint(x::Bool) = x
baseint(x::Bool) = x

for (S,I) in (
    (:SafeInt8, :Int8), (:SafeInt16, :Int16), (:SafeInt32, :Int32), (:SafeInt64, :Int64), (:SafeInt128, :Int128),
    (:SafeUInt8, :UInt8), (:SafeUInt16, :UInt16), (:SafeUInt32, :UInt32), (:SafeUInt64, :UInt64), (:SafeUInt128, :UInt128) )
  @eval begin
    @inline safeint(::Type{$I}) = $S
    @inline safeint(::Type{$S}) = $S
    @inline baseint(::Type{$I}) = $I
    @inline baseint(::Type{$S}) = $I

    @inline safeint(x::$I) = reinterpret($S, x)
    @inline baseint(x::$S) = reinterpret($I, x)
    @inline safeint(x::$S) = x
    @inline baseint(x::$I) = x

    @inline $S(x::$I) = reinterpret($S, x)
    @inline $I(x::$S) = reinterpret($I, x)

    @inline $S(x::Bool) = ifelse(x, one($S), zero($S))

    @inline SafeInteger(::Type{$I}) = $S
    @inline SafeInteger(x::$I) = reinterpret($S, x)
  end
end

for (S,I) in (
    (:SafeInt8, :Int8), (:SafeInt16, :Int16), (:SafeInt32, :Int32), (:SafeInt64, :Int64), (:SafeInt128, :Int128) )
  @eval begin
    @inline SafeSigned(::Type{$I}) = $S
    @inline SafeSigned(x::$I) = reinterpret($S, x)
  end
end

for (S,I) in (
    (:SafeUInt8, :UInt8), (:SafeUInt16, :UInt16), (:SafeUInt32, :UInt32), (:SafeUInt64, :UInt64), (:SafeUInt128, :UInt128) )
  @eval begin
    @inline SafeUnsigned(::Type{$I}) = $S
    @inline SafeUnsigned(x::$I) = reinterpret($S, x)
  end
end

SafeInt128(x::T) where T<:Union{SafeInt8,SafeInt16,SafeInt32,SafeInt64} =
    safeint(Int128(baseint(x)))
SafeInt64(x::T) where T<:Union{SafeInt8,SafeInt16,SafeInt32} =
    safeint(Int64(baseint(x)))
SafeInt32(x::T) where T<:Union{SafeInt8,SafeInt16} =
    safeint(Int32(baseint(x)))
SafeInt16(x::T) where T<:SafeInt8 =
    safeint(Int16(baseint(x)))

SafeUInt128(x::T) where T<:Union{SafeUInt8,SafeUInt16,SafeUInt32,SafeUInt64} =
    safeint(UInt128(baseint(x)))
SafeUInt64(x::T) where T<:Union{SafeUInt8,SafeUInt16,SafeUInt32} =
    safeint(UInt64(baseint(x)))
SafeUInt32(x::T) where T<:Union{SafeUInt8,SafeUInt16} =
    safeint(UInt32(baseint(x)))
SafeUInt16(x::T) where T<:SafeUInt8 =
    safeint(UInt16(baseint(x)))

for (SS,SU, IS, IU) in (
    (:SafeInt8, :SafeUInt8, :Int8, :UInt8),
    (:SafeInt16, :SafeUInt16, :Int16, :UInt16),
    (:SafeInt32, :SafeUInt32, :Int32, :UInt32),
    (:SafeInt64, :SafeUInt64, :Int64, :UInt64),
    (:SafeInt128, :SafeUInt128, :Int128, :UInt128) )
   @eval begin
     $SS(x::T) where T<:Signed   = safeint(baseint($SS)(x))
     $SS(x::T) where T<:Unsigned = safeint(baseint($SS)(baseint($SU)(x)))
     $SU(x::T) where T<:Signed   = safeint(baseint($SU)(x))
     $SU(x::T) where T<:Unsigned = safeint(baseint($SU)(baseint($SS)(x)))
     $SS(x::T) where T<:SafeSigned   = safeint(baseint($SS)(x))
     $SS(x::T) where T<:SafeUnsigned = safeint(baseint($SS)(baseint($SU)(x)))
     $SU(x::T) where T<:SafeSigned   = safeint(baseint($SU)(x))
     $SU(x::T) where T<:SafeUnsigned = safeint(baseint($SU)(baseint($SS)(x)))

     $IS(x::T) where T<:SafeSigned   = $IS(baseint(x))
     $IS(x::T) where T<:SafeUnsigned = $IS(baseint(x))
     $IU(x::T) where T<:SafeUnsigned = $IU(baseint(x))
     $IU(x::T) where T<:SafeSigned   = $IU(baseint(x))
  end
end

safeint(x::UnitRange{I}) where {I<:Integer}     = safeint(x.start):safeint(x.stop)
baseint(x::UnitRange{S}) where {S<:SafeInteger} = baseint(x.start):baseint(x.stop)

function Rational{T}(num::I, den::I) where {I<:Integer, T<:SafeInteger}
    num == den == zero(T) && __throw_rational_argerror(T)
    num2, den2 = (sign(den) < 0) ? (den !== typemin(I) ? divgcd(-num, -den) : divgcd(num, den) ) : divgcd(num, den)
    Rational{T}(SaferIntegers.safeint(num2) , SaferIntegers.baseint(den2))
end

for (S,I) in (
    (:SafeInt8, :Int8), (:SafeInt16, :Int16), (:SafeInt32, :Int32), (:SafeInt64, :Int64), (:SafeInt128, :Int128),
    (:SafeUInt8, :UInt8), (:SafeUInt16, :UInt16), (:SafeUInt32, :UInt32), (:SafeUInt64, :UInt64), (:SafeUInt128, :UInt128) )
  @eval begin
    Base.Float64(x::$S) = Float64(reinterpret($I, x))
    Base.Float32(x::$S) = Float32(reinterpret($I, x))
    Base.Float16(x::$S) = Float16(reinterpret($I, x))
    Base.BigFloat(x::$S) = BigFloat(reinterpret($I, x))
    Base.BigInt(x::$S) = BigInt(reinterpret($I, x))
  end
end

for (S,I) in (
    (:SafeInt8, :Int8), (:SafeInt16, :Int16), (:SafeInt32, :Int32), (:SafeInt64, :Int64), (:SafeInt128, :Int128) )
    for F in (:Float64, :Float32, :Float16)
        @eval function $S(x::$F)
            !isinteger(x) && throw(DomainError("$x"))
            typemin($I) < x <typemax($I) ? $S($I(x)) : throw(OverflowError("$x"))
        end
    end
end

for (S,I) in (
    (:SafeUInt8, :UInt8), (:SafeUInt16, :UInt16), (:SafeUInt32, :UInt32), (:SafeUInt64, :UInt64), (:SafeUInt128, :UInt128) )
    for F in (:Float64, :Float32, :Float16)
        @eval function $S(x::$F)
            !isinteger(x) && throw(DomainError("$x"))
            typemin($I) < x <typemax($I) ? $S($I(x)) : throw(OverflowError("$x"))
        end
    end
end

for T in (:SafeInt8, :SafeInt16, :SafeInt32, :SafeInt64, :SafeInt128,
          :SafeUInt8, :SafeUInt16, :SafeUInt32, :SafeUInt64, :SafeUInt128)
    @eval $T(x::$T) = x
end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

import Pkg, Libdl, PackageCompiler

function compile(sysimage_path = "JuliaSysimage.$(Libdl.dlext)")
    env_to_precompile = joinpath(@__DIR__, "..")
    precompile_execution_file = joinpath(@__DIR__, "precompile.jl")
    project_filename = joinpath(env_to_precompile, "Project.toml")
    project = Pkg.API.read_project(project_filename)
    used_packages = Symbol.(collect(keys(project.deps)))
    # don't need these ever after building this so no need to have them in the sysimage
    filter!(x -> x ∉ [:Libdl, :PackageCompiler, :Pkg], used_packages)
    if Libdl.dlext == "dll"
        @warn "Ignoring packages which use gl dlls on Windows as these cause build errors"
        # see https://github.com/JuliaLang/PackageCompiler.jl/issues/365
        # recently the build leaves a corrupt dll if we include these rather than the error in the issue above
        used_packages = filter(x -> x ∉ [:Makie, :ImageView], used_packages)
    end
    @info "Building a custom sysimage for OpticSim.jl."
    PackageCompiler.create_sysimage(used_packages, sysimage_path = sysimage_path, project = env_to_precompile, precompile_execution_file = precompile_execution_file)
end


text:
using MDP

const rate        = [0.3 0.7 0.9]
const arrivalRate = 0.5

const serviceCost = [6 10 20]
const holdingCost = 1
const dropPenalty = 100_000

const M = 100
const A = size(rate,2)

P = [ spzeros(Float64, M+1, M+1) for u = 1:A]
C = zeros(Float64, M+1, A)

# Initialize cost matrix
C[1,:] = 0

for x = 2:M
    for u = 1:A
        C[x,u] = (x-1) * holdingCost + serviceCost[u]
    end
end

# Add expected cost for dopping packets
for u = 1:A
    C[M+1,u] = M * holdingCost + serviceCost[u] + arrivalRate * dropPenalty
end

# Initialize Probability matrix
for x = 2:M
    for u = 1:A
        P[u][x, x-1] = (1 - arrivalRate) * rate[u]
        P[u][x, x]   = (1 - arrivalRate) * (1 - rate[u]) + arrivalRate * rate[u]
        P[u][x, x+1] = arrivalRate * ( 1 - rate[u])
    end
end

for u = 1:A
    P[u][1,1] = (1 - arrivalRate)
    P[u][1,2] = arrivalRate

    P[u][M+1, M+1] = (1 - rate[u])
    P[u][M+1, M  ] = rate[u]
end

model = ProbModel(C, P; objective=:Min)

@time (v,g) = valueIteration(model; discount=0.9)

using StatsBase

println(rle(g))


text:
function FieldPVCoordinates(arg0::Field, arg1::PVCoordinates)
    return FieldPVCoordinates((Field, PVCoordinates), arg0, arg1)
end

function FieldPVCoordinates(arg0::FieldPVCoordinates, arg1::FieldPVCoordinates)
    return FieldPVCoordinates((FieldPVCoordinates, FieldPVCoordinates), arg0, arg1)
end

function FieldPVCoordinates(arg0::FieldVector3D)
    return FieldPVCoordinates((FieldVector3D,), arg0)
end

function FieldPVCoordinates(arg0::FieldVector3D, arg1::FieldVector3D)
    return FieldPVCoordinates((FieldVector3D, FieldVector3D), arg0, arg1)
end

function FieldPVCoordinates(arg0::FieldVector3D, arg1::FieldVector3D, arg2::FieldVector3D)
    return FieldPVCoordinates((FieldVector3D, FieldVector3D, FieldVector3D), arg0, arg1, arg2)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::FieldPVCoordinates)
    return FieldPVCoordinates((RealFieldElement, FieldPVCoordinates), arg0, arg1)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::FieldPVCoordinates, arg2::RealFieldElement, arg3::FieldPVCoordinates)
    return FieldPVCoordinates((RealFieldElement, FieldPVCoordinates, RealFieldElement, FieldPVCoordinates), arg0, arg1, arg2, arg3)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::FieldPVCoordinates, arg2::RealFieldElement, arg3::FieldPVCoordinates, arg4::RealFieldElement, arg5::FieldPVCoordinates)
    return FieldPVCoordinates((RealFieldElement, FieldPVCoordinates, RealFieldElement, FieldPVCoordinates, RealFieldElement, FieldPVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::FieldPVCoordinates, arg2::RealFieldElement, arg3::FieldPVCoordinates, arg4::RealFieldElement, arg5::FieldPVCoordinates, arg6::RealFieldElement, arg7::FieldPVCoordinates)
    return FieldPVCoordinates((RealFieldElement, FieldPVCoordinates, RealFieldElement, FieldPVCoordinates, RealFieldElement, FieldPVCoordinates, RealFieldElement, FieldPVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::PVCoordinates)
    return FieldPVCoordinates((RealFieldElement, PVCoordinates), arg0, arg1)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::PVCoordinates, arg2::RealFieldElement, arg3::PVCoordinates)
    return FieldPVCoordinates((RealFieldElement, PVCoordinates, RealFieldElement, PVCoordinates), arg0, arg1, arg2, arg3)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::PVCoordinates, arg2::RealFieldElement, arg3::PVCoordinates, arg4::RealFieldElement, arg5::PVCoordinates)
    return FieldPVCoordinates((RealFieldElement, PVCoordinates, RealFieldElement, PVCoordinates, RealFieldElement, PVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldPVCoordinates(arg0::RealFieldElement, arg1::PVCoordinates, arg2::RealFieldElement, arg3::PVCoordinates, arg4::RealFieldElement, arg5::PVCoordinates, arg6::RealFieldElement, arg7::PVCoordinates)
    return FieldPVCoordinates((RealFieldElement, PVCoordinates, RealFieldElement, PVCoordinates, RealFieldElement, PVCoordinates, RealFieldElement, PVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function FieldPVCoordinates(arg0::jdouble, arg1::FieldPVCoordinates)
    return FieldPVCoordinates((jdouble, FieldPVCoordinates), arg0, arg1)
end

function FieldPVCoordinates(arg0::jdouble, arg1::FieldPVCoordinates, arg2::jdouble, arg3::FieldPVCoordinates)
    return FieldPVCoordinates((jdouble, FieldPVCoordinates, jdouble, FieldPVCoordinates), arg0, arg1, arg2, arg3)
end

function FieldPVCoordinates(arg0::jdouble, arg1::FieldPVCoordinates, arg2::jdouble, arg3::FieldPVCoordinates, arg4::jdouble, arg5::FieldPVCoordinates)
    return FieldPVCoordinates((jdouble, FieldPVCoordinates, jdouble, FieldPVCoordinates, jdouble, FieldPVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldPVCoordinates(arg0::jdouble, arg1::FieldPVCoordinates, arg2::jdouble, arg3::FieldPVCoordinates, arg4::jdouble, arg5::FieldPVCoordinates, arg6::jdouble, arg7::FieldPVCoordinates)
    return FieldPVCoordinates((jdouble, FieldPVCoordinates, jdouble, FieldPVCoordinates, jdouble, FieldPVCoordinates, jdouble, FieldPVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function shifted_by(obj::FieldPVCoordinates, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldPVCoordinates, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldPVCoordinates, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldPVCoordinates, (jdouble,), arg0)
end

function to_string(obj::FieldPVCoordinates)
    return jcall(obj, "toString", JString, ())
end


text:
@testset "Utility Functions" begin
    @testset "Dimension Checking" begin
        F = reshape([1., 0.], 1, 2)
        G = [1. 1.; 0. 1.]
        V = reshape([1.], 1, 1)
        W = [1. 0.; 0. 0.2]

        @test_throws TypeError check_dimensions(F, G, V=V, W=W)

        F = reshape([1., 0.], 1, 2)
        G = [1. 1.; 0. 1.]
        V = Symmetric(reshape([1.], 1, 1))
        W = Symmetric([1. 0.; 0. 0.2])

        @test check_dimensions(F, G, V=V, W=W) == (1, 2)
        @test check_dimensions(F, G) == (1, 2)

        F = reshape([1., 0.], 1, 2)
        G = [1. 1.; 0. 1.]
        V = Symmetric(reshape([1.], 1, 1))
        W = Symmetric([1. 0. 0.; 0. 2. 0.; 0. 0. 0.2])

        @test_throws DimensionMismatch check_dimensions(F, G, V=V, W=W)

        F = reshape([1., 0.], 1, 2)
        G = [1. 1.; 0. 1.]
        V = Symmetric([1. 0.; 0. 1.])
        W = Symmetric([1. 0.; 0. 2.])

        @test_throws DimensionMismatch check_dimensions(F, G, V=V, W=W)

        F = reshape([1., 0.], 1, 2)
        G = [1. 1. 1.; 0. 1. 1.; 0. 0. 0.]
        V = Symmetric(reshape([1.], 1, 1))
        W = Symmetric([1. 0.; 0. 2.])

        @test_throws DimensionMismatch check_dimensions(F, G, V=V, W=W)

        # Cover exception for when G is not square
        F = reshape([1., 0.], 1, 2)
        G = [1. 1. 1.; 0. 1. 0.]
        V = Symmetric(reshape([1.], 1, 1))
        W = Symmetric([1. 0.; 0. 0.2])
        @test_throws DimensionMismatch check_dimensions(F, G, V=V, W=W)

        G = [1. 1.; 0. 1.]
        y = reshape([1., 1., 1., 1.], 4, 1)

        @test check_dimensions(F, G, V=V, W=W, Y=y) == (1, 2)

        y = [[1.], [1.], [1.], [1.]]

        @test_throws TypeError check_dimensions(F, G, V=V, W=W, Y=y) == (1, 2)
        @test_throws TypeError check_dimensions(F, G, Y=y) == (1,2)
    end

    @testset "Prior Computation" begin
        y_a = hcat(1.)
        F_a = hcat(1., 0.)
        m_a = ones(2)
        C_a = Symmetric(ones(2, 2))

        @test compute_prior(y_a, F_a, nothing, nothing) == (zeros(2), I(2))
        @test compute_prior(y_a, F_a, m_a, nothing) == (m_a, zeros(2, 2))
        @test compute_prior(y_a, F_a, nothing, C_a) == (zeros(2), ones(2, 2))
        @test compute_prior(y_a, F_a, m_a, C_a) == (m_a, C_a)

        # Make sure only Symmetric is allowed
        @test_throws MethodError compute_prior(y_a, F_a, nothing, ones(2, 2))

        y_b = hcat(2., 3.)
        F_b = [1. 0. .5; 0. 1. .5]
        m_b = [1., 1., .5]
        C_b = Symmetric(ones(3, 3))

        @test compute_prior(y_b, F_b, nothing, nothing) == (zeros(3), 9 * I(3))
        @test compute_prior(y_b, F_b, m_b, nothing) == (m_b, 3.0625 * I(3))
        @test compute_prior(y_b, F_b, nothing, C_b) == (zeros(3), ones(3, 3))
        @test compute_prior(y_b, F_b, m_b, C_b) == (m_b, C_b)
    end

    @testset "Plot Recipe" begin
        F = reshape([1., 0.], 1, 2)
        G = [1. 1.; 0. 1.]
        V = Symmetric(reshape([12.], 1, 1))
        W = Symmetric([1. 0.; 0. 0.2])
        θ₀ = zeros(2)
        T = 36
        θ, y = simulate(F, G, V, W, θ₀, T)
        a, R, m, C = kfilter(y, F, G, V, W)
        s, S = ksmoother(G, a, R, m, C)
        f, Q = fitted(F, V, s, S)
        fh, Qh = forecast(F, G, V, W, m[end,:], C[end], 10)

        # When applying the recipe, this function is required to check validity
        # of the recipe keys. To avoid importing Plots, extend it to always return
        # true.
        function RecipesBase.is_key_supported(x)
            return true
        end

        apply_recipe = RecipesBase.apply_recipe

        @test size(apply_recipe(Dict{Symbol,Any}(), DLMPlot(), y, f, Q
                               )[1].plotattributes[:color]) == (1,4)
        @test size(apply_recipe(Dict{Symbol,Any}(), DLMPlot(), y, f, Q, fh, Qh
                               )[1].plotattributes[:color]) == (1,7)

        @test_throws ArgumentError apply_recipe(Dict{Symbol,Any}(), DLMPlot(), y, f, Q, nothing, Qh)
        @test_throws ArgumentError apply_recipe(Dict{Symbol,Any}(), DLMPlot(), y, f, Q, fh, nothing)
    end
end


text:
struct SimulatableContext
    # map from original variable to popped variable
    popmap::Dict{Ghost.Variable, Ghost.Variable}
    # map from signatures to variable bindings
    opmap::Dict{Tuple, Ghost.Variable}
end
SimulatableContext() = SimulatableContext(Dict{Ghost.Variable, Vector{Ghost.Variable}}(),
                                          Dict{Tuple, Ghost.Variable}())

function Ghost.rebind_context!(tape::Ghost.Tape{SimulatableContext}, subs::Dict)
    replace!(tape.c.popmap) do kv
        old_var, new_var = kv
        old_var = get(subs, old_var, old_var)
        new_var = get(subs, new_var, new_var)

        return old_var => new_var
    end
    replace!(tape.c.opmap) do kv
        sig, var = kv
        newout = get(subs, var, var)
        newargs = [get(subs, arg, arg) for arg in Base.tail(sig)]

        return (sig[1], newargs...) => newout
    end
end

function _update_ctx!(ctx::SimulatableContext, vars)
    if _issimulatable(vars[1])
        if vars[1] isa Ghost.Input && length(vars) > 1
            ctx.popmap[Ghost.Variable(vars[1])] = Ghost.Variable(vars[2])
        elseif vars[1] isa Ghost.Call && _isbcast(vars[1].fn) && length(vars) > 1
            ctx.popmap[Ghost.Variable(vars[1])] = Ghost.Variable(vars[end - 1])
            ctx.popmap[Ghost.Variable(vars[2])] = Ghost.Variable(vars[end - 1])
        elseif vars[1] isa Ghost.Call && length(vars) > 1
            ctx.popmap[Ghost.Variable(vars[1])] = Ghost.Variable(vars[end - 1])
        # elseif vars[1] isa Ghost.Call # identity replacement
        #     ctx.popmap[Ghost.Variable(vars[1])] = vars[1].args[1]
        end
    end
    foreach(vars) do var
        (var isa Ghost.Call) && push!(ctx.opmap, (var.fn, var.args...) => Ghost.Variable(var))
    end
end

function _popcalls(ctx, args)
    calls = []
    new_args = []
    foreach(args) do arg
        if _gettapeval(arg) isa SBitstreamLike
            if haskey(ctx.popmap, arg)
                push!(new_args, ctx.popmap[arg])
            else
                popcall = Ghost.mkcall(getbit, arg)
                push!(calls, popcall)
                push!(new_args, Ghost.Variable(popcall))
            end
        else
            push!(new_args, arg)
        end
    end

    return calls, new_args
end

function _wrap_bcast_bits(popbits)
    wrap_calls = []
    wrap_bits = []
    for bit in popbits
        if _gettapeval(bit) isa SBit
            call = Ghost.mkcall(Ref, bit)
            push!(wrap_calls, call)
            push!(wrap_bits, Ghost.Variable(call))
        else
            push!(wrap_bits, bit)
        end
    end

    return wrap_calls, wrap_bits
end

function _unbroadcasted_transform(ctx, call, sim)
    # insert calls to pop bits from the args
    popcalls, popbits = _popcalls(ctx, call.args)
    # evaluate simulator on popped bits
    bit = Ghost.mkcall(sim, popbits...)
    # push resulting bits onto bitstream
    psh = Ghost.mkcall(setbit!, Ghost.Variable(call), Ghost.Variable(bit))

    return [call, popcalls..., bit, psh], 1
end

function _broadcasted_transform(ctx, call, sim)
    # ignore first (function) arg of broadcasted
    args = call.args[2:end]
    # insert calls to pop bits from the args
    popcalls, popbits = _popcalls(ctx, args)
    # materialize broadcasted
    mat = Ghost.mkcall(Base.materialize, Ghost.Variable(call))
    # evaluate simulator on popped bits
    bit = Ghost.mkcall(sim, popbits...)
    # push resulting bits onto bitstream
    psh = Ghost.mkcall(setbit!, Ghost.Variable(mat), Ghost.Variable(bit))

    return [call, mat, popcalls..., bit, psh], 1
end

function _broadcasted_transform(ctx, call, sims::AbstractArray)
    # ignore first (function) arg of broadcasted
    args = call.args[2:end]
    # insert calls to pop bits from the args
    popcalls, popbits = _popcalls(ctx, args)
    # materialize broadcasted
    mat = Ghost.mkcall(Base.materialize, Ghost.Variable(call))
    # evaluate simulators element-wise on popped bits
    wrapcalls, wrapbits = _wrap_bcast_bits(popbits)
    bits = Ghost.mkcall(Base.broadcasted, (f, a...) -> f(a...), sims, wrapbits...)
    matbits = Ghost.mkcall(Base.materialize, Ghost.Variable(bits))
    # push resulting bits onto bitstreams
    psh = Ghost.mkcall(setbit!, Ghost.Variable(mat), Ghost.Variable(matbits))

    return [call, mat, popcalls..., wrapcalls..., bits, matbits, psh], 1
end

_handle_bcast_and_transform(ctx, call, sim) =
    _isbcast(call.fn) ? _broadcasted_transform(ctx, call, sim) :
                        _unbroadcasted_transform(ctx, call, sim)

_simtransform(ctx, input::Ghost.Input) =
    _gettapeval(Ghost.Variable(input)) isa SBitstreamLike ?
        ([input, Ghost.mkcall(getbit, Ghost.Variable(input))], 1) :
        ([input], 1)
function _simtransform(ctx, call::Ghost.Call)
    # if the args don't contain SBitstreamLike, then skip
    sig = Ghost.get_type_parameters(Ghost.call_signature(call.fn, _gettapeval.(call.args)...))
    is_simulatable_primitive(sig...) || return [call], 1

    # get the simulator for this call signature
    sim = getsimulator(call.fn, map(arg -> _gettapeval(arg), call.args)...)

    # skip nosims
    isnothing(sim) && return [call], 1

    # if call has already been transformed,
    #  then delete this call and rebind to the transformed call
    haskey(ctx.opmap, (call.fn, call.args...)) && return [], ctx.opmap[(call.fn, call.args...)].id

    # otherwise, transform this call while handling broadcasting
    return _handle_bcast_and_transform(ctx, call, sim)
end

getbit(x) = x
getbit(x::SBitstreamLike) = pop!.(x)

setbit!(x::SBitstream, bit) = push!(x, bit)
setbit!(x::AbstractArray{<:SBitstream}, bits) = push!.(x, bits)

is_simulatable_primitive(sig...) = is_trace_primitive(sig...)

function simulator(f, args...)
    # if f itself is a primitive, do a manual tape
    if is_simulatable_primitive(Ghost.get_type_parameters(Ghost.call_signature(f, args...))...)
        tape = Ghost.Tape(SimulatableContext())
        inputs = Ghost.inputs!(tape, f, args...)
        if _isstruct(f)
            tape.result = push!(tape, Ghost.mkcall(inputs...))
        else
            tape.result = push!(tape, Ghost.mkcall(f, inputs[2:end]...))
        end
    else
        tape = trace(f, args...;
                     isprimitive = is_simulatable_primitive,
                     ctx = SimulatableContext())
    end
    transform!(_squash_binary_vararg, tape)
    transform!(_simtransform, _update_ctx!, tape)

    return tape
end
simulatable(f, args...) = Ghost.compile(simulator(f, args...))
show_simulatable(f, args...) = Ghost.to_expr(simulator(f, args...))

macro nosim(ex, kws...)
    kwargs = isempty(kws) ? false :
             (kws[1].args[1] == :kwargs) ? kws[1].args[2] :
             error("Unrecognized options: $kws")
    if @capture(ex, f_(args__))
        argtypes = map(args) do arg
            @capture(arg, x_::T_) ? T : :Any
        end
        primitive_sig = [:(::Type{<:$(esc(T))}) for T in argtypes]

        if kwargs
            return quote
                BitSAD.is_simulatable_primitive(::Type{Core.kwftype(typeof($(esc(f))))},
                                                ::Type{<:Any},
                                                ::Type{typeof($(esc(f)))},
                                                $(primitive_sig...)) = true
                BitSAD.getsimulator(::Core.kwftype(typeof($(esc(f)))),
                                    ::Any,
                                    ::typeof($(esc(f))),
                                    $(esc.(args)...)) = nothing
            end
        else
            return quote
                BitSAD.is_simulatable_primitive(::Type{typeof($(esc(f)))}, $(primitive_sig...)) = true
                BitSAD.getsimulator(::typeof($(esc(f))), $(esc.(args)...)) = nothing
            end
        end
    else
        error("Cannot parse @nosim $ex (expects @nosim f(arg1, arg2::T, ...).")
    end
end


text:
"""
    SeisKolmogoroff(in)

Transform a wavelet into its minimum phase equivalent.

# Arguments
- `in::Array{Real,1}`: input wavelet.

# Example
```julia
julia> using PyPlot
julia> w = Ricker()
julia> wmin = SeisKolmogoroff(w)
julia> plot(w); plot(wmin)
```

# Reference
* Claerbout, Jon F., 1976, Fundamentals of geophysical data processing.
McGraw-Hill Inc.
"""
function SeisKolmogoroff(w::Array{T,1}) where T<:Real

    nw = length(w)
    nf = 8*nextpow(2,nw)
    W = fft(cat(w,zeros(nf-nw),dims=1))
    A = log.(abs.(W) .+ 1.e-8)
    a = 2*ifft(A)
    n2 = floor(Int, nf/2)
    a[n2+2:nf] .= 0.0
    a[1] = a[1]/2
    A = exp.(fft(a))
    a = real(ifft(A))
    wmin = real(a[1:nw])

end