text:
module EigenBase

@enum Eigen4JuliaWrapperPrecondType begin
  EIGENBASE_PRECOND_NONE
  EIGENBASE_PRECOND_ILU
  EIGENBASE_PRECOND_ICC
end

@enum Eigen4JuliaWrapperPrecisionType begin
  EIGENBASE_FLOAT_SINGLE = 4
  EIGENBASE_FLOAT_DOUBLE = 8
end

@enum Eigen4JuliaWrapperCompressionType begin
  EIGENBASE_CRS_FORMAT = 1
  EIGENBASE_CCS_FORMAT = 2
end

mutable struct PrecondDescr
  # preconditioner handle
  pc_args::Ptr{Cvoid}
  # preconditioner type
  pctype::Eigen4JuliaWrapperPrecondType
  # problem matrix dimensions
  rows::Cint
  cols::Cint
  # precision of the preconditioner
  precision::Eigen4JuliaWrapperPrecisionType
  # complex flag
  is_complex::Cint
  # sparse format
  ctype::Eigen4JuliaWrapperCompressionType
  # nnz of the original matrix
  nnz::Cint
  # initialization flag
  is_init::Cint

  PrecondDescr() = new(C_NULL, EIGENBASE_PRECOND_NONE, 0,0,EIGENBASE_FLOAT_DOUBLE,0,EIGENBASE_CRS_FORMAT,0,0)

end #PrecondDescr

CreateEigenPreconditioner_qj(vals,innerIndices,outerIndexPtr,descr) = ccall((:CreatePreconditioner, "libeigenpc"), Cvoid, (Ptr{Cvoid},Ptr{Cint},Ptr{Cint},Ref{PrecondDescr}),vals,innerIndices,outerIndexPtr,descr)

ApplyEigenPreconditioner_qj(out, in, descr) = ccall((:ApplyPreconditioner, "libeigenpc"), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{PrecondDescr}),out, in, descr)

DestroyEigenPreconditioner_qj(descr) = ccall((:DestroyPreconditioner, "libeigenpc"), Cvoid, (Ref{PrecondDescr},),descr)

end #EigenBase

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
mutable struct Snailfish
    parent::Union{Snailfish, Nothing}
    depth::Int
    left::Union{Snailfish, Nothing}
    right::Union{Snailfish, Nothing}
    value::Union{Int, Nothing}

    Snailfish(parent, depth; left = nothing, right = nothing, value = nothing) = new(parent, depth, left, right, value)
end

Base.show(io::IO, node::Snailfish) = if (hasValue(node)) print(io, node.value) else print(io, "[", node.left, ",", node.right, "]") end

hasValue(node::Snailfish) = isa(node.value, Number)

splitCondition(node::Snailfish) = hasValue(node) && node.value >= 10
explodeCondition(node::Snailfish) = node.depth > 4 && !hasValue(node) && hasValue(node.left) && hasValue(node.right)

function setContents(node::Snailfish, left::Union{Snailfish, Nothing}, right::Union{Snailfish, Nothing}, value::Union{Int, Nothing})
    node.left, node.right, node.value = left, right, value
end

function parseSnailfish(input::String) :: Snailfish
    nesting = 1
    current = root = Snailfish(nothing, nesting)

    for token in collect(m.match for m = eachmatch(r"[\[\]\,\d]", input))
        if token == "["
            nesting = nesting + 1
            setContents(current, Snailfish(current, nesting), Snailfish(current, nesting), nothing)
            current = current.left
        elseif token == ","
            current = current.parent.right
        elseif token == "]"
            nesting = nesting - 1
            current = current != root ? current.parent : current
        else
            current.value = parse(Int, token)
        end
    end
    root
end

function getDepthFirstList(result::Array{Snailfish}, node::Snailfish) :: Array{Snailfish}
    push!(result, node)
    node.left != nothing && getDepthFirstList(result, node.left)
    node.right != nothing && getDepthFirstList(result, node.right)
    result
end

getDepthFirstList(node::Snailfish) :: Array{Snailfish} = getDepthFirstList(Array{Snailfish, 1}(), node)

function explode(root::Snailfish) :: Bool
    flat = getDepthFirstList(root)
    for i = 1 : length(flat)
        node = flat[i]
        !explodeCondition(node) && continue

        for nd in filter(hasValue, reverse(flat[1 : i - 1]))
            nd.value = nd.value + node.left.value
            break
        end
        for nd in filter(hasValue, flat[i + 3 : end])
            nd.value = nd.value + node.right.value
            break
        end
        setContents(node, nothing, nothing, 0)
        return true
    end
    false
end

function splitValue(root::Snailfish) :: Bool
    createChild(node::Snailfish, value::Int) = Snailfish(node, node.depth + 1; value = value)

    for node in filter(splitCondition, getDepthFirstList(root))
        setContents(node, createChild(node, floor(Int, node.value / 2)), createChild(node, ceil(Int, node.value / 2)), nothing)
        return true
    end
    false
end

function reduce(node::Snailfish)
    while true
        explode(node) || splitValue(node) || break
    end
end

function add(first::Snailfish, second::Snailfish) :: Snailfish
    for n in [getDepthFirstList(first); getDepthFirstList(second)]
        n.depth = n.depth + 1
    end
    first.parent = second.parent = Snailfish(nothing, 1; left = first, right = second)
end

function run()
    magnitude(node::Snailfish) = hasValue(node) ? node.value : magnitude(node.left) * 3 + magnitude(node.right) * 2

    lines = readlines("input.txt")

    current = parseSnailfish(lines[1])
    for line in lines[2 : end]
        current = add(current, parseSnailfish(line))
        reduce(current)
    end
    println("Part 1: ", magnitude(current))

    maxMagnitude = 0
    for a in lines, b in lines
        a == b && continue

        sum = add(parseSnailfish(a), parseSnailfish(b))
        reduce(sum)
        maxMagnitude = max(maxMagnitude, magnitude(sum))
    end
    println("Part 2: ", maxMagnitude)
end

run()

ext: jl

max_stars_repo_name: skrim/AdventOfCode

max_issues_repo_name: skrim/AdventOfCode

max_forks_repo_name: skrim/AdventOfCode


text:
module RxJulia

    include("./utils.jl")
    include("./kernel.jl")
    include("./operators.jl")

end # module


text:
@testset "Support function" begin
  t = Triangle(P2[(0,0), (1,0), (0,1)])
  @test supportfun(t, V2( 1,  0)) == P2(1, 0)
  @test supportfun(t, V2( 0,  1)) == P2(0, 1)
  @test supportfun(t, V2(-1, -1)) == P2(0, 0)
  @test supportfun(t, V2(-1,  1)) == P2(0, 1)

  b = Ball(P2(0,0), T(2))
  @test supportfun(b, V2( 1, 1)) ≈ P2(√2, √2)
  @test supportfun(b, V2( 1, 0)) ≈ P2(2, 0)
  @test supportfun(b, V2( 0, 1)) ≈ P2(0, 2)
  @test supportfun(b, V2(-1, 1)) ≈ P2(-√2, √2)
end


text:
using Turing
using AbstractGPs, Random

using LinearAlgebra
using VegaLite, DataFrames, StatsPlots, StatsBase
using RDatasets

Random.seed!(1789);

data = dataset("datasets", "iris")
species = data[!, "Species"]
index = shuffle(1:150)
# we extract the four measured quantities,
# so the dimension of the data is only d=4 for this toy example
dat = Matrix(data[index, 1:4])
labels = data[index, "Species"]

# non-linearize data to demonstrate ability of GPs to deal with non-linearity
dat[:, 1] = 0.5 * dat[:, 1] .^ 2 + 0.1 * dat[:, 1] .^ 3
dat[:, 2] = dat[:, 2] .^ 3 + 0.2 * dat[:, 2] .^ 4
dat[:, 3] = 0.1 * exp.(dat[:, 3]) - 0.2 * dat[:, 3] .^ 2
dat[:, 4] = 0.5 * log.(dat[:, 4]) .^ 2 + 0.01 * dat[:, 3] .^ 5

# normalize data
dt = fit(ZScoreTransform, dat; dims=1);
StatsBase.transform!(dt, dat);

@model function pPCA(x, ::Type{TV}=Array{Float64}) where {TV}
    # Dimensionality of the problem.
    N, D = size(x)
    # latent variable z
    z ~ filldist(Normal(), D, N)
    # weights/loadings W
    w ~ filldist(Normal(), D, D)
    mu = (w * z)'
    for d in 1:D
        x[:, d] ~ MvNormal(mu[:, d], I)
    end
end;

linear_kernel(α) = LinearKernel() ∘ ARDTransform(α)
sekernel(α, σ) = σ * SqExponentialKernel() ∘ ARDTransform(α);

@model function GPLVM_linear(Y, K=4, ::Type{T}=Float64) where {T}

    # Dimensionality of the problem.
    N, D = size(Y)
    # K is the dimension of the latent space
    @assert K <= D
    noise = 1e-3

    # Priors
    α ~ MvLogNormal(MvNormal(zeros(K), I))
    Z ~ filldist(Normal(), K, N)
    mu ~ filldist(Normal(), N)

    kernel = linear_kernel(α)

    gp = GP(mu, kernel)
    cv = cov(gp(ColVecs(Z), noise))
    return Y ~ filldist(MvNormal(mu, cv), D)
end;

@model function GPLVM(Y, K=4, ::Type{T}=Float64) where {T}

    # Dimensionality of the problem.
    N, D = size(Y)
    # K is the dimension of the latent space
    @assert K <= D
    noise = 1e-3

    # Priors
    α ~ MvLogNormal(MvNormal(zeros(K), I))
    σ ~ LogNormal(0.0, 1.0)
    Z ~ filldist(Normal(), K, N)
    mu ~ filldist(Normal(), N)

    kernel = sekernel(α, σ)

    gp = GP(mu, kernel)
    cv = cov(gp(ColVecs(Z), noise))
    return Y ~ filldist(MvNormal(mu, cv), D)
end;

# Standard GPs don't scale very well in n, so we use a small subsample for the purpose of this tutorial
n_data = 40
# number of features to use from dataset
n_features = 4
# latent dimension for GP case
ndim = 4;

ppca = pPCA(dat[1:n_data, 1:n_features])
chain_ppca = sample(ppca, NUTS(), 1000);

# we extract the posterior mean estimates of the parameters from the chain
w = reshape(mean(group(chain_ppca, :w))[:, 2], (n_features, n_features))
z = reshape(mean(group(chain_ppca, :z))[:, 2], (n_features, n_data))
X = w * z

df_pre = DataFrame(z', :auto)
rename!(df_pre, Symbol.(["z" * string(i) for i in collect(1:n_features)]))
df_pre[!, :type] = labels[1:n_data]
p_ppca = @vlplot(:point, x = :z1, y = :z2, color = "type:n")(df_pre)

gplvm_linear = GPLVM_linear(dat[1:n_data, 1:n_features], ndim)

chain_linear = sample(gplvm_linear, NUTS(), 500)
# we extract the posterior mean estimates of the parameters from the chain
z_mean = reshape(mean(group(chain_linear, :Z))[:, 2], (n_features, n_data))
alpha_mean = mean(group(chain_linear, :α))[:, 2]

df_gplvm_linear = DataFrame(z_mean', :auto)
rename!(df_gplvm_linear, Symbol.(["z" * string(i) for i in collect(1:ndim)]))
df_gplvm_linear[!, :sample] = 1:n_data
df_gplvm_linear[!, :labels] = labels[1:n_data]
alpha_indices = sortperm(alpha_mean; rev=true)[1:2]
println(alpha_indices)
df_gplvm_linear[!, :ard1] = z_mean[alpha_indices[1], :]
df_gplvm_linear[!, :ard2] = z_mean[alpha_indices[2], :]

p_linear = @vlplot(:point, x = :ard1, y = :ard2, color = "labels:n")(df_gplvm_linear)
p_linear

gplvm = GPLVM(dat[1:n_data, 1:n_features], ndim)

chain_gplvm = sample(gplvm, NUTS(), 500)
# we extract the posterior mean estimates of the parameters from the chain
z_mean = reshape(mean(group(chain_gplvm, :Z))[:, 2], (ndim, n_data))
alpha_mean = mean(group(chain_gplvm, :α))[:, 2]

df_gplvm = DataFrame(z_mean', :auto)
rename!(df_gplvm, Symbol.(["z" * string(i) for i in collect(1:ndim)]))
df_gplvm[!, :sample] = 1:n_data
df_gplvm[!, :labels] = labels[1:n_data]
alpha_indices = sortperm(alpha_mean; rev=true)[1:2]
println(alpha_indices)
df_gplvm[!, :ard1] = z_mean[alpha_indices[1], :]
df_gplvm[!, :ard2] = z_mean[alpha_indices[2], :]

p_gplvm = @vlplot(:point, x = :ard1, y = :ard2, color = "labels:n")(df_gplvm)
p_gplvm

let
    @assert abs(
        mean(z_mean[alpha_indices[1], labels[1:n_data] .== "setosa"]) -
        mean(z_mean[alpha_indices[1], labels[1:n_data] .!= "setosa"]),
    ) > 1.4
end

using Stheno
@model function GPLVM_sparse(Y, K, ::Type{T}=Float64) where {T}

    # Dimensionality of the problem.
    N, D = size(Y)
    # dimension of latent space
    @assert K <= D
    # number of inducing points
    n_inducing = 25
    noise = 1e-3

    # Priors
    α ~ MvLogNormal(MvNormal(zeros(K), I))
    σ ~ LogNormal(1.0, 1.0)
    Z ~ filldist(Normal(), K, N)
    mu ~ filldist(Normal(), N)

    kernel = σ * SqExponentialKernel() ∘ ARDTransform(α)

    ## Standard
    # gpc = GPC()
    # f = atomic(GP(kernel), gpc)
    # gp = f(ColVecs(Z), noise)
    # Y ~ filldist(gp, D)

    ## SPARSE GP
    #  xu = reshape(repeat(locations, K), :, K) # inducing points
    #  xu = reshape(repeat(collect(range(-2.0, 2.0; length=20)), K), :, K) # inducing points
    lbound = minimum(Y) + 1e-6
    ubound = maximum(Y) - 1e-6
    #  locations ~ filldist(Uniform(lbound, ubound), n_inducing)
    #  locations = [-2., -1.5 -1., -0.5, -0.25, 0.25, 0.5, 1., 2.]
    #  locations = collect(LinRange(lbound, ubound, n_inducing))
    locations = quantile(vec(Y), LinRange(0.01, 0.99, n_inducing))
    xu = reshape(locations, 1, :)
    gp = atomic(GP(kernel), GPC())
    fobs = gp(ColVecs(Z), noise)
    finducing = gp(xu, 1e-12)
    sfgp = SparseFiniteGP(fobs, finducing)
    cv = cov(sfgp.fobs)
    return Y ~ filldist(MvNormal(mu, cv), D)
end

n_data = 50
gplvm_sparse = GPLVM_sparse(dat[1:n_data, :], ndim)

chain_gplvm_sparse = sample(gplvm_sparse, NUTS(), 500)
# we extract the posterior mean estimates of the parameters from the chain
z_mean = reshape(mean(group(chain_gplvm_sparse, :Z))[:, 2], (ndim, n_data))
alpha_mean = mean(group(chain_gplvm_sparse, :α))[:, 2]

df_gplvm_sparse = DataFrame(z_mean', :auto)
rename!(df_gplvm_sparse, Symbol.(["z" * string(i) for i in collect(1:ndim)]))
df_gplvm_sparse[!, :sample] = 1:n_data
df_gplvm_sparse[!, :labels] = labels[1:n_data]
alpha_indices = sortperm(alpha_mean; rev=true)[1:2]
df_gplvm_sparse[!, :ard1] = z_mean[alpha_indices[1], :]
df_gplvm_sparse[!, :ard2] = z_mean[alpha_indices[2], :]
p_sparse = @vlplot(:point, x = :ard1, y = :ard2, color = "labels:n")(df_gplvm_sparse)
p_sparse

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
using DashBootstrapComponents

popover = html_div([
    dbc_button(
        "Toggle",
        id = "toggle",
        color = "success",
        className = "me-4",
        n_clicks = 0,
    ),
    dbc_button("Target", id = "target", color = "danger", n_clicks = 0),
    dbc_popover(
        [
            dbc_popoverheader("Popover header"),
            dbc_popoverbody("And here's some amazing content. Cool!"),
        ],
        id = "popover",
        is_open = false,
        target = "target",
    ),
]);

callback!(
    app,
    Output("popover", "is_open"),
    Input("toggle", "n_clicks"),
    State("popover", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open == 1
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using Documenter, CUDAdrv

const src = "https://github.com/JuliaGPU/CUDAdrv.jl"
const dst = "https://juliagpu.gitlab.io/CUDAdrv.jl/"

function main()
    makedocs(
        sitename = "CUDAdrv.jl",
        authors = "Tim Besard",
        repo = "$src/blob/{commit}{path}#{line}",
        format = Documenter.HTML(
            # Use clean URLs on CI
            prettyurls = get(ENV, "CI", nothing) == "true",
            canonical = dst,
            assets = ["assets/favicon.ico"],
            analytics = "UA-154489943-5",
        ),
        doctest = false,
        pages = Any[
            "Home"    => "index.md",
            "APIs" => [
                "driver.md",
            ]
        ]
    )
end

isinteractive() || main()


text:
using AxisArrays
using MetaArrays
using ImageInTerminal
using Colors

struct AxisMeta{Axs}
    axes::Axs # map from :axis -> metadata
end
AxisMeta(;kwds...) = AxisMeta(kwds.data)
Base.getproperty(x::AxisMeta,field::Symbol) =
  getproperty(getfield(x,:axes),field)

function validate(meta::AxisMeta,data::AxisArray)
    for field in fieldnames(typeof(Base.getfield(meta,:axes)))
        if field ∉ axisnames(data)
            error("Extra field `$field` in axis metadata.")
        end
    end

    for axis in axisnames(data)
        if axis ∉ fieldnames(typeof(Base.getfield(meta,:axes)))
            error("Missing field `$axis` in axis metadata.")
        end
    end
end

axismeta(data::AxisArray;keys...) = axismeta(keys.data,data)
axismeta(meta::NamedTuple,data::AxisArray) = axismeta(AxisMeta(meta),data)
function axismeta(meta::AxisMeta,data::AxisArray)
    validate(meta,data)
    MetaArray(meta,data)
end

const MetaAxisArray{A,M} = MetaArray{A,M} where {A<:AxisArray,M<:AxisMeta}
const MetaAxisLike = Union{AxisMeta,MetaAxisArray}
function MetaAxisArray(meta::AxisMeta,data::AxisArray)
    validate(meta,data)
    MetaArray(meta,data)
end

function addaxes(meta::AxisMeta;kwds...)
    AxisMeta(merge(Base.getfield(meta,:axes),kwds.data))
end

function removeaxes(meta::AxisMeta,rem...)
  axes = Base.getfield(meta,:axes)
  fields = setdiff(fieldnames(typeof(axes)),rem)
  vals = Tuple(map(f -> axes[f],fields))
  AxisMeta(NamedTuple{Tuple(fields),typeof(vals)}(vals))
end

function describe_axes(io::IO,x)
  for ax in AxisArrays.axes(x)
    println(io,string(AxisArrays.axisname(ax))," ",
            string(round(ustrip(ax.val[1]),digits=2)),
            " - ",string(round(ustrip(ax.val[end]),digits=2)),
            string(unit(ax.val[1])))
  end
end

resultname(x) = "Data with axes: "
function Base.show(io::IO,mime::MIME"text/plain",x::MetaAxisArray)
  if get(io, :compact, false)
    println(io,resultname(x))
  else
    if hastimes(x) isa HasTimes
      println(io,string(duration(x))," ",resultname(x))
    else
      println(io,resultname(x))
    end
    describe_axes(io,x)
    if ndims(x) <= 2
      if eltype(x) <: Complex
        @warn "Ignoring phase in display of complex values."
        x = real.(abs.(x))
      end
      lo,hi = extrema(x)
      x .= (x .- lo) ./ (hi - lo)
      if ndims(x) == 2
        image = reverse(Gray.(getcontents(x)),dims=2)'
      else
        image = Gray.(getcontents(x))
      end
      ImageInTerminal.imshow(io,image,ImageInTerminal.colormode[1])
    end
  end
end

Base.show(io::IO, m::MIME"img/png", x::MetaAxisArray) = show(io,m,plotaxes(x))


text:
function setticks(y::Real)
    if y == 10
        return 0:3:10
    elseif 10 < y <= 20
        return 0:5:y
    elseif 20 < y <= 50
        return 0:10:y
    else
        s = div(y, 4)
        d, r = divrem(s, 10)
        return r < 5 ? (0:10d:y) : (0:((d + 1) * 10):y)
    end
end

function calculateld!(
    gwas::DataFrame,
    ref::SnpData;
    snp::Union{AbstractString, Tuple{AbstractString, Int}} = "index"
)

    gwas.ind = findsnps(gwas, ref; matchalleles = false)
    dropmissing!(gwas, "ind")
    n = size(gwas, 1)
    gwas.LD = fill(0.0, n)
    if snp == "index"
        i = argmax(gwas.P)
        snp = gwas.SNP[i]
        geno = convert(Matrix{Float64}, ref.snparray[:, gwas.ind])
        for j in 1:n
            gwas.LD[j] = cor(geno[:, i], geno[:, j])^2
        end
        gwas.index = fill(snp, n)
        return
    else
        if snp isa AbstractString
            (chr, bp) = getsnpinfo(snp, ref)
        else
            (chr, bp) = snp
        end
        i = findfirst((gwas.CHR .== chr) .& (gwas.BP .== bp))
        if isnothing(i)
            snp = ""
            i = findfirst((ref.snp_info.chromosome .== chr) .& (ref.snp_info.position .== bp))
            geno = convert(Matrix{Float64}, ref.snparray[:, [i; gwas.ind]])
            for j in 1:n
                gwas.LD[j] = cor(geno[:, 1], geno[:, j + 1])^2
            end
        else
            snp = gwas.SNP[i]
            geno = convert(Matrix{Float64}, ref.snparray[:, gwas.ind])
            for j in 1:n
                gwas.LD[j] = cor(geno[:, i], geno[:, j])^2
            end
        end
        gwas.index = fill(snp, n)
        return
    end
end

"""
    plotlocus!(ax::Axis, chromosome::AbstractString, range1::Real, range2::Real, gwas::DataFrame; kwargs)
    plotlocus!(ax::Axis, chromosome::AbstractString, bp::Real, gwas::DataFrame; kwargs)
    plotlocus!(ax::Axis, gene::AbstractString, gwas::DataFrame, gencode::DataFrame; kwargs)

Plot `gwas` results within a given `chromosome` and genomic range between `range1`
and `range2`.

Alternatively, plot within a given `chromosome` and a certain `window` around a
genomic coordinate `bp` or plot within a certain `window` around `gene`.

# Arguments
- `ld::Union{Nothing, SnpData, Tuple{SnpData, Union{AbstractString, Tuple{AbstractString, Int}}}} = nothing`:
    the reference panel for which LD is calculated.
- `ymax::Real`: the maximum value for y axis.
- `window::Real = 1e6`: the window around `bp` or `gene`.
- `showsnplabel::Bool = true`: draw SNP ID text label.
"""
function plotlocus!(
    ax::Axis,
    chromosome::AbstractString,
    range1::Real,
    range2::Real,
    gwas::DataFrame;
    ld::Union{Nothing, SnpData, Tuple{SnpData, Union{AbstractString, Tuple{AbstractString, Int}}}} = nothing,
    ymax::Real = 0,
    showsnplabel::Bool = true
)

    df = filter(x -> (x.CHR == chromosome) && (x.BP >= range1) && (x.BP <= range2), gwas)
    if nrow(df) == 0
        ymax == 0 ? ymax = 10 : nothing
        ax.spinewidth = 0.75
        ax.ytickwidth = 0.75
        ax.ylabelsize = 6
        ax.yticklabelsize = 6
        ax.yticksize = 3
        ax.yticks = setticks(ymax)
        xlims!(ax, range1, range2)
        ylims!(ax, 0, ymax)
        hidespines!(ax, :t, :r)
        hidexdecorations!(ax)
        hideydecorations!(ax, ticks = false, label = false, ticklabels = false)
        return
    end
    df.P = -log.(10, df.P)
    if ymax == 0
        ymax = maximum(df.P) / 4 * 5
        ymax <= 10 ? ymax = 10 : nothing
    end
    if !isnothing(ld)
        typeof(ld) == SnpData ? calculateld!(df, ld) : calculateld!(df, ld[1]; snp = ld[2])
        scatter!(ax, df.BP, df.P, color = df.LD, colorrange = (0, 1),
            colormap = (:gray60, :red2), markersize = 1.5)
        if typeof(ld) == SnpData
            ind = argmax(df.P)
            bp = df.BP[ind]
            p = df.P[ind]
            scatter!(ax, [bp], [p], color = :purple1, markersize = 4.0, marker = '◆')
            if showsnplabel
                text!(ax, "$(df.index[1])", position = (bp, p), textsize = 6, align = (:center, :bottom))
            end
        elseif length(df.index[1]) > 0
            ind = findfirst(df.SNP .== df.index[1])
            bp, p = df.BP[ind], df.P[ind]
            scatter!(ax, [bp], [p], color = :purple1, markersize = 4.0, marker = '◆')
            if showsnplabel
                text!(ax, "$(df.index[1])", position = (bp, p), textsize = 6, align = (:center, :bottom))
            end
        end
    else
        scatter!(ax, df.BP, df.P, color = :gray60, markersize = 1.5)
    end
    ax.spinewidth = 0.75
    ax.ytickwidth = 0.75
    ax.ylabelsize = 6
    ax.yticklabelsize = 6
    ax.yticksize = 3
    ax.yticks = setticks(ymax)
    xlims!(ax, range1, range2)
    ylims!(ax, 0, ymax)
    hidespines!(ax, :t, :r)
    hidexdecorations!(ax)
    hideydecorations!(ax, ticks = false, label = false, ticklabels = false)
end

plotlocus!(ax::Axis, chromosome::AbstractString, bp::Real, gwas::DataFrame; window::Real = 1e6, kwargs...) =
    plotlocus!(ax, chromosome, bp - window, bp + window, gwas; kwargs...)

function plotlocus!(ax::Axis, gene::AbstractString, gwas::DataFrame, gencode::DataFrame; window::Real = 1e6, kwargs...)
    chr, start, stop = findgene(gene, gencode)
    plotlocus!(ax, chr, start - window, stop + window, gwas; kwargs...)
end


text:
"""
**SeisWrite**

*Write seismic data in .seis format*

**IN**

* filename
* d: data
* h: headers as 1d array with elements of type Header
* extent: extent of the data (try _names(Extent)_ to see the information this contains)
* itrace=1

**OUT**

*Credits: AS, 2015*

"""
function SeisWrite(filename,d,h::Array{Header,1},extent::Extent;itrace=1)

    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_d = join([DATAPATH filename "@data@"])
    filename_h = join([DATAPATH filename "@headers@"])
    if (itrace==1)
	WriteTextHeader(filename,extent,"native_float",4,filename_d,filename_h)
	stream_dout = open(filename_d,"w")
	stream_hout = open(filename_h,"w")
    else
	stream_dout = open(filename_d,"a")
	stream_hout = open(filename_h,"a")
    end
    write(stream_dout,convert(Array{Float32,1},d[:]))
    close(stream_dout)
    nx = size(d[:,:],2)
    h1 = Header32Bits[]
    for j = itrace : itrace + nx - 1
	h[j - itrace + 1].tracenum = j
	h2 = HeaderToBits(h[j - itrace + 1])
	append!(h1,h2)
    end
    write(stream_hout,h1)
    close(stream_hout)
end

function SeisWrite(filename,d,extent::Extent)

    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_d = join([DATAPATH filename "@data@"])
    WriteTextHeader(filename,extent,"native_float",4,filename_d,"NULL")
    stream_dout = open(filename_d,"w")
    write(stream_dout,convert(Array{Float32,1},d[:]))
    close(stream_dout)
end


text:
# ------------------------------------------------------------------------------
#                        VML ARITHMETIC FUNCTIONS
# ------------------------------------------------------------------------------

if debug println("MKL : VML ARITHMETIC") end
for fname in (  :mkl_add!,
                :mkl_sub!,
                :mkl_sqr!,
                :mkl_mul!,
                :mkl_abs!,
                :mkl_linearfrac!)

                if debug println("export : ", fname) end
                @eval export $fname
end

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

## ?Add
fname = :mkl_add!
for (fmkl, elty) in ((:vdAdd,:Float64),
                      (:vsAdd,:Float32))

    @eval begin
        function $fname(n::Integer, a::Union{Ptr{$elty},DenseArray{$elty}}, b::Union{Ptr{$elty},DenseArray{$elty}}, y::Union{Ptr{$elty},DenseArray{$elty}})
            ccall(($(string(fmkl)), libmkl), Void,
                (Int, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}),
                 n, a, b, y)
            # vml_check_error()
            y
        end
    end
end
# ------------------------------------------------------------------------------

## ?Sub
fname = :mkl_sub!
for (fmkl, elty) in ((:vdSub,:Float64),
                      (:vsSub,:Float32))

    @eval begin
        function $fname(n::Integer, a::Union{Ptr{$elty},DenseArray{$elty}}, b::Union{Ptr{$elty},DenseArray{$elty}}, y::Union{Ptr{$elty},DenseArray{$elty}})
            ccall(($(string(fmkl)), libmkl), Void,
                (Int, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}),
                 n, a, b, y)
            # vml_check_error()
            y
        end
    end
end
# ------------------------------------------------------------------------------

## ?Sqr
fname = :mkl_sqr!
for (fmkl, elty) in ((:vdSqr,:Float64),
                      (:vsSqr,:Float32))

    @eval begin
        function $fname(n::Integer, a::Union{Ptr{$elty},DenseArray{$elty}}, y::Union{Ptr{$elty},DenseArray{$elty}})
            ccall(($(string(fmkl)), libmkl), Void,
                (Int, Ptr{$elty}, Ptr{$elty}),
                 n, a, y)
            # vml_check_error()
            y
        end
    end
end

## ?Mul
fname = :mkl_mul!
for (fmkl, elty) in ((:vdMul,:Float64),
                      (:vsMul,:Float32))

  @eval begin
    function $fname(n::Integer, a::Union{Ptr{$elty},DenseArray{$elty}}, b::Union{Ptr{$elty},DenseArray{$elty}}, y::Union{Ptr{$elty},DenseArray{$elty}})
        ccall(($(string(fmkl)), libmkl), Void,
            (Int, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}),
             n, a, b, y)
        # vml_check_error()
        y
    end
  end
end
# ------------------------------------------------------------------------------

## ?Abs
fname = :mkl_abs!
for (fmkl, elty) in ((:vdAbs,:Float64),
                      (:vsAbs,:Float32))

    @eval begin
      function $fname(n::Integer, a::Union{Ptr{$elty},DenseArray{$elty}}, y::Union{Ptr{$elty},DenseArray{$elty}})
          ccall(($(string(fmkl)), libmkl), Void,
              (Int, Ptr{$elty}, Ptr{$elty}),
               n, a, y)
          # vml_check_error()
          y
      end
    end
end
# ------------------------------------------------------------------------------

# ?LinearFrac
fname = :mkl_linearfrac
for (fmkl, elty) in ((:vdLinearFrac,:Float64),
                      (:vsLinearFrac,:Float32))

    @eval begin
        function $fname(n::Integer, a::Union{Ptr{$elty},DenseArray{$elty}}, b::Union{Ptr{$elty},DenseArray{$elty}},
          scalea::$elty, shifta::$elty, scaleb::$elty, shiftb::$elty, y::Union{Ptr{$elty},DenseArray{$elty}})
          ccall(($(string(fmkl)), libmkl), Void,
              (Int, Ptr{$elty}, Ptr{$elty}, $elty, $elty, $elty, $elty, Ptr{$elty}),
               n, a, b, scalea, shifta, scaleb, shiftb, y)
          # vml_check_error()
          y
        end
    end
end
# ------------------------------------------------------------------------------


text:
# test use of arbitrary arrays with Wood example from MINPACK

@testset "abstractarray" begin

c3 = 2e2
c4 = 2.02e1
c5 = 1.98e1
c6 = 1.8e2

function f!(F, x)
    temp1 = x[2] - x[1]^2
    temp2 = x[4] - x[3]^2
    F[1] = -c3*x[1]*temp1 - (1 - x[1])
    F[2] = c3*temp1 + c4*(x[2] - 1) + c5*(x[4] - 1)
    F[3] = -c6*x[3]*temp2 - (1 - x[3])
    F[4] = c6*temp2 + c4*(x[4] - 1) + c5*(x[2] - 1)
end

function j!(J, x)
    fill!(J, 0)
    temp1 = x[2] - 3x[1]^2
    temp2 = x[4] - 3x[3]^2
    J[1,1] = -c3*temp1 + 1
    J[1,2] = -c3*x[1]
    J[2,1] = -2*c3*x[1]
    J[2,2] = c3 + c4
    J[2,4] = c5
    J[3,3] = -c6*temp2 + 1
    J[3,4] = -c6*x[3]
    J[4,2] = c5
    J[4,3] = -2*c6*x[3]
    J[4,4] = c6 + c4
end

initial_x_matrix = [-3. -3; -1 -1]
initial_x = vec(initial_x_matrix)
#initial_x_wrapped = WrappedArray(initial_x_matrix)

for method in (:trust_region, :newton, :anderson, :broyden)
    r = nlsolve(f!, j!, initial_x, method = method)
    r_matrix = nlsolve(f!, j!, initial_x_matrix, method = method)
    #r_wrapped = nlsolve(f!, j!, initial_x_wrapped, method = method)

    @test r.zero == vec(r_matrix.zero)
    #@test r_matrix.zero == r_wrapped.zero
    @test r.residual_norm == r_matrix.residual_norm
    #@test r.residual_norm == r_wrapped.residual_norm
    @test typeof(r.zero) == typeof(initial_x)
    @test typeof(r_matrix.zero) == typeof(initial_x_matrix)
    #@test typeof(r_wrapped.zero) == typeof(initial_x_wrapped)
    r_AD = nlsolve(f!, initial_x, method = method, autodiff = :forward)
    r_matrix_AD = nlsolve(f!, initial_x_matrix, method = method, autodiff = :forward)
    #r_wrapped_AD = nlsolve(f!, initial_x_wrapped, method = method, autodiff = :forward)

    @test r_AD.zero == vec(r_matrix_AD.zero)
    #@test r_matrix_AD.zero == r_wrapped_AD.zero
    @test r_AD.residual_norm == r_matrix_AD.residual_norm
    #@test r_AD.residual_norm == r_wrapped_AD.residual_norm
    @test typeof(r_AD.zero) == typeof(initial_x)
    @test typeof(r_matrix_AD.zero) == typeof(initial_x_matrix)
    #@test typeof(r_wrapped_AD.zero) == typeof(initial_x_wrapped)
end
end
#=
when BandedMatrices are not super slow to precompile, we can add this test
remember to add REQUIRE in test and using BandedMatrices.
@testset "banded matrices" begin
function f!(fstor, x)
    fstor.=(1/3).*x.^3
end
function j!(jstor, x)
   for i = 1:3
       jstor[i,i] = x[i]^2
   end
end

dv = OnceDifferentiable(f!, j!, ones(3); J = bones(Float64, 3, 0, 0))

zer = nlsolve(dv, ones(3))
F = zeros(3)
f!(F, zer.zero)
@test norm(F, Inf) < 1e-8
@test typeof(NLsolve.jacobian(dv)) <: BandedMatrix
end
=#


text:
function draw_burn_arrow(brn::Burn; normalize=nothing, scale=nothing, color=(255,0,0), name="Burn")
    if isnothing(normalize)
        normalize = norm(orbital_velocity(0, Orbit(brn.orbit.primary.eqradius+100000, brn.orbit.primary)))
    end
    if isnothing(scale)
        scale = brn.orbit.primary.eqradius
    end

    length = norm(brn.Δv̄)/normalize*scale
    arrowdir = length*brn.Δv̄/norm(brn.Δv̄)
    tail_pos = time_orbital_position(brn.time, brn.orbit)
    head_pos = tail_pos .+ arrowdir
    components = orientation_components(brn)

    tail_trace = scatter3d(;x=[tail_pos[1], head_pos[1]],
                            y=[tail_pos[2], head_pos[2]],
                            z=[tail_pos[3], head_pos[3]],
                            mode="lines",
                            line=attr(color="rgb$color", ),
                            showlegend=false, hoverinfo="skip",
    )
    head_trace = cone(;x=[head_pos[1]],
                       y=[head_pos[2]],
                       z=[head_pos[3]],
                       u=[arrowdir[1]],
                       v=[arrowdir[2]],
                       w=[arrowdir[3]],
                       colorscale = [[0, "rgb$color"],
                                     [1, "rgb$color"]],
                       showscale = false,
                       showlegend = false,
                       name = name,
                       hovertemplate = join(["$(round(components[1],digits=2)) m/s prograde
",
                                             "$(round(components[3],digits=2)) m/s normal
",
                                             "$(round(components[2],digits=2)) m/s radial

",
                                             "UT: $(round(brn.time,digits=2)) s"]),
                       hoverlabel = orbit_hoverlabel(color),

    )
    return [tail_trace, head_trace]
end

function draw_ref_line(r; color=(255,255,255), dash="longdash", width=1)
    trace = scatter3d(;x=[0.,r], y=[0.,0.], z=[0.,0.,],
                       mode="lines",
                       line=attr(color="rgb$color",
                                 dash=dash,
                                 width=width,
                       ),
                       showlegend=false,
                       hoverinfo="skip",
    )
    return trace
end

function draw_angle_in_plane(vec::AbstractArray{<:Real}, vecref::AbstractArray{<:Real}, basis; npts=100, color=(255,255,255), linedash="dot", arcdash="dash", width=1, rscale=1.25, wrap_min=-π)
    r = rscale*max(norm(vec), norm(vecref))
    θ = angle_in_plane(vecref, basis)
    Δθ = wrap_angle(angle_in_plane(vec, vecref, basis), wrap_min)
    vecplane, vecrefplane = basis\vec, basis\vecref
    arcplane = fill(NaN, 3, npts)
    angles = collect(range(θ, stop=θ+Δθ, length=npts))
    for i=1:npts
        arcplane[:,i] = r*[cos(angles[i]), sin(angles[i]), 0.]
    end
    vecproj, vecrefproj = basis*@SVector([vecplane[1], vecplane[2], 0.]), basis*@SVector([vecrefplane[1], vecrefplane[2], 0.])
    arc = basis*arcplane

    traces = AbstractTrace[]
    # dashed lines
    xs = arc[1,:]
    prepend!(xs, [vecrefproj[1], arc[1,1]])
    append!(xs, vecproj[1], arc[1,end])
    ys = arc[2,:]
    prepend!(ys, [vecrefproj[2], arc[2,1]])
    append!(ys, vecproj[2], arc[2,end])
    zs = arc[3,:]
    prepend!(zs, [vecrefproj[3], arc[3,1]])
    append!(zs, vecproj[3], arc[3,end])
    push!(traces, scatter3d(;x=xs,
                             y=ys,
                             z=zs,
                             mode="lines",
                             line=attr(color="rgb$color",
                                       dash=linedash,
                                       width=width,
                             ),
                             hoverinfo="skip",
                             showlegend=false,
        )
    )
    # text
    sqrtrscale = sqrt(rscale)
    xmid, ymid, zmid = arc[:,Int(round(npts/2))]
    push!(traces, scatter3d(;x=[sqrtrscale*xmid],
                             y=[sqrtrscale*ymid],
                             z=[sqrtrscale*zmid],
                             mode="text",
                             text="$(round(rad2deg(Δθ), digits=2))°",
                             textfont = attr(color="rgb$color",
                                             size=12,
                                             family="Courier New, monospace",
                             ),
                             textposition="middle center",
                             hoverinfo="skip",
                             showlegend=false,
    )
)
    return traces
end
draw_angle_in_plane(vec::AbstractArray{<:Real}, orbref::Orbit, t) = draw_angle_in_plane(vec, time_orbital_position(t,orbref), orbref.basis)
draw_angle_in_plane(orb::Orbit, orbref::Orbit, t) = draw_angle_in_plane(time_orbital_position(t, orb), time_orbital_position(t,orbref), orbref.basis)

function draw_orbit_position(orb::Orbit, time; color=(255,255,255), name="", symbol="diamond")
    pos = time_orbital_position(time, orb)
    trace = scatter3d(;x=[pos[1]], y=[pos[2]], z=[pos[3]],
                       mode="markers",
                       marker=attr(color="rgb$color",
                                   symbol=symbol,
                                   size = 3,
                       ),
                       showlegend=false, name=name,
                       hovertemplate=name,
                       hoverlabel = orbit_hoverlabel(color),
    )
    return trace
end


text:
module FitnessPredictionsScript

using SynapseClient
using DISSEQT
using DISSEQT.SynapseTools
using DISSEQT.Plots
using DISSEQT.AnnotatedArrays
using LinearAlgebra
using DataFrames
using CSV
using Statistics
using JLD
using Gadfly
using Colors

# put mutagen doses into low,medium,high categories
function addlevel!(metadata)
    levelDict = Dict(("5FU","50uM")=>"Low", ("5FU","100uM")=>"Medium", ("5FU","150uM")=>"Medium", ("5FU","200uM")=>"High", ("AML","50uM")=>"Low", ("AML","100uM")=>"Medium", ("AML","200uM")=>"High", ("AZC","50uM")=>"Low", ("AZC","100uM")=>"Medium", ("AZC","200uM")=>"Medium", ("AZC","300uM")=>"High", ("Mn","0.25mM")=>"Medium", ("Mn","0.33mM")=>"Medium", ("Mn","0.5mM")=>"High", ("Mock","0uM")=>"Low", ("RIBA","50uM")=>"Low", ("RIBA","100uM")=>"Medium", ("RIBA","200uM")=>"Medium", ("RIBA","300uM")=>"High")
    metadata[:Level] = map( (m,d)->get(levelDict, (m,d), missing), metadata[:Mutagen], metadata[:Dose] )
end

function main()
    syn = SynapseClient.login()

    # projectFolder should point to MyProject folder (containing subfolders Metadata and Analysis)
    projectFolder  = "syn11639899" # VignuzziLabPublic/Projects/FitnessLandscapes
    analysisFolder = getchildbyname(syn, projectFolder, "Analysis")
    fitnessFolder  = getchildbyname(syn, projectFolder, "Fitness")

    isdir("MutantSwarms") || mkdir("MutantSwarms")
    isdir("plots") || mkdir("plots")

    # set to true to upload files
    doUpload = true
    uploadName = "fitness landscapes invitro"

    metadataID = getchildbyname(syn, projectFolder, "Metadata", "fitness landscapes invitro.csv")
    metadata = getsamplemetadata(syn, metadataID, :Rejected=>"No") # Always remove rejected samples
    metadata = metadata[ .~ismissing.(metadata[:Dose]), : ] # skip samples with unknown dose

    seqSpaceRepID = getchildbyname(syn, analysisFolder, "SequenceSpace", uploadName, "SequenceSpaceRepresentation", "seqspacerepresentation.jld")
    seqSpaceRepDict = load(localpath(syn, seqSpaceRepID))
    U = seqSpaceRepDict["U"]
    Σ = seqSpaceRepDict["Σ"]
    V = seqSpaceRepDict["V"]
    sampleID  = seqSpaceRepDict["SampleID"]
    positions = seqSpaceRepDict["positions"]
    codons    = seqSpaceRepDict["codons"]
    consensus = permutedims(seqSpaceRepDict["consensus"],(2,1)) # transpose to get Samples x Codons

    nbrDims = length(Σ)

    Y = V.*Σ' # sample representation scaled by singular values

    # also load the Consensus-based PCA representation
    seqSpaceRepConsensusID = getchildbyname(syn, analysisFolder, "SequenceSpace", uploadName, "SequenceSpaceRepresentation", "seqspacerepresentationConsensusPCA.jld")
    seqSpaceRepConsensusDict = load(localpath(syn, seqSpaceRepConsensusID))
    UConsensus = seqSpaceRepConsensusDict["U"]
    ΣConsensus = seqSpaceRepConsensusDict["Σ"]
    VConsensus = seqSpaceRepConsensusDict["V"]
    @assert sampleID == seqSpaceRepConsensusDict["SampleID"] # make sure the same samples were used
    nbrDimsConsensus = length(ΣConsensus)

    YConsensus = VConsensus.*ΣConsensus' # sample representation scaled by singular values

    # filter metadata by samples included in the Sequence Space Representation
    ind = indexin(sampleID, metadata[:SampleID])
    @assert all(ind.!=0) "Unknown SampleID"
    metadata = metadata[ind,:]

    addlevel!(metadata) # categorize dosages as low, medium, high

    # get fitness
    fitnessID = getchildbyname(syn,fitnessFolder,"fitness landscapes invitro fitness.csv")
    appendfitness!(metadata, getsamplefitness(syn, fitnessID))

    # setup dependencies
    dependencies = unique(vcat(seqSpaceRepID, seqSpaceRepConsensusID, metadata[:MetadataID], collect(skipmissing(metadata[:FitnessID]))))

    uploads = UploadList()

    virusColors = Dict("WT"  =>colorant"black",
                       "Less"=>colorant"blue",
                       "More"=>colorant"green",
                       "Stop"=>colorant"red")
    virusColorMap = groupcolors(metadata[:Reference],virusColors)

    mutagenColors = Dict("5FU" =>RGB(1,0,0),
                         "AML" =>RGB(0,1,0),
                         "AZC" =>RGB(0,0,1),
                         "Mn"  =>RGB(0.8,0.8,0),
                         "Mock"=>RGB(1,0,1),
                         "RIBA"=>RGB(0,1,1))
    mutagenColorMap = groupcolors(metadata[:Mutagen],mutagenColors)

    # Isomap representation
    K = Y*Y'
    D = sqrt.(diag(K) .+ diag(K)' .- 2*K) # sample distance matrix
    Z,stress,converged,dists = kruskalisomap(D, 2, 2, maximum(D)/20)
    converged || error("Kruskal's stress MDS did not converge.")

    # Isomap representation (Consensus PCA)
    KConsensus = YConsensus*YConsensus'
    DConsensus = sqrt.(max.(0., diag(KConsensus) .+ diag(KConsensus)' .- 2*KConsensus)) # sample distance matrix (make sure it's nonnegative to handle rounding errors for colocated points)
    ZConsensus,_,converged,_ = kruskalisomap(DConsensus, 2, 2, maximum(DConsensus)/20)
    converged || error("Kruskal's stress MDS did not converge for Consensus PCA.")

    # prediction
    fitnessMask = .~ismissing.(metadata[:Fitness])

    # make a table with all predictors
    predictionResults = DataFrame(Name=String[], Type=String[], VarianceExplained=Float64[], NbrSamples=Int[])

    # run twice, the second time excluding samples that constitute singleton groups in the "Lineage, Mutagen, Dose" group predictor
    mask = copy(fitnessMask)
    for i=1:2
        fitness = collect(skipmissing(metadata[mask,:Fitness]))

        # Isomap landscape
        σValues = 10.0.^range(-3,stop=1,length=1000) * sqrt(mean(sum(Z[mask,:].^2; dims=2))) # scale by point size cloud
        @time σPerSample,_,_ = leaveoneoutlandscapekernelwidth(Z[mask,:], fitness, σValues, nbrIter=1000)
        f = leaveoneoutpredict(LandscapeModel, Z[mask,:], fitness, perSampleModelArgs=σPerSample)
        push!(predictionResults, ("[GKS] Isomap (2d)", "Gaussian Kernel Smoother", 1.0.-varianceunexplained(f,fitness), count(mask)) )
        # println(predictionResults)

        # Reduced space landscape
        σValues = 10.0.^range(-3,stop=1,length=1000) * sqrt(mean(sum(Y[mask,:].^2; dims=2))) # scale by point size cloud
        @time σPerSample,_,_ = leaveoneoutlandscapekernelwidth(Y[mask,:], fitness, σValues, nbrIter=1000)
        f = leaveoneoutpredict(LandscapeModel, Y[mask,:], fitness, perSampleModelArgs=σPerSample)
        push!(predictionResults, ("[GKS] SMSSVD ($(nbrDims)d)", "Gaussian Kernel Smoother", 1.0.-varianceunexplained(f,fitness), count(mask)) )
        # println(predictionResults)

        # Isomap landscape (Consensus PCA)
        σValues = 10.0.^range(-3,stop=1,length=1000) * sqrt(mean(sum(ZConsensus[mask,:].^2; dims=2))) # scale by point size cloud
        @time σPerSample,_,_ = leaveoneoutlandscapekernelwidth(ZConsensus[mask,:], fitness, σValues, nbrIter=1000)
        f = leaveoneoutpredict(LandscapeModel, ZConsensus[mask,:], fitness, perSampleModelArgs=σPerSample)
        push!(predictionResults, ("[GKS] Consensus Isomap (2d)", "Gaussian Kernel Smoother", 1.0.-varianceunexplained(f,fitness), count(mask)) )
        # println(predictionResults)

        # Reduced space landscape  (Consensus PCA)
        σValues = 10.0.^range(-3,stop=1,length=1000) * sqrt(mean(sum(YConsensus[mask,:].^2; dims=2))) # scale by point size cloud
        @time σPerSample,_,_ = leaveoneoutlandscapekernelwidth(YConsensus[mask,:], fitness, σValues, nbrIter=1000)
        f = leaveoneoutpredict(LandscapeModel, YConsensus[mask,:], fitness, perSampleModelArgs=σPerSample)
        push!(predictionResults, ("[GKS] Consensus PCA ($(nbrDimsConsensus)d)", "Gaussian Kernel Smoother", 1.0.-varianceunexplained(f,fitness), count(mask)) )

        # Nearest neighbor predictor: Reduced dimension
        f = leaveoneoutpredict(NearestNeighborModel, Y[mask,:], fitness)
        push!(predictionResults, ("[NN] SMSSVD ($(nbrDims)d)", "Nearest Neighbor", 1.0.-varianceunexplained(f,fitness), count(mask)) )

        # Nearest neighbor predictor: Consensus
        f = leaveoneoutpredict(NearestNeighborModel, consensus[mask,:], fitness)
        push!(predictionResults, ("[NN] Consensus", "Nearest Neighbor", 1.0.-varianceunexplained(f,fitness), count(mask)) )

        # Group predictor: Lineage, Drug
        f = leaveoneoutpredict(GroupModel, metadata[mask,[:Reference,:Mutagen]], fitness)
        push!(predictionResults, ("[G] Lineage/Mutagen", "Group", 1.0.-varianceunexplained(f,fitness), count(mask)) )

        # Group predictor: Lineage, Level
        f = leaveoneoutpredict(GroupModel, metadata[mask,[:Reference,:Level]], fitness)
        push!(predictionResults, ("[G] Lineage/Dose", "Group", 1.0.-varianceunexplained(f,fitness), count(mask)) )

        # Group predictor: Lineage, Drug, Dose
        #f = leaveoneoutpredict(GroupModel, metadata[mask,[:Reference,:Mutagen,:Dose]], fitness)
        f = leaveoneoutpredict(GroupModel, metadata[mask,[:Reference,:Mutagen,:Level]], fitness)
        push!(predictionResults, ("[G] Lineage/Mutagen/Dose", "Group", 1.0.-varianceunexplained(f,fitness), count(mask)) )
        newMask = i==1 ? .~isnan.(f) : trues(size(f)) # For next iteration, exclude singleton groups
        i==1 && println(count(.~newMask), " singleton groups removed.")

        mask[mask] = newMask
    end

    println(predictionResults)
    CSV.write("plots/fitnesspredictions.csv", predictionResults; missingstring="NA")
    markforupload!(uploads, "FitnessPredictionPlots", "plots/fitnesspredictions.csv", dependencies)

    N = div(size(predictionResults,1),2)
    xMax = maximum(predictionResults[.~isnan.(predictionResults[:VarianceExplained]), :VarianceExplained])

    # with all samples
    df = predictionResults[1:N,:]
    df = df[.~isnan.(df[:VarianceExplained]),:] # remove predictors that couldn't be used ([G] Lineage/Mutagen/Dose)
    plotArgs = (layer(df, x=:VarianceExplained, color=:Name, Geom.bar(orientation=:horizontal)),
                Coord.Cartesian(xmin=0, xmax=xMax, yflip=true),
                Guide.yticks(ticks=nothing),
                Guide.xlabel("Variance Explained"), Guide.ylabel(""))
    themeArgs = (:key_title_font_size=>2*11pt, :key_label_font_size=>2*8pt, :major_label_font_size=>2*11pt, :minor_label_font_size=>2*8pt, :point_label_font_size=>2*8pt)

    pl = plot(plotArgs..., Theme(;themeArgs...))
    markforupload!(uploads, "FitnessPredictionPlots", saveplot(pl, [:png,:svg,:pdf], "plots/fitnesspredictions_allsamples"), dependencies)

    pl = plot(plotArgs..., Theme(key_position=:none; themeArgs...)) # no legend
    markforupload!(uploads, "FitnessPredictionPlots", saveplot(pl, [:png,:svg,:pdf], "plots/fitnesspredictions_allsamples_nolegend"), dependencies)

    # with singleton groups removed
    df = predictionResults[N+1:end,:]
    df = df[.~isnan.(df[:VarianceExplained]),:] # remove predictors that couldn't be used (none)
    plotArgs = (layer(df, x=:VarianceExplained, color=:Name, Geom.bar(orientation=:horizontal)),
                Coord.Cartesian(xmin=0, xmax=xMax, yflip=true),
                Guide.yticks(ticks=nothing),
                Guide.xlabel("Variance Explained"), Guide.ylabel(""))
    themeArgs = (:key_title_font_size=>2*11pt, :key_label_font_size=>2*8pt, :major_label_font_size=>2*11pt, :minor_label_font_size=>2*8pt, :point_label_font_size=>2*8pt)

    pl = plot(plotArgs..., Theme(;themeArgs...))
    markforupload!(uploads, "FitnessPredictionPlots", saveplot(pl, [:png,:svg,:pdf], "plots/fitnesspredictions"), dependencies)

    pl = plot(plotArgs..., Theme(key_position=:none; themeArgs...)) # no legend
    markforupload!(uploads, "FitnessPredictionPlots", saveplot(pl, [:png,:svg,:pdf], "plots/fitnesspredictions_nolegend"), dependencies)

    doUpload && uploadfiles(syn, createchildfolder(syn, analysisFolder, "SequenceSpace", uploadName), uploads, @__FILE__, "Sequence Space Plots")

    # seqSpaceRepDict, metadata, U, Σ, V, Y, Z
    nothing
end

end
FitnessPredictionsScript.main();
# seqSpaceRepDict, metadata, U, Σ, V, Y, Z = FitnessPredictionsScript.main();


text:
"""model definition of arm33"""

struct StateArm33 <: AbsState
    θ1::Float64
    ϕ1::Float64
    θ2::Float64
    ϕ2::Float64
    θ3::Float64
    ϕ3::Float64
end

function get_mid_status(p::StateArm33, q::StateArm33)::StateArm33
    return StateArm33(
        diff_angles(q.θ1, p.θ1) / 2 + p.θ1,
        diff_angles(q.ϕ1, p.ϕ1) / 2 + p.ϕ1,
        diff_angles(q.θ2, p.θ2) / 2 + p.θ2,
        diff_angles(q.ϕ2, p.ϕ2) / 2 + p.ϕ2,
        diff_angles(q.θ3, p.θ3) / 2 + p.θ3,
        diff_angles(q.ϕ3, p.ϕ3) / 2 + p.ϕ3,
    )
end

function dist(q1::StateArm33, q2::StateArm33)
    return norm([
        diff_angles(q1.θ1, q2.θ1) / π,
        diff_angles(q1.ϕ1, q2.ϕ1) / π,
        diff_angles(q1.θ2, q2.θ2) / π,
        diff_angles(q1.ϕ2, q2.ϕ2) / π,
        diff_angles(q1.θ3, q2.θ3) / π,
        diff_angles(q1.ϕ3, q2.ϕ3) / π,
    ])
end

function get_arm33_positions(q::StateArm33, pos::Vector{Float64}, rad::Float64)
    a = pos
    b = a + [sin(q.θ1) * cos(q.ϕ1), sin(q.θ1) * sin(q.ϕ1), cos(q.θ1)] * rad
    c = b + [sin(q.θ2) * cos(q.ϕ2), sin(q.θ2) * sin(q.ϕ2), cos(q.θ2)] * rad
    d = c + [sin(q.θ3) * cos(q.ϕ3), sin(q.θ3) * sin(q.ϕ3), cos(q.θ3)] * rad
    return [a, b, c, d]
end

function gen_connect(
    q::StateArm33,
    obstacles::Vector{CircleObstacle3D},
    positions::Vector{Vector{Float64}},
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
    max_dist::Union{Nothing,Float64} = nothing,
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Function

    # check: q \in C_free
    f(q::StateArm33, i::Int64)::Bool = begin
        P = get_arm33_positions(q, positions[i], rads[i])

        any(x -> (x < 0 || 1 < x), vcat(P...)) && return false

        # obstacle
        for k = 2:4
            any(o -> dist(P[k-1], P[k], [o.x, o.y, o.z]) < o.r, obstacles) && return false
        end

        # self collision
        dist(P[1], P[2], P[3], P[4]) < safety_dist && return false

        return true
    end

    f(q_from::StateArm33, q_to::StateArm33, i::Int64)::Bool = begin
        D = dist(q_from, q_to)
        !isnothing(max_dist) && D > max_dist && return false

        dt1_θ = diff_angles(q_to.θ1, q_from.θ1)
        dt1_ϕ = diff_angles(q_to.ϕ1, q_from.ϕ1)
        dt2_θ = diff_angles(q_to.θ2, q_from.θ2)
        dt2_ϕ = diff_angles(q_to.ϕ2, q_from.ϕ2)
        dt3_θ = diff_angles(q_to.θ3, q_from.θ3)
        dt3_ϕ = diff_angles(q_to.ϕ3, q_from.ϕ3)

        for e in vcat(collect(0:step_dist:D) / D, 1.0)
            # angles
            arr_θ = [q_from.θ1 + e * dt1_θ, q_from.θ2 + e * dt2_θ, q_from.θ3 + e * dt3_θ]
            arr_ϕ =
                [q_from.ϕ1 + e * dt1_ϕ, q_from.ϕ2 + e * dt2_ϕ, q_from.ϕ3 + e * dt3_ϕ]
            # positions
            P = [positions[i]]
            for k = 2:4
                push!(
                    P,
                    P[k-1] +
                    rads[i] * [
                        sin(arr_θ[k-1]) * cos(arr_ϕ[k-1]),
                        sin(arr_θ[k-1]) * sin(arr_ϕ[k-1]),
                        cos(arr_θ[k-1]),
                    ],
                )
            end

            # outside
            any(x -> (x < 0 || 1 < x), vcat(P...)) && return false

            # obstacles
            for k = 2:4
                any(o -> dist(P[k-1], P[k], [o.x, o.y, o.z]) < o.r, obstacles) &&
                    return false
            end

            # self collision
            dist(P[1], P[2], P[3], P[4]) < safety_dist && return false
        end

        return true
    end
    return f
end

function gen_collide(
    q::StateArm33,
    positions::Vector{Vector{Float64}},
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Function

    N = length(rads)

    f(
        q_i_from::StateArm33,
        q_i_to::StateArm33,
        q_j_from::StateArm33,
        q_j_to::StateArm33,
        i::Int64,
        j::Int64,
    ) = begin
        # check each pair of step
        D_i = dist(q_i_from, q_i_to)
        D_j = dist(q_j_from, q_j_to)

        dt1_θ_i = diff_angles(q_i_to.θ1, q_i_from.θ1)
        dt1_ϕ_i = diff_angles(q_i_to.ϕ1, q_i_from.ϕ1)
        dt2_θ_i = diff_angles(q_i_to.θ2, q_i_from.θ2)
        dt2_ϕ_i = diff_angles(q_i_to.ϕ2, q_i_from.ϕ2)
        dt3_θ_i = diff_angles(q_i_to.θ3, q_i_from.θ3)
        dt3_ϕ_i = diff_angles(q_i_to.ϕ3, q_i_from.ϕ3)

        dt1_θ_j = diff_angles(q_j_to.θ1, q_j_from.θ1)
        dt1_ϕ_j = diff_angles(q_j_to.ϕ1, q_j_from.ϕ1)
        dt2_θ_j = diff_angles(q_j_to.θ2, q_j_from.θ2)
        dt2_ϕ_j = diff_angles(q_j_to.ϕ2, q_j_from.ϕ2)
        dt3_θ_j = diff_angles(q_j_to.θ3, q_j_from.θ3)
        dt3_ϕ_j = diff_angles(q_j_to.ϕ3, q_j_from.ϕ3)

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # intermediate angles & positions for agent-i
            arr_θ_i = [
                q_i_from.θ1 + e_i * dt1_θ_i,
                q_i_from.θ2 + e_i * dt2_θ_i,
                q_i_from.θ3 + e_i * dt3_θ_i,
            ]
            arr_ϕ_i = [
                q_i_from.ϕ1 + e_i * dt1_ϕ_i,
                q_i_from.ϕ2 + e_i * dt2_ϕ_i,
                q_i_from.ϕ3 + e_i * dt3_ϕ_i,
            ]
            # positions
            P_i = [positions[i]]
            for k = 2:4
                push!(
                    P_i,
                    P_i[k-1] +
                    rads[i] * [
                        sin(arr_θ_i[k-1]) * cos(arr_ϕ_i[k-1]),
                        sin(arr_θ_i[k-1]) * sin(arr_ϕ_i[k-1]),
                        cos(arr_θ_i[k-1]),
                    ],
                )
            end

            for e_j in vcat(collect(0:step_dist:D_j) / D_j, 1.0)
                # intermediate angles & positions for agent-j
                arr_θ_j = [
                    q_j_from.θ1 + e_j * dt1_θ_j,
                    q_j_from.θ2 + e_j * dt2_θ_j,
                    q_j_from.θ3 + e_j * dt3_θ_j,
                ]
                arr_ϕ_j = [
                    q_j_from.ϕ1 + e_j * dt1_ϕ_j,
                    q_j_from.ϕ2 + e_j * dt2_ϕ_j,
                    q_j_from.ϕ3 + e_j * dt3_ϕ_j,
                ]
                # positions
                P_j = [positions[j]]
                for k = 2:4
                    push!(
                        P_j,
                        P_j[k-1] +
                        rads[j] * [
                            sin(arr_θ_j[k-1]) * cos(arr_ϕ_j[k-1]),
                            sin(arr_θ_j[k-1]) * sin(arr_ϕ_j[k-1]),
                            cos(arr_θ_j[k-1]),
                        ],
                    )
                end

                # check_collision
                for k = 2:4, l = 2:4
                    dist(P_i[k-1], P_i[k], P_j[l-1], P_j[l]) < safety_dist && return true
                end
            end
        end
        return false
    end

    f(q_i::StateArm33, q_j::StateArm33, i::Int64, j::Int64) = begin
        P_i = get_arm33_positions(q_i, positions[i], rads[i])
        P_j = get_arm33_positions(q_j, positions[j], rads[j])

        # check_collision
        for k = 2:4, l = 2:4
            dist(P_i[k-1], P_i[k], P_j[l-1], P_j[l]) < safety_dist && return true
        end
        return false
    end

    f(Q_from::Vector{Node{StateArm33}}, Q_to::Vector{Node{StateArm33}}) = begin
        for i = 1:N, j = i+1:N
            if f(Q_from[i].q, Q_to[i].q, Q_from[j].q, Q_to[j].q, i, j)
                return true
            end
        end
        return false
    end

    f(Q::Vector{Node{StateArm33}}, q_i_to::StateArm33, i::Int64) = begin
        q_i_from = Q[i].q
        D_i = dist(q_i_from, q_i_to)

        dt1_θ_i = diff_angles(q_i_to.θ1, q_i_from.θ1)
        dt1_ϕ_i = diff_angles(q_i_to.ϕ1, q_i_from.ϕ1)
        dt2_θ_i = diff_angles(q_i_to.θ2, q_i_from.θ2)
        dt2_ϕ_i = diff_angles(q_i_to.ϕ2, q_i_from.ϕ2)
        dt3_θ_i = diff_angles(q_i_to.θ3, q_i_from.θ3)
        dt3_ϕ_i = diff_angles(q_i_to.ϕ3, q_i_from.ϕ3)

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # intermediate angles & positions for agent-i
            arr_θ_i = [
                q_i_from.θ1 + e_i * dt1_θ_i,
                q_i_from.θ2 + e_i * dt2_θ_i,
                q_i_from.θ3 + e_i * dt3_θ_i,
            ]
            arr_ϕ_i = [
                q_i_from.ϕ1 + e_i * dt1_ϕ_i,
                q_i_from.ϕ2 + e_i * dt2_ϕ_i,
                q_i_from.ϕ3 + e_i * dt3_ϕ_i,
            ]
            # positions
            P_i = [positions[i]]
            for k = 2:4
                push!(
                    P_i,
                    P_i[k-1] +
                    rads[i] * [
                        sin(arr_θ_i[k-1]) * cos(arr_ϕ_i[k-1]),
                        sin(arr_θ_i[k-1]) * sin(arr_ϕ_i[k-1]),
                        cos(arr_θ_i[k-1]),
                    ],
                )
            end

            for j in filter(k -> k != i, 1:N)
                P_j = get_arm33_positions(Q[j].q, positions[j], rads[j])

                # check collision
                for k = 2:4, l = 2:4
                    dist(P_i[k-1], P_i[k], P_j[l-1], P_j[l]) < safety_dist && return true
                end
            end
        end

        return false
    end

    return f
end

function gen_uniform_sampling(q::StateArm33)::Function
    () -> begin
        return StateArm33((rand(6) * 2π)...)
    end
end

function plot_motion!(
    q_from::StateArm33,
    q_to::StateArm33,
    pos::Vector{Float64},
    rad::Float64,
    params,
)
    p1 = get_arm33_positions(q_from, pos, rad)[end]
    p2 = get_arm33_positions(q_to, pos, rad)[end]
    plot!([p1[1], p2[1]], [p1[2], p2[2]], [p1[3], p2[3]]; params...)
end

function plot_start_goal!(
    q_init::StateArm33,
    q_goal::StateArm33,
    pos::Vector{Float64},
    rad::Float64,
    params,
)
    pos_init = get_arm33_positions(q_init, pos, rad)[end]
    pos_goal = get_arm33_positions(q_goal, pos, rad)[end]
    plot!([pos_init[1]], [pos_init[2]], [pos_init[3]]; markershape = :hex, params...)
    plot!([pos_goal[1]], [pos_goal[2]], [pos_goal[3]]; markershape = :star, params...)
end

function plot_agent!(q::StateArm33, pos::Vector{Float64}, rad::Float64, color::String)
    P = get_arm33_positions(q, pos, rad)
    plot!(
        map(p -> p[1], P),
        map(p -> p[2], P),
        map(p -> p[3], P),
        color = color,
        lw = 5,
        markershape = :circle,
        markersize = 3,
        label = nothing,
    )
    scatter!(
        [pos[1]],
        [pos[2]],
        [pos[3]],
        markersize = 10,
        markershape = :rect,
        color = color,
        label = nothing,
    )
end

function gen_random_instance_StateArm33(;
    N_min::Int64 = 2,
    N_max::Int64 = 8,
    N::Int64 = rand(N_min:N_max),
    num_obs_min::Int64 = 0,
    num_obs_max::Int64 = 10,
    num_obs::Int64 = rand(num_obs_min:num_obs_max),
    rad::Float64 = 0.025,
    rad_obs::Float64 = 0.05,
    rad_min::Float64 = rad,
    rad_max::Float64 = rad,
    rad_obs_min::Float64 = rad_obs,
    rad_obs_max::Float64 = rad_obs,
    TIME_LIMIT::Float64 = 0.5,
)::Tuple{
    Vector{StateArm33},  # initial configuration
    Vector{StateArm33},  # goal configuration
    Vector{CircleObstacle3D},
    Vector{Vector{Float64}},  # positions
    Vector{Float64},  # radius
}

    _q = StateArm33(0, 0, 0, 0, 0, 0)
    while true
        t_s = now()
        timeover() = elapsed_sec(t_s) > TIME_LIMIT

        # generate obstacles
        obstacles = gen_obstacles(3, num_obs, rad_obs_min, rad_obs_max)

        # determine positions
        positions = map(e -> rand(3), 1:N)

        # determine rads
        rads = map(e -> rand() * (rad_max - rad_min) + rad_min, 1:N)

        # generate starts & goals
        connect = gen_connect(_q, obstacles, positions, rads)
        collide = gen_collide(_q, positions, rads)
        config_init, config_goal = gen_config_init_goal(_q, N, connect, collide, timeover)

        timeover() && continue
        return (config_init, config_goal, obstacles, positions, rads)
    end
end


text:
struct SequenceGraph
    nodes::Vector{SequenceGraphNode}
    links::Vector{Vector{SequenceGraphLink}}
end

nodes(sg::SequenceGraph) = sg.nodes
node(sg::SequenceGraph, i::NodeID) = nodes(sg)[abs(i)]
links(sg::SequenceGraph) = sg.links

"""
    links(sg::SequenceGraph, node::NodeID)

Get all of the links of a Node of a sequence graph.
"""
function links(sg::SequenceGraph, node::NodeID)
    l = links(sg)
    @inbounds return l[abs(node)]
end

"""
    add_node!(sg::SequenceGraph, n::SequenceGraphNode)

Add a SequenceGraphNode `n` to a SequenceGraph `sg`.

Returns the node ID used to access the node added from the graph.
"""
function add_node!(sg::SequenceGraph, n::SequenceGraphNode)
    newlen = length(push!(nodes(sg)))
    resize!(links(sg), newlen)
    return newlen
end

function remove_node!(sg::SequenceGraph, n::NodeID)
    oldlinks = copy(links(sg, n))
    for oldlink in oldlinks
        remove_link!(sg, source(oldlink), destination(oldlink))
    end
end

"""
    add_link!(sg::SequenceGraph, source::NodeID, dest::NodeID, dist::Int)

Construct a link between two nodes in a sequence Graph.
"""
function add_link!(sg::SequenceGraph, source::NodeID, dest::NodeID, dist::Int)
    push!(links(sg, source), Link(source, dest, dist))
    push!(links(sg, dest), Link(dest, source, dist))
end

function remove_link!(sg::SequenceGraph, source::NodeID, dest::NodeID)
    slinks = links(sg, source)
    slinkslen = length(slinks)
    filter!(!isequal(SequenceGraphLink(source, dest, 0)), slinks)
    dlinks = links(sg, dest)
    dlinkslen = length(dlinks)
    filter!(!isequal(SequenceGraphLink(dest, source, 0)), dlinks)
    return slinkslen != length(slinks) || dlinkslen != length(dlinks)
end

"""
    forward_links(sg::SequenceGraph, n::NodeID)

Get the list of links leaving `n`.
"""
function forward_links(sg::SequenceGraph, n::NodeID)
    r = Vector{Link}(0)
    nodelinks = links(sg, n)
    sizehint!(r, length(nodelinks))
    for link in nodelinks
        if is_forward_link(link, n)
            push!(r, link)
        end
    end
    return r
end

backward_links(sg::SequenceGraph, n::NodeID) = forward_links(sg, -n)


text:
@testset "Polygon" begin
    @testset "Triangle2D" begin
        for T in [Float32, Float64, BigFloat]
            p₁ = Point2D{T}(0, 0)
            p₂ = Point2D{T}(1, 0)
            p₃ = Point2D{T}(0, 1)
            tri = Triangle(p₁, p₂, p₃)
            @test tri.points == SVector(p₁, p₂, p₃)

            @test isconvex(tri)
        end
    end

    @testset "Triangle3D" begin
        for T in [Float32, Float64, BigFloat]
            p₁ = Point3D{T}(0, 0, 0)
            p₂ = Point3D{T}(0, 1, 0)
            p₃ = Point3D{T}(0, 0, 1)
            tri = Triangle(p₁, p₂, p₃)
            @test tri.points == SVector(p₁, p₂, p₃)

            @test isplanar(tri)
        end
    end

    @testset "Quadrilateral2D" begin
        for T in [Float32, Float64, BigFloat]
            p₁ = Point2D{T}(0, 0)
            p₂ = Point2D{T}(1, 0)
            p₃ = Point2D{T}(1, 1)
            p₄ = Point2D{T}(0, 1)
            quad = Quadrilateral(p₁, p₂, p₃, p₄)
            @test quad.points == SVector(p₁, p₂, p₃, p₄)

            @test isconvex(quad)
            @test !isconvex(Quadrilateral(p₄, p₃, p₂, p₁))
        end
    end

    @testset "Quadrilateral3D" begin
        for T in [Float32, Float64, BigFloat]
            p₁ = Point3D{T}(0, 0, 0)
            p₂ = Point3D{T}(0, 1, 0)
            p₃ = Point3D{T}(0, 1, 1)
            p₄ = Point3D{T}(0, 0, 1)
            quad = Quadrilateral(p₁, p₂, p₃, p₄)
            @test quad.points == SVector(p₁, p₂, p₃, p₄)

            @test isplanar(quad)
            @test !isplanar(Quadrilateral(p₁, p₂, p₃, Point3D{T}(1//10, 0, 1)))
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using NumericalMethodsforEngineers

old = pwd()
ProjDir = Pkg.dir("NumericalMethodsforEngineers", "examples", "NumericalMethodsforEngineers", "Ch07", "BVP")
cd(ProjDir)

function f(x::Float64, y::Vector{Float64})
  [y[2], 3x^2+4y[1]]
end

(iter, nsteps, converged, res) = shootingmethod(f, [0.0, 0.0, 1.0, 1.0], [0.0, 1.0], 5)
@show round.(res, 4)

cd(old)


text:
#####
##### Plotting Utility Functions
#####

rectangle(x, y, w, h) = (x .+ [0, w, w, 0]), (y .+ [0, 0, h, h])
getname(v, s::Symbol) = getindex.(v, s)

load_save_files(f, formulations::String) = first(load_save_files(f, (formulations,)))
function load_save_files(f, formulations)
    savefiles = [joinpath(savedir(f), join((name(f), i), "_") * ".jls") for i in formulations]
    data = deserialize.(savefiles)
    for d in data
        sort!(d.runs; rev = true, by = x -> get(x, :dram_limit, 0))
    end

    return data
end

# For drawing a vertical asymptote on a graph
vasymptote() = """
\\pgfplotsset{vasymptote/.style={
    before end axis/.append code={
        \\draw[densely dashed] ({rel axis cs:0,0} -| {axis cs:#1,0})
        -- ({rel axis cs:0,1} -| {axis cs:#1,0});
}
}}
"""

hasymptote() = """
\\pgfplotsset{hasymptote/.style={
    before end axis/.append code={
        \\draw[densely dashed] ({rel axis cs:0,0} -| {axis cs:0,#1})
        -- ({rel axis cs:1,0} -| {axis cs:0,#1});
}
}}
"""

hline(y; xl = 0, xu = 1, color = "red") = """
\\draw[$color, densely dashed, ultra thick]
    ({axis cs:$xl,$y} -| {rel axis cs:0,0}) --
    ({axis cs:$xu,$y} -| {rel axis cs:1,0});
""" |> rm_newlines

vline(x; yl = 0, yu = 1, color = "red") = """
\\draw[$color, sharp plot] ($x, $yl) -- ($x, $yu);
""" |> rm_newlines

comment(x, y, str; kw...) = @pgf(["\\node[align = center$(format(kw))] at ", Coordinate(x, y), "{$str};"])
format(kw) = ", " * join(["$a = $b" for (a,b) in kw], ", ")

rm_newlines(str) = join(split(str, "\n"))

# Node - must sort data before hand
# using load_save_files does this automatically
get_dram_performance(data) = minimum(get_dram_performance.(data))
get_dram_performance(data::NamedTuple) = first(getname(data.runs, :actual_runtime))

get_pmm_performance(data) = maximum(get_pmm_performance.(data))
get_pmm_performance(data::NamedTuple) = last(getname(data.runs, :actual_runtime))

function findabsmin(f, x)
    _, ind = findmin(abs.(f.(x)))
    return ind
end

#####
##### ILP Solution Times
#####

function pgf_ilp_solution_times(
        fns,
        cache;
        file = "plot.tex",
        formulations = ("synchronous",),
        backend = nGraph.Backend("CPU"),
        width = "10cm",
        height = "6cm",
    )

    # Collect the symbolic coordinates
    coordinates = String[]
    for f in fns
        push!(coordinates, "$(AutoTM.Experiments.titlename(f))")
    end

    @show coordinates

    # Construct bar plots
    plots = []
    for formulation in formulations
        x = []
        y = []

        for f in fns
            data = deserialize(canonical_path(f, formulation, cache, backend))

            # Extract the extrema of optimization.
            #
            # We map the `sum` function because of multiple solution attempts due to
            # fragmentation.
            min_opt, max_opt = extrema(sum, getname(data.runs, :optimization_times))

            push!(x, "$(AutoTM.Experiments.titlename(f))")
            push!(y, max_opt)
        end

        append!(plots, [
            @pgf(PlotInc(Coordinates(x, y))),
            @pgf(LegendEntry(get(_speedup_formulation_lut, formulation, formulation)))
        ])
    end

    document = TikzDocument()
    push!(document, "\\pgfplotsset{width=$width,height=$height}")

    picture = TikzPicture()
    axis = @pgf Axis(
        {
            # Style
            ybar,

            # Axis Setup
            ymin = 0,
            ymode = "log",
            bar_width = "18pt",
            xtick = "data",
            symbolic_x_coords = coordinates,

            # Labels
            ylabel = "ILP Solution Time (s)\\\\\\textbf{Logarithmic}",
            ylabel_style = {
                align = "center",
                rotate = -90,
            },

            # Legend Setup
            legend_style =
            {
                 at = Coordinate(0.05, 1.05),
                 anchor = "south west",
                 legend_columns = -1
            },
        },
        plots...,
    )

    push!(picture, axis)
    push!(document, picture)
    pgfsave(file, document)
    return document
end

#####
##### Statistics on what AutoTM has done.
#####

_wh() = ("8cm", "6cm")
_ms() = 1.5

function pgf_movement_plot(f, cache, suffix; file = "plot.tex", formulation = "synchronous")
    savefile = canonical_path(f, formulation, cache, nGraph.Backend("CPU"), suffix)
    data = deserialize(savefile)

    # Plot the number of move nodes.
    io_size = data.io_size[]
    dram_sizes = (getname(data.runs, :dram_limit) ./ 1E3) .+ (io_size ./ 1E9)

    x = dram_sizes

    plot = TikzDocument()
    push!(plot, vasymptote())

    plots = [
        @pgf(PlotInc(
             {
                thick,
                mark_options = {
                    scale = _ms(),
                },
             },
             Coordinates(x, getname(data.runs, :bytes_moved_pmem) ./ 1E9)
        )),
        @pgf(LegendEntry("sync DRAM to PMEM")),
        @pgf(PlotInc(
             {
                thick,
                mark_options = {
                    scale = _ms(),
                },
             },
             Coordinates(x, getname(data.runs, :bytes_moved_dram) ./ 1E9)
        )),
        @pgf(LegendEntry("sync PMEM to DRAM")),
    ]

    width, height = _wh()
    axs = @pgf Axis(
        {
            grid = "major",
            xlabel = "DRAM Limit (GB)",
            ylabel = "Memory Moved (GB)",
            width = width,
            height = height,
            vasymptote = data.default_alloc_size[] / 1E9,
        },
        plots...,
    )

    push!(plot, TikzPicture(axs))

    empty!(PGFPlotsX.CUSTOM_PREAMBLE)
    push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{colorbrewer}")

    pgfsave(file, plot)
    return plot
end

function pgf_io_plot(f, cache, suffix; file = "plot.tex", formulations = ("synchronous",))
    savefiles = canonical_path.(
        Ref(f),
        formulations,
        Ref(cache),
        Ref(nGraph.Backend("CPU")),
        Ref(suffix)
    )
    data = deserialize.(savefiles)

    # Plot the number of move nodes.
    io_size = first(data).io_size[]
    plots = []
    for (d, formulation) in zip(data, formulations)
        dram_sizes = (getname(d.runs, :dram_limit) ./ 1E3) .+ (io_size ./ 1E9)
        push!(plots, @pgf(
            PlotInc(
                 {
                    thick,
                    mark_options = {
                        scale = _ms(),
                    },
                 },
                 Coordinates(
                    dram_sizes,
                    getname(d.runs, :bytes_dram_input_tensors) ./ getname(d.runs, :bytes_input_tensors)
                )
            ),
        ))
        push!(plots, @pgf(LegendEntry("$formulation: input tensors")))
        push!(plots, @pgf(PlotInc(
                 {
                    thick,
                    mark_options = {
                        scale = _ms(),
                    },
                 },
                 Coordinates(
                    dram_sizes,
                    getname(d.runs, :bytes_dram_output_tensors) ./ getname(d.runs, :bytes_output_tensors)
                 )
            ),
        ))
        push!(plots, @pgf(LegendEntry("$formulation: output tensors")))
    end

    plot = TikzDocument()
    push!(plot, vasymptote())
    width, height = _wh()

    axs = @pgf Axis(
        {
            grid = "major",
            width = width,
            height = height,
            xlabel = "DRAM Limit (GB)",
            ylabel = "Percent of Kernel IO in DRAM",
            # put legend on the bottom right
            legend_style = {
                at = Coordinate(1.0, 0.0),
                anchor = "south east",
            },
            vasymptote = first(data).default_alloc_size[] ./ 1E9,
        },
        plots...
    )

    push!(plot, TikzPicture(axs))

    empty!(PGFPlotsX.CUSTOM_PREAMBLE)
    push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{colorbrewer}")

    pgfsave(file, plot)
    return plot
end

function pgf_plot_performance(f, cache, suffix;
        file = "plot.tex",
        formulations = ("static", "synchronous",),
        filter_function = (x,y) -> (x,y),
        ylabel = "Slow Down Relative\\\\to all DRAM",
        xlabel = "DRAM Limit (GB)",
    )

    paths = canonical_path.(
        Ref(f),
        formulations,
        Ref(cache),
        Ref(nGraph.Backend("CPU")),
        Ref(suffix),
    )
    data = deserialize.(paths)

    pmm_performance = get_pmm_performance(data)
    dram_performance = get_dram_performance(data)

    xhi = typemin(Float64)
    yhi = typemin(Float64)

    plots = []
    for (datum, formulation) in zip(data, formulations)
        io_size = datum.io_size[]
        dram_sizes = (getname(datum.runs, :dram_limit) ./ 1E3) .+ (io_size ./ 1E9)

        x = dram_sizes
        y = getname(datum.runs, :actual_runtime) ./ dram_performance

        # Update bounds on plot
        xhi = max(xhi, maximum(x))
        yhi = max(yhi, maximum(y))

        x, y = filter_function(x, y)

        append!(plots, [
            @pgf(PlotInc(
                {
                   thick,
                   mark_options = {
                       scale = _ms(),
                   },
                },
                Coordinates(x, y))
            ),
            @pgf(LegendEntry(formulation))
        ])
    end

    @show xhi
    @show yhi

    plot = TikzDocument()
    width, height = _wh()
    ymax = ceil(Int, yhi)
    axs = @pgf Axis(
        {
            width = width,
            height = height,
            xmin = -10,
            xmax = 20 * ceil(Int, xhi / 20),
            ymin = 0.5,
            ymax = ymax,
            ytick = 1:ymax,
            ylabel_style = {
                rotate = -90,
                align = "center",
            },
            xlabel_style = {
                align = "center",
            },
            grid = "major",
            xlabel = xlabel,
            ylabel = ylabel,
            # put legend on the bottom right
            legend_style = {
                at = Coordinate(1.0, 1.0),
                anchor = "north east",
            },
        },
        plots...
    )

    push!(plot, TikzPicture(axs))

    empty!(PGFPlotsX.CUSTOM_PREAMBLE)
    push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{colorbrewer}")

    pgfsave(file, plot)
    return plot
end

#####
##### Speedup over all PMM
#####

# Plot ratios of PMEM to DRAM on the x-axis.
function pgf_speedup(fns, ratios::Vector{<:Rational}, cache;
        file = joinpath(FIGDIR, "speedup.pdf"),
        formulations = ("numa", "static", "synchronous")
    )

    doc = TikzDocument()
    push!(doc, hasymptote())

    count = 1
    axes = map(fns) do f
        axs = _speedup(f, ratios, cache, formulations;
            legend_entries = (count == 1),
            top = (count == 1),
            bot = (count == length(fns)),
            left_label = in(count, (1, 3)),
            bottom_label = in(count, (3, 4)),
        )
        count += 1
        return axs
    end

    gp = @pgf GroupPlot(
        {
            height = "4cm",
            width = "6cm",
            group_style = {
                group_size = "2 by 2",
                vertical_sep = "1.5cm",
            },
        },
        axes...
    )

    pic = TikzPicture(gp)

    # Post processing for legends and stuff
    push!(pic, raw"\path (top)--(bot) coordinate[midway] (group center);")
    push!(pic, raw"\node[below=4.5cm] at(group center -| current bounding box.south) {\pgfplotslegendfromname{legend}};")

    push!(doc, pic)
    pgfsave(file, doc)
    return doc
end

const _speedup_formulation_lut = Dict(
    "static" => "static-AutoTM",
    "synchronous" => "sync-AutoTM",
)

# Inner plot for each `f`
function _speedup(f, ratios, cache, formulations;
        legend_entries = false,
        top = false,
        bot = false,
        left_label = false,
        bottom_label = false,
    )
    paths = canonical_path.(Ref(f), formulations, Ref(cache), Ref(nGraph.Backend("CPU")))
    data = deserialize.(paths)

    pmm_performance = get_pmm_performance(data)
    dram_performance = get_dram_performance(data)

    plots = []
    for (datum, formulation) in zip(data, formulations)
        # This x and y data point
        x = []
        y = []

        for ratio in ratios
            # Check if we have the "ratio" key
            runs = datum.runs
            if haskey(first(runs), :ratio)
                ind = findfirst(x -> x[:ratio] == ratio, runs)
            else
                @info "Performing Ratio Search Fallback"
                ind = findabsmin(x -> compare_ratio(getratio(x), ratio), datum.runs)
            end

            perf = pmm_performance / datum.runs[ind][:actual_runtime]
            push!(x, ratio_string(ratio))
            push!(y, perf)
        end

        # Emit the plot for this series.
        push!(plots,
            @pgf(PlotInc(
                Coordinates(x, y),
            ))
        )

        if legend_entries
            push!(
                plots,
                @pgf(LegendEntry(get(_speedup_formulation_lut, formulation, formulation)))
            )
        end
    end

    symbolic_coords = ratio_string.(ratios)
    dline = pmm_performance / dram_performance

    options = []
    if legend_entries
        push!(options, "legend to name = legend")
    end

    if left_label
        push!(options, "ylabel = Speedup over all PMM")
    end

    if bottom_label
        push!(options, "xlabel = PMM to DRAM Ratio")
    end

    axs = @pgf Axis(
        {
            title = AutoTM.Experiments.titlename(f),
            enlarge_x_limits=0.30,
            ybar,
            ymin=0,
            symbolic_x_coords = symbolic_coords,
            nodes_near_coords_align={vertical},
            ymajorgrids,
            xtick="data",
            ytick = 1:(ceil(Int, pmm_performance / dram_performance)+1),
            ymax = ceil(Int, pmm_performance / dram_performance),
            options...,
            legend_style = {
                legend_columns = -1,
            },
        },
        plots...,
        # Line for All DRAM performance
        hline(dline;
              xl = first(symbolic_coords),
              xu = last(symbolic_coords)
             ),
    )

    if legend_entries
        push!(axs, raw"\addlegendimage{line legend, red, densely dashed, ultra thick}")
        push!(axs, @pgf(LegendEntry("All DRAM")))
    end

    if top
        push!(axs, "\\coordinate (top) at (rel axis cs:0,1);")
    end

    if bot
        push!(axs, "\\coordinate (bot) at (rel axis cs:1,0);")
    end

    return axs
end

#####
##### Performance of AutoTM for the large workloads.
#####

function pgf_large_performance(fns, cache, cache_2lm;
        file = joinpath(FIGDIR, "large_performance.tex"),
        formulations = ("static", "synchronous"),
        ylabel = "Speedup over 2LM",
    )

    # Step through each function, get the 2lm performance
    bar_plots = []
    coords = String[]
    backend = nGraph.Backend("CPU")

    baseline_runtime = Dict{Any,Float64}()
    for f in fns
        baseline = deserialize(canonical_path(f, "2lm", cache_2lm, backend))
        baseline_runtime[f] = minimum(getname(baseline.runs, :actual_runtime))

        push!(coords, "$(AutoTM.Experiments.titlename(f)) ($(f.batchsize))")
    end

    @show coords

    for formulation in formulations
        x = []
        y = []

        for (i, f) in enumerate(fns)
            data = deserialize(canonical_path(f, formulation, cache, backend))
            speedup = baseline_runtime[f] / minimum(getname(data.runs, :actual_runtime))
            push!(x, "$(AutoTM.Experiments.titlename(f)) ($(f.batchsize))")
            push!(y, speedup)
        end

        append!(bar_plots, [
            @pgf(PlotInc(
                Coordinates(x, y),
            ))
            @pgf(LegendEntry(get(_speedup_formulation_lut, formulation, formulation)))
        ])
    end

    plt = TikzDocument()
    push!(plt, """
    \\pgfplotsset{width=10cm,height=4cm}
    """)

    # Bar axis
    tikz = TikzPicture()
    axs = @pgf Axis(
        {
            ybar,
            bar_width = "18pt",
            enlarge_x_limits=0.20,
            symbolic_x_coords = coords,
            legend_style =
            {
                 at = Coordinate(0.05, 1.05),
                 anchor = "south west",
                 legend_columns = -1,
            },
            ymin=0,
            ymajorgrids,
            ylabel_style={
                align = "center",
                rotate = -90,
            },
            xticklabel_style={
                rotate = 10,
            },
            xtick = "data",

            # Lables
            ylabel = ylabel,
        },
        bar_plots...,
    )
    push!(tikz, axs)
    push!(plt, tikz)

    pgfsave(file, plt)
    return plt
end

#####
##### Plots for GPU code
#####

function pgf_gpu_performance_plot(
        funcs, cache;
        file = joinpath(FIGDIR, "gpu_performance.pdf"),
        formulations = ("synchronous", "asynchronous")
    )

    backend = nGraph.GPU

    bar_plots = []
    memory_plots = []
    annotations = []

    coords = []
    batchsize_offset = -0.5
    network_offset = -1.5

    # First step, compute the managed runtime performance for each function as a baseline
    baseline_runtime = Dict{Any,Float64}()
    for (i, f) in enumerate(funcs)
        for formulation in formulations
            data = deserialize(canonical_path(f, formulation, cache, backend))
            baseline_runtime[f] = min(
                get(baseline_runtime, f, typemax(Float64)),
                data.gpu_managed_runtime[]
            )
        end
        push!(annotations, comment(i, batchsize_offset, f.batchsize))
    end

    # Generate function labels
    for title in unique(AutoTM.Experiments.titlename.(funcs))
        # Get the start and stop indices for this title
        start = findfirst(x -> AutoTM.Experiments.titlename(x) == title, funcs)
        stop = findlast(x -> AutoTM.Experiments.titlename(x) == title, funcs)
        center = (start + stop) / 2

        push!(annotations, comment(center, network_offset, title))
    end

    # Generate data for the formulations
    for formulation in formulations
        y = []
        x = []

        for (i, f) in enumerate(funcs)
            @show name(f)
            data = deserialize(canonical_path(f, formulation, cache, backend))
            speedup = baseline_runtime[f] / minimum(getname(data.runs, :actual_runtime))

            println("Runtime - $formulation - $f:", minimum(getname(data.runs, :actual_runtime)))

            push!(x, i)
            push!(y, speedup)

        end
        append!(bar_plots, [
            @pgf(PlotInc(
                Coordinates(x, y),
            ))
            @pgf(LegendEntry(formulation))
        ])
    end

    # Finally, plot the ideal implmenetation
    x = []
    y = []
    for (i, f) in enumerate(funcs)
        data = deserialize.(canonical_path.(Ref(f), formulations, Ref(cache), Ref(backend)))

        # Conovert milliseconds to seconds
        ideal = minimum(minimum.(getname.(getname(data, :runs), :oracle_time))) / 1E6
        @show ideal

        push!(x, i)
        push!(y, baseline_runtime[f] / ideal)
    end

    append!(bar_plots, [
        @pgf(PlotInc(
            Coordinates(x, y),
        ))
        @pgf(LegendEntry("oracle"))
    ])

    # Generate the lefthand bar plot
    plt = TikzDocument()
    push!(plt, """
    \\pgfplotsset{
        width=12cm,
        height=5cm
    }
    """)

    tikz = TikzPicture()
    axs = @pgf Axis(
        {
            ybar,
            xmin = 0.5,
            xmax = length(funcs) + 0.5,
            xmajorticks = "false",
            bar_width = "5pt",
            clip = "false",
            legend_style =
            {
                 at = Coordinate(0.05, 1.05),
                 anchor = "south west",
                 legend_columns = -1
            },
            ymin=0,
            ymajorgrids,
            ytick = 1:8,
            ylabel_style={
                align = "center",
            },

            # Lables
            ylabel = "Speedup over\\\\CudaMallocManaged",
        },
        bar_plots...,
        annotations...,
    )
    push!(tikz, axs)
    push!(plt, tikz)

    pgfsave(file, plt)
    return plt
end

#####
##### The front plot.
#####

iterations_per_hour(runtime) = 3600 / (runtime)

const LEGEND_LUT = Dict(
    "synchronous" => "\\textbf{AutoTM}",
    "numa" => "NUMA",
)
_legend(x) = get(LEGEND_LUT, x, x)

function plot_front(f, ratio::Rational, cache;
        formulations = ("numa", "synchronous"),
        file = joinpath(FIGDIR, "front.pdf")
    )

    bar_plots = []

    backend = nGraph.Backend("CPU")
    data = deserialize.(canonical_path.(Ref(f), formulations, Ref(cache), Ref(backend)))

    pmm_iterations = iterations_per_hour(get_pmm_performance(data))
    dram_iterations = iterations_per_hour(get_dram_performance(data))

    @show pmm_iterations
    @show dram_iterations

    tick_offset = -10
    horizontal = 1
    barwidth = "20pt"

    # Add the PMM performance
    append!(bar_plots, [
        @pgf(PlotInc(
            {
                bar_shift = 0,
            },
            Coordinates([horizontal], [pmm_iterations]))
        ),
        comment(horizontal, tick_offset, "All PMM"),
    ])

    # Add performance for the inner formuliations
    horizontal += 1
    for formulation in formulations
        x = []
        y = []

        data = deserialize(canonical_path(f, formulation, cache, backend))
        ind = findabsmin(x -> compare_ratio(getratio(x), ratio), data.runs)
        perf = iterations_per_hour(data.runs[ind][:actual_runtime])

        @show formulation
        @show perf

        append!(bar_plots, [
            @pgf(PlotInc(
                {
                    bar_shift = 0,
                },
                Coordinates([horizontal], [perf]))
            ),
            comment(horizontal, tick_offset, _legend(formulation)),
        ])
        horizontal += 1
    end

    # Add DRAM performance
    append!(bar_plots, [
        @pgf(PlotInc(
            {
               bar_shift = 0,
            },
            Coordinates([horizontal], [dram_iterations]))
        ),
        comment(horizontal, tick_offset, "All DRAM")
    ])

    plt = TikzDocument()
    push!(plt, """
    \\pgfplotsset{width=12cm,height=5cm}
    """)

    #xticklabels = collect(titlename.(fns))
    #@show xticklabels

    # Bar axis
    comment_level = 110
    tikz = TikzPicture()
    font = "\\scriptsize"
    ymax = 130

    axs = @pgf Axis(
        {
            ybar,
            bar_width = barwidth,
            clip = "false",
            xmajorticks = "false",
            enlarge_x_limits=0.20,

            legend_style =
            {
                 at = Coordinate(0.15, 1.05),
                 anchor = "south west",
                 legend_columns = -1
            },
            ymin=0,
            ymax = ymax,
            ymajorgrids,
            ytick = [0, 40, 80],
            ylabel_style={
                align = "center",
            },

            # Lables
            ylabel = "Iterations Per Hour",
        },
        bar_plots...,
        vline(1.5, yl = 0, yu = ymax, color = "black"),
        vline(1.5 + length(formulations), yl = 0, yu = ymax, color = "black"),

        # System Configuration Comments
        comment(
            1,
            comment_level,
            "System with\\\\160 GB PMM\\\\\\textbf{(lowest cost)}";
            font = font,
        ),
        comment(
            1 + (length(formulations) + 1) / 2,
            comment_level,
            "System with 128 GB PMM\\\\and 32 GB DRAM\\\\\\textbf{(mid cost)}";
            font = font,
        ),
        comment(
            2 + length(formulations),
            comment_level,
            "System with\\\\160 GB DRAM\\\\\\textbf{(highest cost)}";
            font = font,
        )
    )
    push!(tikz, axs)
    push!(plt, tikz)

    pgfsave(file, plt)
    return plt
end

#####
##### Error from profiling
#####

function pgf_error_plot(fns, ratios, caches;
        file = joinpath(FIGDIR, "error.pdf"),
        formulations = ("static", "synchronous"),
        # Older versions of the ILP formulation measured time in micro-seconds instead of
        # 100ths of a second.
        #
        # This flag adjusts that.
        legacy = false,
        suffix = nothing
    )

    plots = []
    backend = nGraph.Backend("CPU")
    for formulation in formulations
        for f in fns
            for cache in caches
                x = []
                y = []
                data = deserialize(canonical_path(f, formulation, cache, backend, suffix))

                for ratio in ratios
                    ind = findabsmin(x -> compare_ratio(getratio(x), ratio), data.runs)
                    actual_runtime = data.runs[ind][:actual_runtime]

                    if legacy
                        predicted_runtime = data.runs[ind][:predicted_runtime] / 1E6
                    else
                        predicted_runtime = data.runs[ind][:predicted_runtime] / 1E2
                    end

                    push!(x, ratio_string(ratio))
                    push!(y, 100 * (predicted_runtime - actual_runtime) / actual_runtime)
                end

                append!(plots, [
                    @pgf(PlotInc(Coordinates(x, y)))
                    @pgf(LegendEntry("$(AutoTM.Experiments.titlename(f)) - $formulation"))
                ])
            end
        end
    end

    # Make the plot itself
    plt = TikzDocument()
    push!(plt, """
    \\pgfplotsset{
        width=8cm,
        height=4cm}
    """)

    axs = @pgf Axis(
        {
            ybar,
            bar_width = "2pt",
            grid = "major",
            xlabel = "DRAM Limit (GB)",
            ylabel = "Relative Predicted\\\\Runtime Error \\%",
            ylabel_style={
                align = "center",
            },
            # Put the legend outside on the right
            legend_style = {
                at = Coordinate(0.50, 1.20),
                anchor = "south",
                legend_columns = 2,
            },
            legend_cell_align = {"left"},

            # Setup x coordinates
            symbolic_x_coords = ratio_string.(ratios),
            xtick = "data",
        },
        plots...
    )

    push!(plt, TikzPicture(axs))
    pgfsave(file, plt)

    return plt
end

#####
##### Price Performance across a range of DRAM/PMM costs.
#####

function pgf_price_performance(
        pairs::Vector{<:Pair},
        ratios::Vector{<:Rational},
        cache;

        base_cost = 0,
        dram_cost_per_gb = 0,
        pmm_cost_per_gb = 0,
        as_ratio = false,
        # Can be `nothing` if no file to be emitted
        file = "plot.tex",
        width = 7,
        height = 4,
    )

    # Initialize a vector of PGFPlotsX structs that will be used to generate the actual
    # plot.
    plots = []
    backend = nGraph.Backend("CPU")

    for (f, formulation) in pairs
        #data = load_save_files(f, formulation)
        data = deserialize(canonical_path(f, formulation, cache, backend))
        dram_performance = get_dram_performance(data)

        # Compute the dram cost
        all_dram_cost = base_cost +
            (dram_cost_per_gb * first(data.runs)[:dram_alloc_size] / 1E9) +
            (pmm_cost_per_gb * first(data.runs)[:pmem_alloc_size] / 1E9)

        @show first(data.runs)[:dram_alloc_size] / 1E9
        @show all_dram_cost

        # This x and y data point
        x = []
        y = []

        for (i, ratio) in enumerate(ratios)
            ind = findabsmin(x -> compare_ratio(getratio(x), ratio), data.runs)
            perf = dram_performance / data.runs[ind][:actual_runtime]
            pmm_amount = data.runs[ind][:pmem_alloc_size]
            dram_amount = data.runs[ind][:dram_alloc_size]
            cost = base_cost +
                (dram_cost_per_gb * dram_amount / 1E9) +
                (pmm_cost_per_gb * pmm_amount / 1E9)

            push!(x, cost / all_dram_cost)
            push!(y, perf)
        end

        # Emit the plot for this series.
        append!(plots, [
            @pgf(PlotInc(
                Coordinates(x, y),
            )),
            @pgf(LegendEntry(replace(AutoTM.Experiments.titlename(f), "_" => " "))),
        ])
    end

    plt = TikzDocument()
    push!(plt, """
    \\pgfplotsset{width=$(width)cm,height=$(height)cm}
    """)

    symbolic_coords = ratio_string.(ratios)

    push!(plots, @pgf(PlotInc(Coordinates([0, 1], [0, 1]))))

    # Bar axis
    tikz = TikzPicture()
    axs = @pgf Axis(
        {
            #ybar,
            #bar_width = "5pt",
            legend_style =
            {
                 at = Coordinate(1.05, 1.05),
                 anchor = "north west",
                 legend_columns = 1,
            },
            ymin = 0,
            #symbolic_x_coords = symbolic_coords,
            nodes_near_coords_align={vertical},
            ymajorgrids,
            xmajorgrids,
            ylabel_style={
                align = "center",
            },
            title = "Base Cost: $base_cost. DRAM Cost: $dram_cost_per_gb. PMM Cost: $pmm_cost_per_gb",
            #ytick = [0,1],
            # Lables
            xlabel = "Cost Ratio to all DRAM",
            ylabel = "Perforamnce Ratio to all DRAM",

        },
        plots...,
    )
    push!(tikz, axs)
    push!(plt, tikz)

    isnothing(file) || pgfsave(file, plt)
    return plt
end

#####
##### The cost-performance plot from the ASPLOS 20 submission.
#####

"""
`pairs`: Vector of Pairs, first element is a model, second element is a formulation string.
"""
function pgf_cost(
        pairs::Vector{<:Pair},
        ratios::Vector{<:Rational},
        cache;
        # Can either define cost_ratio or the costs of PMM and DRAM separately
        cost_ratio = 2.1,
        # location of output file - if `nothing`, no file will be emitted
        file = joinpath(FIGDIR, "costs.pdf"),
    )

    backend = nGraph.Backend("CPU")
    plots = []
    for (f, formulation) in pairs
        data = deserialize(canonical_path(f, formulation, cache, backend))
        dram_performance = get_dram_performance(data)

        # This x and y data point
        x = []
        y = []

        for ratio in ratios
            ind = findabsmin(x -> compare_ratio(getratio(x), ratio), data.runs)
            perf = dram_performance / data.runs[ind][:actual_runtime]

            push!(x, ratio_string(ratio))
            push!(y, perf)
        end

        # Emit the plot for this series.
        append!(plots, [
            @pgf(PlotInc(
                Coordinates(x, y),
            )),
            @pgf(LegendEntry(replace(AutoTM.Experiments.titlename(f), "_" => " "))),
        ])
    end

    plt = TikzDocument()
    push!(plt, """
    \\pgfplotsset{width=7cm,height=4cm}
    """)

    symbolic_coords = ratio_string.(ratios)

    # Bar axis
    tikz = TikzPicture()
    axs = @pgf Axis(
        {
            "axis_y_line*=left",
            "scale_only_axis",
            ybar,
            bar_width = "5pt",
            legend_style =
            {
                 at = Coordinate(0.05, 1.05),
                 anchor = "south west",
                 legend_columns = 2,
            },
            ymin = 0,
            ymax = 1,
            symbolic_x_coords = symbolic_coords,
            nodes_near_coords_align={vertical},
            ymajorgrids,
            ylabel_style={
                align = "center",
            },
            xtick="data",
            #ytick = [0,1],
            # Lables
            xlabel = "PMM to DRAM Ratio",
            ylabel = "Performance Relative to all DRAM",

        },
        plots...,
    )
    push!(tikz, axs)

    # Cost Axis
    y_cost = map(ratios) do ratio
        # PMM
        num = ratio.num

        # DRAM
        den = ratio.den

        # Total approximately sums to 1
        return (num / cost_ratio + den) / (num + den)
    end

    axs = @pgf Axis(
        {
            "axis_y_line*=right",
            "scale_only_axis",
            axis_x_line = "none",
            only_marks,
            ymin = 0,
            ymax = 1,
            symbolic_x_coords = symbolic_coords,
            ylabel = "Memory cost\\\\relative to all DRAM"
        },
        PlotInc(
            {
                mark = "text",
                text_mark = raw"\$",
                mark_options = {
                    color = "black",
                    scale = 1.5,
                },
            },
            Coordinates(symbolic_coords, y_cost),
        ),
    )
    push!(tikz, axs)
    push!(plt, tikz)

    !isnothing(file) && pgfsave(file, plt)
    return plt
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using Redpitaya
using PyPlot

rp = RedPitaya("192.168.1.26")

USE_TRIGGER_VERSION = true

freqs  = linspace( 1000, 1e6, 400)
for freq in freqs
  dec = optimalDecimation(rp,freq)

  numPeriods = 10
  freqR = roundFreq(rp,dec,freq)
  numSampPerPeriod = numSamplesPerPeriod(rp,dec,freqR)
  numSamp = numSampPerPeriod*numPeriods

  println("Frequency = $freqR Hz")
  println("Number Sampling Points per Period: $numSampPerPeriod")

  # start sending
  send(rp,"GEN:RST")
  if !USE_TRIGGER_VERSION
    sendAnalogSignal(rp,1,"SINE",freqR,0.4)
  else
    sendAnalogSignal(rp,1,"SINE",freqR,0.4, numPeriods*2)
  end

  # receive data
  if !USE_TRIGGER_VERSION
    u1 = receiveAnalogSignal(rp, 1, 0, numSamp, dec=dec, delay=0.2)
  else
    trigger = "AWG_PE" #"CH1_PE"
    u1 = receiveAnalogSignalWithTrigger(rp, 1, 0, numSamp, dec=dec,
                delay=0.01, typ="OLD", trigger=trigger, triggerLevel=-0.1,
                triggerDelay=numSampPerPeriod)

  end

  figure(1)
  clf()
  subplot(2,1,1)
  plot(u1)
  subplot(2,1,2)
  semilogy(abs(rfft(u1))[1:numPeriods*2],"o-b",lw=2)
end