text:
using Random, Turing, PGFPlotsX, NNlib, SpecialFunctions, LaTeXStrings
using StatsFuns: logistic, norminvcdf

using Memoization
using ReverseDiff

using Zygote

N = 50
X = vcat(-rand(N) .- 0.5, rand(N) .+ 0.5)
l = vcat(ones(N), zeros(N))
y = l + 0.1*randn(2*N)

xt = range(-5., 5., length = 100);

const triangleconst = π/(2*sqrt(2))
const preluconst = π/4
_trianglewave(x) = (x - π * floor(x/π + 0.5)) * (-1)^floor(x/π + 0.5)
_trianglewave(x::ReverseDiff.TrackedReal) = ReverseDiff.ForwardOptimize(_trianglewave)(x)

heaviside(x) = x < zero(x) ? zero(x) : x >= zero(x) ? one(x) : oftype(x,0.5)
lsmatern(x, q, ν, λ) = q*heaviside(x)*x^(ν-1/2)*exp(-λ*x)
heaviside(x::ReverseDiff.TrackedReal) = ReverseDiff.ForwardOptimize(heaviside)(x)

abstract type PeriodicFunction end

struct SinActivation <: PeriodicFunction end
const sineconst = sqrt(2)
(σ::SinActivation)(x::T) where T = T(sineconst) * sin(x)

struct SinCosActivation <: PeriodicFunction end
(σ::SinCosActivation)(x) = sin(x) + cos(x)

struct TriangleWave <: PeriodicFunction end
(σ::TriangleWave)(x) = triangleconst * _trianglewave(x)

struct PeriodicReLU <: PeriodicFunction end
(σ::PeriodicReLU)(x) = preluconst * (_trianglewave(x) + _trianglewave(x + π/2))

abstract type ActivationFunction end

struct ReLU <: ActivationFunction end
(σ::ReLU)(x) = relu(x)
prior_weights(::ReLU) = Normal()
prior_bias(::ReLU) = Normal()
getname(::ReLU) = "relu"

struct Step <: ActivationFunction end
(σ::Step)(x) = heaviside(x)
prior_weights(::Step) = Normal()
prior_bias(::Step) = Normal()
getname(::Step) = "step"

struct Sigmoid <: ActivationFunction end
(σ::Sigmoid)(x) = sigmoid(x)
prior_weights(::Sigmoid) = Normal()
prior_bias(::Sigmoid) = Normal()
getname(::Sigmoid) = "sigmoid"

struct RBFNN <: ActivationFunction end
(σ::RBFNN)(x) = exp(-x^2)
prior_weights(::RBFNN) = Normal()
prior_bias(::RBFNN) = Normal()
getname(::RBFNN) = "RBF_NN"

struct MaternNN <: ActivationFunction
    q # precomputed constant
    λ # precomputed constant
    ν # smoothness
    function MaternNN(ν)
        λf(ν) = sqrt(2*ν)
        Af(ν) = sqrt(2 * π^(1/2) * λf(ν)^(2*ν) * gamma(ν + 0.5) / gamma(ν)) / gamma(ν + 0.5)
        return new(Af(ν), λf(ν), ν)
    end
end
(σ::MaternNN)(x) = σ.q*heaviside(x)*x^(σ.ν-1/2)*exp(-σ.λ*x)
prior_weights(::MaternNN) = Normal()
prior_bias(::MaternNN) = Normal()
getname(σ::MaternNN) = "Matern_NN_$(Int(σ.ν*10))"

struct RBF <: ActivationFunction
    σ::PeriodicFunction
end
(σ::RBF)(x) = σ.σ(x)
prior_weights(::RBF) = Normal()
prior_bias(::RBF) = Uniform(-π, π)
getname(::RBF) = "RBF"

struct Matern <: ActivationFunction
    σ::PeriodicFunction
    ν # smoothness
end
(σ::Matern)(x) = σ.σ(x)
prior_weights(σ::Matern) = σ.ν == 0.5 ? Cauchy() : TDist(2*σ.ν)
prior_bias(::Matern) = Uniform(-π, π)
getname(σ::Matern) = "Matern-$(Int(σ.ν*10))"

@pgf p = Axis(
    {
        ylabel= raw"$\sigma(x)$",
        no_markers,
        xticklabels={},
        yticklabels={},
        width = "15cm",
    },
    PlotInc(Table(xt, RBF(SinActivation()).(xt))),
    LegendEntry("sinusiod"),
    PlotInc(Table(xt, RBF(SinCosActivation()).(xt))),
    LegendEntry("sine-cosine"),
    PlotInc(Table(xt, RBF(TriangleWave()).(xt))),
    LegendEntry("triangle wave"),
    PlotInc(Table(xt, RBF(PeriodicReLU()).(xt))),
    LegendEntry("periodic ReLU"),
)

display("image/svg+xml", p)

@model function bnn(x, y, afun, K)

    w ~ filldist(prior_weights(afun), 1, K)
    b ~ filldist(prior_bias(afun), 1, K)

    l ~ FlatPos(0)
    z = afun.(l * x * w .+ b)

    wo ~ filldist(Normal(), K)
    bo ~ Normal()

    f = (z * (wo / K) .+ bo) / 2

    s ~ FlatPos(0)
    Turing.@addlogprob! -sum(abs2((yj-fj)/s)/2 + log(s) for (fj,yj) in zip(f,y))
end

Turing.setadbackend(:reversediff)
Turing.setrdcache(true);
Turing.setprogress!(true);

function predict(x, afun, w::Chains, b::Chains, l, wo::Chains, bo, s; burnin = 5000)
    nsamples = length(w)

    latents = map(s -> begin
            w_ = Array(w[s])
            b_ = Array(b[s])

            l_ = l[s]

            wo_ = vec(Array(wo[s]))
            bo_ = bo[s]

            z = afun.(x * l_ * w_ .+ b_)
            f = z * (wo_ / length(wo_)) .+ bo_
            f ./= 2
        end, burnin:nsamples)

    err_o = 1.96*sqrt.(var(latents)) .+ s
    err = 1.96*std(latents)

    return mean(latents), err, err_o
end

function plot(m, sf, so)

    @pgf p = Axis(
        {
            "scatter/classes" = {
                0 = {"black", mark_size="0.5"},
                1 = {"black", mark_size="0.5"},
            },
            axis_lines = "none",
            ymin = minimum(m .- so) - 0.25,
            ymax = maximum(m .+ so) + 0.25
        },
        PlotInc( { no_markers, white, solid }, Table(xt, m) ),
        PlotInc( { no_markers, white, solid }, Table(xt, m .+ so) ),
        PlotInc( { no_markers, white, solid }, Table(xt, m .- so) ),

        PlotInc( { scatter, "only marks", "scatter src" = "explicit symbolic", },
            Table({ meta = "label" }, x = X, y = y, label = Int.(l)) )
        )

    return p
end

function plot_lines(m, sf, so, filename)

        @pgf p_lines = Axis(
        {
            "scatter/classes" = {
                0 = {"black", mark_size="0.1"},
                1 = {"black", mark_size="0.1"},
            },
            ymin = minimum(m .- so) - 0.2,
            ymax = maximum(m .+ so) + 0.2
        },

        PlotInc( Graphics({xmin=minimum(xt),xmax=maximum(xt), ymin=minimum(m .- so), ymax=maximum(m .+ so)}, "$(pwd())/$(filename)") ),
        PlotInc( { no_markers, gray, dotted }, Table(xt, ones(length(xt))*0.5)),
        PlotInc( { no_markers, red, solid }, Table(xt, m) ),
        PlotInc( { no_markers, black, dashed }, Table(xt, m .+ sf) ),
        PlotInc( { no_markers, black, dashed }, Table(xt, m .- sf) ),
        PlotInc( { no_markers, black, solid }, Table(xt, m .+ so) ),
        PlotInc( { no_markers, black, solid }, Table(xt, m .- so) ),
        )

    return p_lines
end

afuns = [
    ReLU(),
    Step(),
    Sigmoid(),
    Matern(SinActivation(), 1/2),
    Matern(SinActivation(), 3/2),
    RBF(SinActivation()),
    MaternNN(1/2, 1),
    MaternNN(3/2, 1),
    RBFNN()
]

K = 10
for afun in afuns#[k:k]

    @info afun

    model = bnn(X, y, afun, K);
    chain = sample(model, NUTS(1000, 0.75), 5000);

    @info "Finished sampling"

    # extract names

    w_names = filter(n -> contains(string(n), "w["), names(chain))
    wo_names = filter(n -> contains(string(n), "wo["), names(chain))
    b_names = filter(n -> contains(string(n), "b["), names(chain))

    m, sf, so = predict(xt, afun,
        chain[w_names],
        chain[b_names],
        chain["l"],
        chain[wo_names],
        chain["bo"],
        mean(chain["s"]),
        burnin = 1000
    )

    p = plot(m, sf, so)

    @info "Saving plots"

    pgfsave("$(getname(afun))_plot.png", p)

    p = plot_lines(m, sf, so, "$(getname(afun))_plot.png")

    pgfsave("$(getname(afun))_plot.pdf", p)
    pgfsave("$(getname(afun))_plot.tex", p, include_preamble = false)
end

max_stars_repo_name: AaltoML/PeriodicBNN

max_issues_repo_name: AaltoML/PeriodicBNN

max_forks_repo_name: AaltoML/PeriodicBNN


text:
module Blender

using ColorTypes
using PrettyTables

include("bpy.jl")
include("materials.jl")
include("object.jl")
include("render.jl")

end # module

ext: jl


text:
"""Returns a matrix of water depth for the whole domain that contains a
Gaussian seamount in the middle. Water depth, heigth and width of the
seamount are adjusted with the constants water_depth, topofeat_height and topofeat_width."""
function seamount()
    xx_T,yy_T = meshgrid(x_T_halo,y_T_halo)
    bumpx = exp.(-((xx_T .- Lx/2).^2)/(2*topofeat_width^2))
    bumpy = exp.(-((yy_T .- Ly/2).^2)/(2*topofeat_width^2))

    H = water_depth .- topofeat_height*bumpx.*bumpy
    return Numtype.(H)
end

"""Returns a matrix of water depth for the whole domain that contains a
meridional Gaussian ridge in the middle. Water depth, heigth and width of the
ridge are adjusted with the constants water_depth, topofeat_height and topofeat_width."""
function ridge()
    xx_T,yy_T = meshgrid(x_T_halo,y_T_halo)
    bumpx = exp.(-((xx_T .- Lx/2).^2)/(2*topofeat_width^2))

    H = water_depth .- topofeat_height*bumpx
    return Numtype.(H)
end

"""Same as ridge() but for 3 ridges at 1/4,1/2,3/4 of the domain."""
function ridges()
    xx_T,yy_T = meshgrid(x_T_halo,y_T_halo)

    # bumps in x direction
    # shift slightly left/right to avoid a symmetric solution
    b0x = exp.(-(xx_T.^2)/(2*topofeat_width^2))
    b1x = exp.(-((xx_T .- 0.99*Lx/4).^2)/(2*topofeat_width^2))
    b2x = exp.(-((xx_T .- 1.01*Lx/2).^2)/(2*topofeat_width^2))
    b3x = exp.(-((xx_T .- 0.99*3*Lx/4).^2)/(2*topofeat_width^2))
    b4x = exp.(-((xx_T .- Lx).^2)/(2*topofeat_width^2))

    th = topofeat_height    # for convenience
    H = water_depth .- th*b0x .- th*b1x .- th*b2x .- th*b3x .- th*b4x
    return Numtype.(H)
end

"""Returns a matrix of constant water depth specified by the constant water_depth."""
function flat_bottom()
    H = fill(water_depth,(nx+2*haloη,ny+2*haloη))
    return Numtype.(H)
end

# rename for convenience
if topography_feature == "ridge"
    topography = ridge
elseif topography_feature == "ridges"
    topography = ridges
elseif topography_feature == "seamount"
    topography = seamount
elseif topography_feature == "flat"
    topography = flat_bottom
else
    throw(error("Topography feature not correctly declared."))
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
module Helper

read_content(fname::String)::String = read(fname, String)

read_lines(fname::String)::Vector{String} = readlines(fname)

read_lines_as_ints(fname::String)::Vector{Int} = [parse(Int, line) for line in readlines(fname)]

end # module


text:
using Documenter, IteratorInterfaceExtensions

makedocs(
	modules = [IteratorInterfaceExtensions],
	sitename = "IteratorInterfaceExtensions.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/IteratorInterfaceExtensions.jl.git"
)


text:
using DataFrames
using CSV
using Lathe.preprocess:TrainTestSplit
using GLM
using StatsBase
using StatsModels
using ROC
using Plots

data = DataFrame(CSV.File("tacke1000a.csv"))

dataTrain, dataTest = TrainTestSplit(data,.80)

fm = @formula(boja ~ x + y)

logisticRegeressor = glm(fm,dataTrain, Binomial(), ProbitLink())

dataPredictBoja = predict(logisticRegeressor, dataTest)

dataPredictTestClass = zeros(length(dataPredictBoja))

for i in 1:length(dataPredictBoja)
    if dataPredictBoja[i] < 0.5
        dataPredictTestClass[i] = 0;
    else
        dataPredictTestClass[i] = 1;
    end
end

FPtest = 0
FNtest = 0
TPtest = 0
TNtest = 0

for i in 1:length(dataPredictTestClass)
    if dataTest.boja[i] == 0 && dataPredictTestClass[i] == 0
        global TNtest +=1
    elseif dataTest.boja[i] == 1 && dataPredictTestClass[i] == 1
        global TPtest +=1
    elseif dataTest.boja[i] == 1 && dataPredictTestClass[i] == 0
        global FNtest +=1
    elseif dataTest.boja[i] == 0 && dataPredictTestClass[i] == 1
        global FPtest +=1
    end
end

accuracyTest = (TPtest + TNtest) / (TPtest + TNtest + FPtest + FNtest)

sensitivity = TPtest / (TPtest + FNtest)

specificity = TNtest / (TNtest + FPtest)

println("TP = $TPtest, FP = $FPtest, TN = $TNtest, FN = $FNtest")
println("accuracy = $accuracyTest\n sensitivity = $sensitivity\n specificity = $specificity\n")

rocTest = ROC.roc(dataPredictBoja, dataTest.boja, true)

aucTest = AUC(rocTest)

if aucTest > 0.9
    println("Klasifikator je jako dobar!")
elseif aucTest > 0.8
    println("Klasifikator je veoma dobar!")
elseif  aucTest > 0.7
    println("Klasifikator je dosta dobar!")
elseif aucTest > 0.5
    println("Klasifikator je relativno dobar!")
else
    println("Klasifikator je los!")
end

#plot(rocTest, label="ROC")

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
# function that generates the ressource flow matrix the
# sim_data is build on
function ressource_flow_matrix(sim_data::DataFrame,
                               incidents::DataFrame,
                               simulation_capacity::Array{Int64,2})
### sim_start:  epoch (in minutes) of the first incident
### sim_length: epoch (in minutes) ends 5 minutes after the incident_minute
###             of the last incidents in our incidents DataFrame
    sim_start  = minimum(sim_data[:,:incident_minute])
    sim_length = maximum(sim_data[:,:incident_minute])

### ressource_flow: the ressource flow matrix of the sim_data
### ressource_flow[:,:,1] = capcacity at district center
### ressource_flow[:,:,2] = capcacity driving to incident
### ressource_flow[:,:,3] = capcacity at incident
### ressource_flow[:,:,4] = capcacity driving back to district center
### ressource_flow[:,:,5] = capcacity at paperwork
### ressource_flow[:,:,6] = capcacity at patrol
### ressource_flow[:,:,7] = minutes of paperwork in the corresponding district
    ressource_flow = Array{Int64,3}(undef,sim_length - sim_start + 1000,
                                    size(simulation_capacity,2), 7) .= 0
    current_weekhour = incidents[1,:weekhour]
    current_minute   = minute(unix2datetime(incidents[1,:epoch]))
    ressource_flow[1,:,1] = simulation_capacity[current_weekhour,:]
    for i = 2:size(ressource_flow,1)
        if current_minute == 60
            current_minute = 0
            if current_weekhour >= size(simulation_capacity,1)
                current_weekhour = 0
            end
            current_weekhour += 1
            if current_weekhour == 1
                ressource_flow[i,:,1] .=
                @view(simulation_capacity[current_weekhour,:]) -
                @view(simulation_capacity[size(simulation_capacity,1),:])
            else
                ressource_flow[i,:,1] .=
                @view(simulation_capacity[current_weekhour,:]) -
                @view(simulation_capacity[current_weekhour-1,:])
            end
        end
        current_minute +=1
    end
    return ressource_flow
end


text:
"""
Generates a `Float64` from a `MRG32k3a` instance.
"""
function rand(rng::MRG32k3a, ::Type{Float64})
    rand(rng)
end

"""
Generates a `Float32` from a `MRG32k3a` instance.
"""
function rand(rng::MRG32k3a, ::Type{Float32})
    convert(Float32, rand(rng))
end

"""
Generates a `Float16` from a `MRG32k3a` instance.
"""
function rand(rng::MRG32k3a, ::Type{Float16})
    convert(Float16, rand(rng))
end

"""
Generates a `Int64` in the given range from a `MRG32k3a` instance.
"""
function rand(rng::MRG32k3a, r::UnitRange{Int64})
    r.start + convert(Int64, (r.stop - r.start) * rand(rng, Float64))
end


text:
"""
Author: Mohammad Dehghani Ashkezari <mdehghan@uw.edu>

Date: 2020-05-08

Function: implemantation of matching (colocalizing) data sets.
"""

using Dates
include("./rest.jl")

"""
    _atomic_match(match::Match)

Colocalizes the source variable (from source table) with a single target variable (from target table).
The tolerance parameters set the matching boundaries between the source and target data sets.
Returns a dataframe containing the source variable joined with the target variable.
This function is meant to be used internally.

Examples
≡≡≡≡≡≡≡≡≡≡

_atomic_match(
             "uspMatch",
             "tblKM1314_Cobalmins",
             "Me_PseudoCobalamin_Particulate_pM",
             "tblDarwin_Phytoplankton",
             "picoprokaryote",
             "2013-08-11",
             "2013-09-05",
             22.25,
             450.25,
             -159.25,
             -127.75,
             -5,
             305,
             1,
             0.25,
             0.25,
             5
             );
"""
function _atomic_match(
                       spName::String,
                       sourceTable::String,
                       sourceVariable::String,
                       targetTable::String,
                       targetVariable::String,
                       dt1::String,
                       dt2::String,
                       lat1::Number,
                       lat2::Number,
                       lon1::Number,
                       lon2::Number,
                       depth1::Number,
                       depth2::Number,
                       timeTolerance::Number,
                       latTolerance::Number,
                       lonTolerance::Number,
                       depthTolerance::Number
                      )

    statement = """
                EXEC $spName
                    '$sourceTable',
                    '$sourceVariable',
                    '$targetTable',
                    '$targetVariable',
                    '$dt1',
                    '$dt2',
                    '$lat1',
                    '$lat2',
                    '$lon1',
                    '$lon2',
                    '$depth1',
                    '$depth2',
                    '$timeTolerance',
                    '$latTolerance',
                    '$lonTolerance',
                    '$depthTolerance'
                """
    api = API();
    status, response = query(api, statement);
    return response
end

"""
Loops through the target data sets and match them with the source data set according to the the accosiated tolerance parameters.
Returns a compiled dataframe of the source and matched target data sets.
The matching results rely on the spatio-temporal tolerance parameters.
Notice the source has to be a single variable and cannot be a climatological variable.
You may pass empty string ('') as source variable if you only want to get the time and location info from the source table.
The target variables (one or more) are matched with the source variable, if any match exists.
Please note that the number of matching entries between each target variable and the source variable might vary depending on the temporal and spatial resolution of the target variable.

# Arguments
- `spname::String`: stored procedure name that executes the matching logic.
- `sourceTable::String`: table name of the source data set.
- `sourceVariable::String`: the source variable. The target variables are matched (colocalized) with this variable.
- `targetTables::Array`: table names of the target data sets to be matched with the source data.
- `targetVariables::Array`: variable names to be matched with the source variable.
- `dt1::String`: start date or datetime.
- `dt2::String`: end date or datetime.
- `lat1::Number`: start latitude [degree N]. This parameter sets the lower bound of the meridional cut. Note latitude ranges from -90 to 90.
- `lat2::Number`: end latitude [degree N]. This parameter sets the upper bound of the meridional cut. Note latitude ranges from -90 to 90.
- `lon1::Number`: start longitude [degree E]. This parameter sets the lower bound of the zonal cut. Note latitude ranges from -180 to 180.
- `lon2::Number`: end longitude [degree E]. This parameter sets the upper bound of the zonal cut. Note latitude ranges from -180 to 180.
- `depth1::Number`: start depth [m]. This parameter sets the lower bound of the vertical cut. Note depth is a positive number (it is 0 at surface and grows towards ocean floor).
- `depth2::Number`: end depth [m]. This parameter sets the upper bound of the vertical cut. Note depth is a positive number (it is 0 at surface and grows towards ocean floor).
- `timeTolerance::Number`: integer array of temporal tolerance values between pairs of source and target datasets. The size and order of values in this array should match those of targetTables. This parameter is in day units except when the target variable represents monthly climatology data in which case it is in month units. Note fractional values are not supported in the current version.
- `latTolerance::Number`: float array of spatial tolerance values in meridional direction [deg] between pairs of source and target data sets.
- `lonTolerance::Number`: float array of spatial tolerance values in zonal direction [deg] between pairs of source and target data sets.
- `depthTolerance::Number`: float array of spatial tolerance values in vertical direction [m] between pairs of source and target data sets.

"""
function compile(;
                 spname::String,
                 sourceTable::String,
                 sourceVariable::String,
                 targetTables::Array,
                 targetVariables::Array,
                 dt1::String,
                 dt2::String,
                 lat1::Number,
                 lat2::Number,
                 lon1::Number,
                 lon2::Number,
                 depth1::Number,
                 depth2::Number,
                 timeTolerance::Array,
                 latTolerance::Array,
                 lonTolerance::Array,
                 depthTolerance::Array
                )

    function shift_dt(dt::String, delta::Int)
        dtFormat = "yyyy-mm-dd";
        if occursin("T", dt)
            dt = replace(dt, "T" => " ");
            if occursin(".000Z", dt)
                dt = replace(dt, ".000Z" => "");
            end
            dtFormat = "yyyy-mm-dd HH:MM:SS";
        end
        dt = DateTime(dt, dtFormat);
        dt = dt + Dates.Day(delta);
        # TODO: Handle monthly climatology data sets
        shifted = Dates.format(dt, "yyyy-mm-dd HH:MM:SS");
    end

    df = DataFrame();
    for i =  1:length(targetTables)
        data = _atomic_match(
                            spname,
                            sourceTable,
                            sourceVariable,
                            targetTables[i],
                            targetVariables[i],
                            shift_dt(dt1, -timeTolerance[i]),
                            shift_dt(dt2, timeTolerance[i]),
                            lat1 - latTolerance[i],
                            lat2 + latTolerance[i],
                            lon1 - latTolerance[i],
                            lon2 + latTolerance[i],
                            depth1 - depthTolerance[i],
                            depth2 + depthTolerance[i],
                            timeTolerance[i],
                            latTolerance[i],
                            lonTolerance[i],
                            depthTolerance[i]
                            )
        if nrow(data) < 1
            println("$i: No matching entry associated with $(targetVariables[i]).")
            continue
        end
        println("$i: $(targetVariables[i]) matched.")

        if nrow(df) < 1
            df = data;
        elseif (
                df[!, names(df)[1]] == data[!, names(df)[1]] &&
                df.lat == data.lat &&
                df.lon == data.lon
                )

            df[!, Symbol(targetVariables[i])] = data[!, Symbol(targetVariables[i])]
            df[!, Symbol(targetVariables[i]*"_std")] = data[!, Symbol(targetVariables[i]*"_std")]
        else
            println("The matched dataframe corresponding to $(targetVariables[i]) does not have the same size as the first targert variable. Please change the tolerance parameters.")
        end
    end
    return df
end


text:
# Structured matrices
using LinearAlgebra: AbstractTriangular

# Matrix wrapper types that we know are square and are thus potentially invertible. For
# these we can use simpler definitions for `/` and `\`.
const SquareMatrix{T} = Union{Diagonal{T},AbstractTriangular{T}}

function rrule(::typeof(/), A::AbstractMatrix{<:Real}, B::T) where {T<:SquareMatrix{<:Real}}
    Y = A / B
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)
    function slash_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        ∂A = @thunk project_A(Ȳ / B')
        ∂B = @thunk project_B(-Y' * (Ȳ / B'))
        return (NoTangent(), ∂A, ∂B)
    end
    return Y, slash_pullback
end

function rrule(
    ::typeof(\), A::T, B::AbstractVecOrMat{<:Real}
) where {T<:SquareMatrix{<:Real}}
    Y = A \ B
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)
    function backslash_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        ∂A = @thunk project_A(-(A' \ Ȳ) * Y')
        ∂B = @thunk project_B(A' \ Ȳ)
        return NoTangent(), ∂A, ∂B
    end
    return Y, backslash_pullback
end

#####
##### `Diagonal`
#####

_diagview(x::Diagonal) = x.diag
_diagview(x::AbstractMatrix) = view(x, diagind(x))
_diagview(x::Tangent{<:Diagonal}) = x.diag
function ChainRulesCore.rrule(::typeof(sqrt), d::Diagonal)
    y = sqrt(d)
    @assert y isa Diagonal
    function sqrt_pullback(Δ)
        Δ_diag = _diagview(unthunk(Δ))
        return NoTangent(), Diagonal(Δ_diag ./ (2 .* y.diag))
    end
    return y, sqrt_pullback
end

# these functions are defined outside the rrule because otherwise type inference breaks
# see https://github.com/JuliaLang/julia/issues/40990
_Diagonal_pullback(ȳ::AbstractMatrix) = return (NoTangent(), diag(ȳ)) # should we emit a warning here? this shouldn't be called if project works right
_Diagonal_pullback(ȳ::Diagonal) = return (NoTangent(), diag(ȳ))
function _Diagonal_pullback(ȳ::Tangent)
    # TODO: Assert about the primal type in the Tangent, It should be Diagonal
    # infact it should be exactly the type of `Diagonal(d)`
    # but right now Zygote loses primal type information so we can't use it.
    # See https://github.com/FluxML/Zygote.jl/issues/603
    return (NoTangent(), ȳ.diag)
end
_Diagonal_pullback(ȳ::AbstractThunk) = return _Diagonal_pullback(unthunk(ȳ))

function rrule(::Type{<:Diagonal}, d::AbstractVector)
    return Diagonal(d), _Diagonal_pullback
end

function rrule(::typeof(diag), A::AbstractMatrix)
    function diag_pullback(ȳ)
        return (NoTangent(), Diagonal(ȳ))
    end
    return diag(A), diag_pullback
end
if VERSION ≥ v"1.3"
    function rrule(::typeof(diag), A::AbstractMatrix, k::Integer)
        function diag_pullback(ȳ)
            return (NoTangent(), diagm(size(A)..., k => ȳ), NoTangent())
        end
        return diag(A, k), diag_pullback
    end

    function rrule(
        ::typeof(diagm), m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...
    )
        function diagm_pullback(ȳ)
            return (NoTangent(), NoTangent(), NoTangent(), _diagm_back.(kv, Ref(ȳ))...)
        end
        return diagm(m, n, kv...), diagm_pullback
    end
end
function rrule(::typeof(diagm), kv::Pair{<:Integer,<:AbstractVector}...)
    function diagm_pullback(ȳ)
        return (NoTangent(), _diagm_back.(kv, Ref(ȳ))...)
    end
    return diagm(kv...), diagm_pullback
end

function _diagm_back(p, ȳ)
    k, v = p
    d = diag(unthunk(ȳ), k)[1:length(v)] # handle if diagonal was smaller than matrix
    return Tangent{typeof(p)}(; second=d)
end

function rrule(::typeof(*), D::Diagonal{<:Real}, V::AbstractVector{<:Real})
    project_D = ProjectTo(D)
    project_V = ProjectTo(V)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        dD = @thunk(project_D(Diagonal(Ȳ .* V)))
        dV = @thunk(project_V(D * Ȳ))
        return (NoTangent(), dD, dV)
    end
    return D * V, times_pullback
end

#####
##### `Adjoint`
#####

# these functions are defined outside the rrule because otherwise type inference breaks
# see https://github.com/JuliaLang/julia/issues/40990
_Adjoint_mat_pullback(ȳ::Tangent, proj) = (NoTangent(), proj(ȳ.parent))
_Adjoint_mat_pullback(ȳ::AbstractVecOrMat, proj) = (NoTangent(), proj(adjoint(ȳ)))
function _Adjoint_mat_pullback(ȳ::AbstractThunk, proj)
    return _Adjoint_mat_pullback(unthunk(ȳ), proj)
end
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_Adjoint_mat_pullback(ȳ::AbstractZero, proj) = (NoTangent(), proj(ȳ))
function rrule(::Type{<:Adjoint}, A::AbstractMatrix{<:Number})
    project_A = ProjectTo(A)
    Adjoint_mat_pullback(ȳ) = _Adjoint_mat_pullback(ȳ, project_A)
    return Adjoint(A), Adjoint_mat_pullback
end

_Adjoint_vec_pullback(ȳ::Tangent) = (NoTangent(), vec(ȳ.parent))
_Adjoint_vec_pullback(ȳ::AbstractMatrix) = (NoTangent(), vec(adjoint(ȳ)))
_Adjoint_vec_pullback(ȳ::AbstractThunk) = return _Adjoint_vec_pullback(unthunk(ȳ))
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_Adjoint_vec_pullback(ȳ::AbstractZero) = (NoTangent(), ȳ)
function rrule(::Type{<:Adjoint}, A::AbstractVector{<:Number})
    return Adjoint(A), _Adjoint_vec_pullback
end

_adjoint_mat_pullback(ȳ::Tangent, proj) = (NoTangent(), proj(ȳ.parent))
_adjoint_mat_pullback(ȳ::AbstractVecOrMat, proj) = (NoTangent(), proj(adjoint(ȳ)))
function _adjoint_mat_pullback(ȳ::AbstractThunk, proj)
    return _adjoint_mat_pullback(unthunk(ȳ), proj)
end
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_adjoint_mat_pullback(ȳ::AbstractZero, proj) = (NoTangent(), proj(ȳ))
function rrule(::typeof(adjoint), A::AbstractMatrix{<:Number})
    project_A = ProjectTo(A)
    adjoint_mat_pullback(ȳ) = _adjoint_mat_pullback(ȳ, project_A)
    return adjoint(A), adjoint_mat_pullback
end

_adjoint_vec_pullback(ȳ::Tangent) = (NoTangent(), vec(ȳ.parent))
_adjoint_vec_pullback(ȳ::AbstractMatrix) = (NoTangent(), vec(adjoint(ȳ)))
_adjoint_vec_pullback(ȳ::AbstractThunk) = return _adjoint_vec_pullback(unthunk(ȳ))
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_adjoint_vec_pullback(ȳ::AbstractZero) = (NoTangent(), ȳ)
function rrule(::typeof(adjoint), A::AbstractVector{<:Number})
    return adjoint(A), _adjoint_vec_pullback
end

#####
##### `Transpose`
#####

# these functions are defined outside the rrule because otherwise type inference breaks
# see https://github.com/JuliaLang/julia/issues/40990
_Transpose_mat_pullback(ȳ::Tangent, proj) = (NoTangent(), proj(ȳ.parent))
_Transpose_mat_pullback(ȳ::AbstractVecOrMat, proj) = (NoTangent(), proj(Transpose(ȳ)))
function _Transpose_mat_pullback(ȳ::AbstractThunk, proj)
    return _Transpose_mat_pullback(unthunk(ȳ), proj)
end
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_Transpose_mat_pullback(ȳ::AbstractZero, proj) = (NoTangent(), proj(ȳ))
function rrule(::Type{<:Transpose}, A::AbstractMatrix{<:Number})
    project_A = ProjectTo(A)
    Transpose_mat_pullback(ȳ) = _Transpose_mat_pullback(ȳ, project_A)
    return Transpose(A), Transpose_mat_pullback
end

_Transpose_vec_pullback(ȳ::Tangent) = (NoTangent(), vec(ȳ.parent))
_Transpose_vec_pullback(ȳ::AbstractMatrix) = (NoTangent(), vec(Transpose(ȳ)))
_Transpose_vec_pullback(ȳ::AbstractThunk) = return _Transpose_vec_pullback(unthunk(ȳ))
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_Transpose_vec_pullback(ȳ::AbstractZero) = (NoTangent(), ȳ)
function rrule(::Type{<:Transpose}, A::AbstractVector{<:Number})
    return Transpose(A), _Transpose_vec_pullback
end

_transpose_mat_pullback(ȳ::Tangent, proj) = (NoTangent(), proj(ȳ.parent))
_transpose_mat_pullback(ȳ::AbstractVecOrMat, proj) = (NoTangent(), proj(transpose(ȳ)))
function _transpose_mat_pullback(ȳ::AbstractThunk, proj)
    return _transpose_mat_pullback(unthunk(ȳ), proj)
end
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_transpose_mat_pullback(ȳ::AbstractZero, proj) = (NoTangent(), proj(ȳ))
function rrule(::typeof(transpose), A::AbstractMatrix{<:Number})
    project_A = ProjectTo(A)
    transpose_mat_pullback(ȳ) = _transpose_mat_pullback(ȳ, project_A)
    return transpose(A), transpose_mat_pullback
end

_transpose_vec_pullback(ȳ::Tangent) = (NoTangent(), vec(ȳ.parent))
_transpose_vec_pullback(ȳ::AbstractMatrix) = (NoTangent(), vec(transpose(ȳ)))
_transpose_vec_pullback(ȳ::AbstractThunk) = return _transpose_vec_pullback(unthunk(ȳ))
# currently needed by Diffractor (ref https://github.com/JuliaDiff/Diffractor.jl/issues/25)
_transpose_vec_pullback(ȳ::AbstractZero) = (NoTangent(), ȳ)
function rrule(::typeof(transpose), A::AbstractVector{<:Number})
    return transpose(A), _transpose_vec_pullback
end

#####
##### Triangular matrices
#####

function rrule(::Type{<:UpperTriangular}, A::AbstractMatrix)
    function UpperTriangular_pullback(ȳ)
        return (NoTangent(), Matrix(ȳ))
    end
    return UpperTriangular(A), UpperTriangular_pullback
end

function rrule(::Type{<:LowerTriangular}, A::AbstractMatrix)
    function LowerTriangular_pullback(ȳ)
        return (NoTangent(), Matrix(ȳ))
    end
    return LowerTriangular(A), LowerTriangular_pullback
end

function rrule(::typeof(triu), A::AbstractMatrix, k::Integer)
    function triu_pullback(ȳ)
        return (NoTangent(), triu(ȳ, k), NoTangent())
    end
    return triu(A, k), triu_pullback
end
function rrule(::typeof(triu), A::AbstractMatrix)
    function triu_pullback(ȳ)
        return (NoTangent(), triu(ȳ))
    end
    return triu(A), triu_pullback
end

function rrule(::typeof(tril), A::AbstractMatrix, k::Integer)
    function tril_pullback(ȳ)
        return (NoTangent(), tril(ȳ, k), NoTangent())
    end
    return tril(A, k), tril_pullback
end
function rrule(::typeof(tril), A::AbstractMatrix)
    function tril_pullback(ȳ)
        return (NoTangent(), tril(ȳ))
    end
    return tril(A), tril_pullback
end

_diag_view(X) = view(X, diagind(X))
_diag_view(X::Diagonal) = parent(X)  #Diagonal wraps a Vector of just Diagonal elements

function rrule(::typeof(det), X::Union{Diagonal,AbstractTriangular})
    y = det(X)
    s = conj!(y ./ _diag_view(X))
    function det_pullback(ȳ)
        return (NoTangent(), Diagonal(ȳ .* s))
    end
    return y, det_pullback
end

function rrule(::typeof(logdet), X::Union{Diagonal,AbstractTriangular})
    y = logdet(X)
    s = conj!(one(eltype(X)) ./ _diag_view(X))
    function logdet_pullback(ȳ)
        return (NoTangent(), Diagonal(ȳ .* s))
    end
    return y, logdet_pullback
end


text:
using FFTW
using PyPlot
using DSP
using WAV

function downsample(data, M)
       order = 35 # filter order
       d = 0.01 # transition band
       coeffs = remez(order, [0, 1/(M*2) - d, 1/(M*2)+d, 0.5], [1, 0])
       return (decimate(filt(coeffs, [1], data), M), coeffs, [1])
end

function decimate(data, M)
       return data[1:M:end]
end

function discriminator(data, M)
       # Limiter
       data = data ./ abs.(data)

       # differentiator
       (data2, b2, a2) = differentiator(data, M)

       # complex conj delay
       (ds, b, a) = delay((conj.(data2)), (M-1)/2.0, 80)
       data_mod = ds.* data2
       return (imag.(data_mod), b, a, b2, a2)
end

function differentiator(data, M)
       b = remez(M, [0, 0.5], [1.0], filter_type=RemezFilterType(2))
       return (filt(b, [1], data), b, [1])
end

function delay(data, D, o)
       a = 0.1
       b = [sin.(a*(n - D))/(a*(n - D) ) + 0im for n = 0:o-1]
       c = filt(b, a, data)
       return (c, b, [1])
end

function deemphasis(data, t, f)
       T = 1/f
       al = 1/tan(T/(2*t))
       b = [1, 1]
       a = [1+al, 1-al]
       return (filt(b, a, data), b, a)
end

function lowPassFilter(data, f, pb, sb)
       fp = pb/f
       ft = sb/f
       order = 35
       coeffs = remez(order, [0, fp, ft, 0.5], [1, 0])
       return (filt(coeffs, [1], data), coeffs, [1])
end

function plotTimeFreq(storeFft, fs, f0)
    w, h = figaspect(0.33)
    figure(figsize=[w,h])
    minF = (f0 - fs/2 )./ 1e6
    maxF = (f0 + fs/2 )./ 1e6
    maxT = 1/fs * size(storeFft)[1] * size(storeFft)[2]
    imshow(storeFft, extent=[minF,maxF,0,maxT], aspect="auto")
    xlabel("Frequnecy (MHz)")
    ylabel("Time (s)")
    cb = colorbar()
    cb[:set_label]("Magnitude (dB)")
end

function plotTime(data, fs)
    w, h = figaspect(0.33)
    figure(figsize=[w,h])
    maxT = 1/fs * size(data)[1]
    t = range(0, maxT, length=size(data)[1])
    plot(t, data, linewidth=0.1)
    xlim([0, maxT])
    xlabel("Time (s)")
    ylabel("Amplitude")
end

function plotIq(data)
    figure()
    scatter(real.(data), imag.(data), s=10)
    xlabel("In-phase")
    ylabel("Quadrature")
    title("IQ data")
end

function fmDemod(data, fs)
    # downsample
    fs2 = 256e3
    M = Int(fs/fs2)
    (data3, b3, a3) = downsample(data, M);

    (data4, b4, a4, b42, a42) = discriminator(data3, 10);

    t = 75e-6
    (data5, b5, a5) = deemphasis(data4, t, fs2);

    (data6, b6, a6) = lowPassFilter(data5, fs2, 15e3, 18e3)

    fs3 = 64e3
    M = Int(fs2 / fs3)
    (data7, b7, a7) = downsample(data6, M)

    maxval = abs(data7[argmax(abs.(data7))])
    data8 = data7./maxval
    return (data8, fs3)
    return (data7, fs3)
end


text:
module ScSTO

using OptiMo

export scstoproblem

include("types.jl")
include("interface.jl")
include("l0norm.jl")
include("evaluator.jl")
include("simulation.jl")

end


text:
# Generate structures that fill distance constraints

export
    fillconstraints,
    errorscore,
    fixchirality!,
    addstructure!,
    removemodulator!,
    generateensemble,
    trimbyscore!

using ProgressMeter

"""
Use stochastic proximity embedding to generate a structure iteratively.
See Agrafiotis et.al., Chapter 14 in Mucherino et. al., 2013 for algorithm details (Algorithm 8).
Returns coordinate list.
"""
function fillconstraints(constraints::Constraints;
                    iters_per_atom::Integer=defaults["iters_per_atom"],
                    cyc_step_ratio::Real=defaults["cyc_step_ratio"],
                    box_size::Real=defaults["box_size"]
    )
    @assert iters_per_atom > 0 "iters_per_atom must be positive"
    @assert cyc_step_ratio > 0 "cyc_step_ratio must be positive"
    constraints_lower = constraints.lower
    constraints_upper = constraints.upper
    present_i = constraints.pres_inds[:,1]
    present_j = constraints.pres_inds[:,2]
    squares_lower = constraints_lower .* constraints_lower
    squares_upper = constraints_upper .* constraints_upper
    n_atoms = length(constraints.atoms)
    # Work out n_cycles and n_steps using n_atoms, iters_per_atom and the ratio of steps to cycles
    n_cycs_steps = n_atoms * iters_per_atom
    n_cycles = round(Int, sqrt(n_cycs_steps / cyc_step_ratio))
    n_steps = round(Int, sqrt(n_cycs_steps * cyc_step_ratio))
    # Randomise atomic coordinates in a cube
    coords = rand!(zeros(3, n_atoms)) * box_size
    n_constraints = length(present_i)
    learning_rate = 1.0
    coord_new_i = zeros(3)
    coord_new_j = zeros(3)
    constraint_near = 0.0
    rand_list = zeros(Int, n_steps)
    @inbounds for cyc in 1:n_cycles
        # List of random index tuples
        rand!(rand_list, 1:n_constraints)
        for rand_ind in rand_list
            # Pick a random constraint
            i = present_i[rand_ind]
            j = present_j[rand_ind]
            # Compute distance between atoms
            sq_dist = 0.0
            for m = 1:3
                sq_dist += abs2(coords[m, i] - coords[m, j])
            end
            # Check if distance is within constraints
            if sq_dist < squares_lower[rand_ind]
                constraint_near = constraints_lower[rand_ind]
            elseif sq_dist > squares_upper[rand_ind]
                constraint_near = constraints_upper[rand_ind]
            else
                continue
            end
            # Adjust coordinates
            dist = sqrt(sq_dist)
            prefactor = learning_rate * 0.5 * (constraint_near - dist) / dist
            for k in 1:3
                coord_new_i[k] = newcoords(k, i, j, prefactor, coords)
                coord_new_j[k] = newcoords(k, j, i, prefactor, coords)
            end
            for l in 1:3
                coords[l, i] = coord_new_i[l]
                coords[l, j] = coord_new_j[l]
            end
        end
        # Decrease learning rate
        learning_rate -= 1 / n_cycles
    end
    return coords
end

# Speeds up fillconstraints
newcoords(k::Integer, i::Integer, j::Integer, prefactor::Float64, coords::Array{Float64,2}) = coords[k, i] + prefactor * (coords[k, i] - coords[k, j])

"Calculate SPE error score."
function errorscore(constraints::Constraints, coords::Array{Float64})
    constraints_lower = constraints.lower
    constraints_upper = constraints.upper
    pres_inds = constraints.pres_inds
    score = 0.0
    for m in 1:size(pres_inds, 1)
        k, l = pres_inds[m,:]
        dist = norm(coords[:,k] - coords[:,l])
        if dist < constraints_lower[m]
            # Arbitrary minimum constraint size used here to prevent division by zero and infinite scores
            score += ((dist-constraints_lower[m])^2) / max((constraints_upper[m]-constraints_lower[m]), 0.001)
        elseif dist > constraints_upper[m]
            score += ((dist-constraints_upper[m])^2) / max((constraints_upper[m]-constraints_lower[m]), 0.001)
        end
    end
    return score
end

"""
Checks the chirality of coordinate set and switches it if it isn't correct.
The coordinates must correspond to the atoms.
Inverts each z coordinate in the coordinate array if required.
"""
function fixchirality!(coords::Array{Float64}, atoms::Array{Atom,1})
    n_atoms = length(atoms)
    @assert n_atoms > 0 "No atoms in atom list"
    @assert n_atoms == size(coords, 2) "Number of atoms and coordinates are different"
    count_l = 0
    count_d = 0
    for outer_id in 1:n_atoms
        if atoms[outer_id].atom_name == "CA"
            res_n = atoms[outer_id].res_n
            chain_id = atoms[outer_id].chain_id
            coord_ca = coords[:, outer_id]
            found_co = false
            found_n = false
            found_cb = false
            vec_oc = zeros(3)
            vec_nc = zeros(3)
            vec_rc = zeros(3)
            for inner_id in 1:n_atoms
                if atoms[inner_id].res_n == res_n && atoms[inner_id].chain_id == chain_id
                    if atoms[inner_id].atom_name == "C"
                        found_co = true
                        vec_oc = coord_ca - coords[:,inner_id]
                    elseif atoms[inner_id].atom_name == "N"
                        found_n = true
                        vec_nc = coord_ca - coords[:,inner_id]
                    elseif atoms[inner_id].atom_name == "CB"
                        found_cb = true
                        vec_rc = coord_ca - coords[:,inner_id]
                    end
                end
            end
            if found_co && found_n && found_cb
                if dot(vec_oc, LinearAlgebra.cross(vec_rc, vec_nc)) > 0
                    count_d += 1
                else
                    count_l += 1
                end
            end
        end
    end
    @assert max(count_d, count_l) > 0 "Could not find any tetrahedral C-alpha systems to determine chirality"
    if count_d > count_l
        coords[3,:] *= -1
    end
end

"Add a structure to an ensemble."
function addstructure!(ensemble::ModelledEnsemble, struc::ModelledStructure)
    n_coords = size(struc.coords, 2)
    @assert n_coords > 0 "The structure to add does not contain any coordinates"
    strucs = ensemble.strucs
    if length(strucs) > 0
        @assert n_coords == size(strucs[1].coords, 2) "Attempting to add a structure with a different number of atoms to the existing structures"
    end
    push!(strucs, struc)
end

"Takes a `ModelledEnsemble` and removes the fake modulator atoms and coords."
function removemodulator!(ensemble::ModelledEnsemble)
    atoms = ensemble.atoms
    n_atoms = length(atoms)
    @assert n_atoms > 0 "There are no atoms in the atom list"
    strucs = ensemble.strucs
    n_strucs = length(strucs)
    @assert n_strucs == 0 || size(strucs[1].coords, 2) == n_atoms "Number of coordinates and number of atoms are different"
    # Find indices corresponding to non-modulator atoms
    atom_inds = Int[]
    for i in 1:n_atoms
        if !ismodulator(atoms[i])
            push!(atom_inds, i)
        end
    end
    if length(atom_inds) != n_atoms
        # Remove modulator atoms from atom list
        ensemble.atoms = getindex(atoms, atom_inds)
        # Remove modulator coords from each structure
        for j in 1:n_strucs
            ensemble.strucs[j].coords = getindex(strucs[j].coords, collect(1:3), atom_inds)
        end
    end
end

"""
Generate sets of coordinates with the correct chirality from distance constraints.
discard_ratio gives the fraction of extra structures to calculate, then the low-scoring structures are discarded.
If `remove_mod` is true then modulator atoms and coords are removed from the ensemble before it is returned.
Returns a `ModelledEnsemble`.
"""
function generateensemble(constraints::Constraints,
                    n_strucs::Integer;
                    iters_per_atom::Integer=defaults["iters_per_atom"],
                    discard_ratio::Real=defaults["discard_ratio"],
                    remove_mod::Bool=true)
    atoms = constraints.atoms
    n_atoms = length(atoms)
    @assert n_atoms > 0 "The Constraints object does not have the atoms set"
    n_constraints = size(constraints.pres_inds, 1)
    @assert n_constraints == length(constraints.lower) && n_constraints == length(constraints.upper) "The number of lower and upper constraints do not correspond to the number of constraints present"
    @assert n_constraints > 0 "There are no constraints present"
    @assert discard_ratio >= 1.0 "discard_ratio must be >= 1"
    ensemble = ModelledEnsemble(atoms)
    n_strucs_init = round(Int, ceil(n_strucs * discard_ratio))
    println("Will generate ", n_strucs_init, " structures and keep best-scoring ", n_strucs)
    p = Progress(n_strucs_init, 1, "Progress: ", 50)
    for struc_n in 1:n_strucs_init
        coords = fillconstraints(constraints, iters_per_atom=iters_per_atom)
        score = errorscore(constraints, coords)
        fixchirality!(coords, atoms)
        addstructure!(ensemble, ModelledStructure(score, coords))
        next!(p)
    end
    trimbyscore!(ensemble, n_strucs)
    if remove_mod
        removemodulator!(ensemble)
    end
    return ensemble
end

"Retains only the n_to_keep top scoring coordinate sets in an ensemble."
function trimbyscore!(ensemble::ModelledEnsemble, n_to_keep::Integer)
    strucs = ensemble.strucs
    @assert n_to_keep > 0 "n_to_keep must be positive"
    @assert n_to_keep <= length(strucs) "Asking to trim an ensemble of size $(length(strucs)) to $n_to_keep structures"
    scores = [struc.score for struc in strucs]
    inds_to_keep = sortperm(scores)[1:n_to_keep]
    ensemble.strucs = getindex(strucs, inds_to_keep)
    println("Kept lowest scoring ", n_to_keep, " out of ", length(scores), " structures")
    println("Scores range from ", round(sort(scores)[1], digits=1), " to ", round(sort(scores)[n_to_keep], digits=1))
end


text:
using Downloads
using GadgetReader
using GadgetIO
using LinearAlgebra
using Test
using Unitful
using UnitfulAstro

D = Downloads
GG = GadgetReader

@testset "GadgetReader.jl" begin

    @info "downloading test data..."
    url = "http://www.usm.uni-muenchen.de/~lboess/GadgetIO/"
    for i in 0:3
        D.download(joinpath(url, "sub_002.$i"), "./sub_002.$i")
        D.download(joinpath(url, "snap_002.$i"), "./snap_002.$i")
        D.download(joinpath(url, "snap_002.$i.key"), "./snap_002.$i.key")
    end
    D.download(joinpath(url, "snap_002.key.index"), "./snap_002.key.index")
    @info "done!"

    @testset "Snapshot" begin
        @test GG.get_snapbase("box", 32) == joinpath("box", "snapdir_032", "snap_032")
        @test GG.get_subbase("box", 2) == joinpath("box", "groups_002", "sub_002")

        snapshot = Snapshot("snapbase", "subbase")
        @test snapshot.snapbase == "snapbase"
        @test snapshot.subbase == "subbase"

        snapshot = Snapshot("snapbase", nothing)
        @test snapshot.snapbase == "snapbase"
        @test isnothing(snapshot.subbase)

        snapshot = Snapshot(nothing, "subbase")
        @test isnothing(snapshot.snapbase)
        @test snapshot.subbase == "subbase"

        @test Snapshot(snapbase="snapbase") == Snapshot("snapbase", nothing)
        @test Snapshot(subbase="subbase") == Snapshot(nothing, "subbase")

        snapbase = joinpath("box", "snapdir_010", "snap_010")
        subbase = joinpath("box", "groups_010", "sub_010")

        snapdir = joinpath("snapdir_010", "snap_010")
        subdir = joinpath("groups_010", "sub_010")
        @test Snapshot("box", snapdir, subdir) == Snapshot(snapbase, subbase)
        @test Snapshot("box", 10) == Snapshot(snapbase, subbase)
        @test Snapshot("box", 10; snapbase=false) == Snapshot(nothing, subbase)
        @test Snapshot("box", 10; subbase=false) == Snapshot(snapbase, nothing)

        snapshot = Snapshot(snapbase, subbase)
        io = IOBuffer()
        show(io, "text/plain", snapshot)
        @test String(take!(io)) == "Snapshot\nSnap: $snapbase\nSub:  $subbase"
    end

    @testset "Galaxy" begin
        # particles
        p = GG.Particles(:stars, Dict("ID" => [1, 2], "MASS" => [1e11, 2e11]))
        p.zs = [1 2; 3 4; 5 6]
        p.iM = [3, 4]

        @test p.type === :stars
        @test p.id == [1, 2]
        @test p.mass == [1e11, 2e11]
        @test p.zs == [1 2; 3 4; 5 6]
        @test p.iM == [3, 4]
        @test p.zs == p[:zs]
        @test p.properties["Zs"] === p.zs
        @test p.properties["MASS"] === p.mass
        @test p.properties["Zs"] === p["Zs"]
        @test p.properties["iM"] === p.iM
        @test issetequal(keys(p), [:id, :mass, :zs, :im])
        @test haskey(p, :id)
        @test haskey(p, "ID")
        @test !haskey(p, :pos)
        @test !haskey(p, "POS")
        @test values(p) === values(p.properties)
        @test issetequal(propertynames(p), [:type, :id, :mass, :zs, :im])

        p[:zs] = [6 5; 4 3; 2 1]
        @test p.Zs == [6 5; 4 3; 2 1]
        @test p.zs === p[:ZS]
        p["Zs"] = [1 2; 3 4; 5 6]
        @test p.Zs == [1 2; 3 4; 5 6]
        @test p.zs === p[:ZS]

        @test_throws ErrorException p[:type] = :dm

        @test particle_type_id(:gas) == 0

        q = copy(p)
        @test q !== p
        @test q.id === p.id
        q.id = p.id .+ 1
        @test q.id !== p.id
        @test q.id[1] != p.id[1]

        # galaxy
        snapshot = Snapshot("box", 10)
        dm = GG.Particles(:dm, Dict())
        g = Galaxy(snapshot, 1234, nothing, Dict(:stars => p))
        g.dm = dm

        @test g.snapshot === snapshot
        @test GG.getind(g) === g.isub === 1234
        @test isnothing(g.subid) && isnothing(GG.getid(g))
        @test g.stars === g[:stars] === p
        @test g.dm === g[:dm] === dm
        @test g.particles[:stars] === g.stars
        @test g.particles[:dm] === g.dm
        @test issetequal(keys(g), [:stars, :dm])
        @test haskey(g, :stars)
        @test !haskey(g, :gas)
        @test values(g) === values(g.particles)
        @test issetequal(propertynames(g), [:snapshot, :isub, :subid, :stars, :dm])

        g[:bh] = dm
        @test g.bh === dm

        h = copy(g)
        @test h !== g
        @test h.stars !== g.stars
        @test h.stars.id === g.stars.id
        h.stars.id = g.stars.id .+ 1
        @test h.stars.id !== g.stars.id

        g = Galaxy(snapshot, 1234, false)
        @test g.snapshot === snapshot
        @test g.isub === 1234
        @test isnothing(g.subid)
        @test isempty(keys(g))
        @test g.isub === 1234

        g.stars = p
        @test g[:stars] === p

        @test_throws ErrorException g[:isub] = 1235

        io = IOBuffer()
        show(io, "text/plain", g)
        @test String(take!(io)) == "Galaxy isub 1234\nstars: 2 Particles\n ID MASS Zs iM\n"

        snapshot = Snapshot("snap_002", "sub_002")
        g = Galaxy(snapshot, 0)
        @test g.subid == HaloID(0, 1)

        # galaxy group
        snapshot = Snapshot("box", 10)

        gr = GalaxyGroup(snapshot, 1234, nothing, Dict(:stars => p))
        gr.dm = dm

        @test gr.snapshot === snapshot
        @test GG.getind(gr) === gr.igroup === 1234
        @test isnothing(gr.groupid) && isnothing(GG.getid(gr))
        @test gr.stars === gr[:stars] === p
        @test gr.dm === gr[:dm] === dm
        @test gr.particles[:stars] === gr.stars
        @test gr.particles[:dm] === gr.dm
        @test issetequal(keys(gr), [:stars, :dm])
        @test values(gr) === values(gr.particles)
        @test issetequal(propertynames(gr), [:snapshot, :igroup, :groupid, :stars, :dm])

        gr[:bh] = dm
        @test gr.bh === dm

        gr = GalaxyGroup(snapshot, 1234, false)
        @test gr.snapshot === snapshot
        @test gr.igroup === 1234
        @test isnothing(gr.groupid)
        @test isempty(keys(gr))
        @test gr.igroup === 1234

        gr.stars = p
        @test gr[:stars] === p

        @test_throws ErrorException gr[:igroup] = 1235

        io = IOBuffer()
        show(io, "text/plain", gr)
        @test String(take!(io)) == "GalaxyGroup igroup 1234\nstars: 2 Particles\n ID MASS Zs iM\n"

        snapshot = Snapshot("snap_002", "sub_002")
        gr = GalaxyGroup(snapshot, 0)
        @test gr.groupid == HaloID(0, 1)
    end

    # create real snapshot for rest of tests
    snapshot = Snapshot("snap_002", "sub_002")
    g = Galaxy(snapshot, 0)
    gr = GalaxyGroup(snapshot, 0)
    h = read_header(snapshot.subbase)

    @testset "Units" begin
        v, w = 1.0, 1.0f0
        vs = ones(Float64, 3)
        ws = ones(Float32, 3)

        a = 0.277777
        b = a * u"kpc"
        @test convert_units_pos(v, h, :sim) == v
        @test convert_units_pos(v, h, :physical) ≈ a rtol = 1e-5
        @test convert_units_physical_pos(v, h) ≈ a rtol = 1e-5
        @test convert_units_physical_pos(v, h) |> typeof === Float64
        @test convert_units_physical_pos(w, h) ≈ a rtol = 1e-5
        @test convert_units_physical_pos(w, h) |> typeof === Float32
        @test convert_units_physical_pos(vs, h)[1] ≈ a rtol = 1e-5
        @test convert_units_physical_pos(vs, h) |> eltype === Float64
        @test convert_units_physical_pos(ws, h)[1] .≈ a rtol = 1e-5
        @test convert_units_physical_pos(ws, h) |> eltype === Float32
        @test convert_units_pos(v, h, :full) ≈ b rtol = 1e-5
        @test convert_units_full_pos(v, h) ≈ b rtol = 1e-5
        @test convert_units_full_pos(v, h) |> typeof <: Quantity{Float64}
        @test convert_units_full_pos(w, h) ≈ b rtol = 1e-5
        @test convert_units_full_pos(w, h) |> typeof <: Quantity{Float32}
        @test convert_units_full_pos(vs, h)[1] ≈ b rtol = 1e-5
        @test convert_units_full_pos(vs, h) |> eltype <: Quantity{Float64}
        @test convert_units_full_pos(ws, h)[1] .≈ b rtol = 1e-5
        @test convert_units_full_pos(ws, h) |> eltype <: Quantity{Float32}

        af = convert(Float32, a)
        bf = af * u"kpc"
        as = [a, a, a]
        asf = [af, af, af]
        bs = [b, b, b]
        bsf = [bf, bf, bf]
        @test simulation_units_pos(a, h) |> typeof === Float64
        @test simulation_units_pos(af, h) |> typeof === Float32
        @test simulation_units_pos(b, h) |> typeof === Float64
        @test simulation_units_pos(bf, h) |> typeof === Float32
        @test simulation_units_pos(a, h) ≈ v rtol = 1e-5
        @test simulation_units_pos(b, h) ≈ v rtol = 1e-5
        @test simulation_units_pos(b |> u"m", h) ≈ v rtol = 1e-5
        @test simulation_units_pos(as, h) |> eltype === Float64
        @test simulation_units_pos(asf, h) |> eltype === Float32
        @test simulation_units_pos(bs, h) |> eltype === Float64
        @test simulation_units_pos(bsf, h) |> eltype === Float32
        @test simulation_units_pos(as, h)[1] ≈ v rtol = 1e-5
        @test simulation_units_pos(asf, h)[1] ≈ v rtol = 1e-5
        @test simulation_units_pos(bs, h)[1] ≈ v rtol = 1e-5
        @test simulation_units_pos(bsf, h)[1] ≈ v rtol = 1e-5
        @test simulation_units_pos(bs .|> u"m", h)[1] ≈ v rtol = 1e-5
        simulation_units_pos!(as, h)
        @test as[1] ≈ v rtol = 1e-5

        vs2, ws2 = copy(vs), copy(ws)
        convert_units_physical_pos!(vs2, h)
        @test vs2[1] ≈ a rtol = 1e-5
        convert_units_physical_pos!(ws2, h)
        @test ws2[1] ≈ a rtol = 1e-5

        # TODO: perhaps also test vel, temp, mass

        @test convert_units_physical(vs, :pos, h) == convert_units_physical_pos(vs, h)
        @test convert_units_physical(vs, :vel, h) == convert_units_physical_vel(vs, h)
        @test convert_units_physical(vs, :temp, h) == convert_units_physical_temp(vs, h)
        @test convert_units_physical(vs, :zs, h) == convert_units_physical_mass(vs, h)

        v, w = 0.1, 0.1
        vs = ones(Float64, 3) / 10
        ws = ones(Float32, 3) / 10

        a = 1.06078
        b = a * Unitful.Gyr
        @test convert_units_age(v, h, :sim) == v
        @test convert_units_age(v, h, :physical) ≈ a rtol = 1e-5
        @test convert_units_physical_age(v, h) ≈ a rtol = 1e-5
        @test convert_units_physical_age(v, h) |> typeof === Float64
        @test convert_units_physical_age(w, h) ≈ a rtol = 1e-5
        @test convert_units_physical_age(w, h) |> typeof === Float64
        @test convert_units_physical_age(vs, h)[1] ≈ a rtol = 1e-5
        @test convert_units_physical_age(vs, h) |> eltype === Float64
        @test convert_units_physical_age(ws, h)[1] ≈ a rtol = 1e-5
        @test convert_units_physical_age(ws, h) |> eltype === Float64
        @test convert_units_age(v, h, :full) ≈ b rtol = 1e-5
        @test convert_units_full_age(v, h) ≈ b rtol = 1e-5
        @test convert_units_full_age(v, h) |> typeof <: Quantity{Float64}
        @test convert_units_full_age(w, h) ≈ b rtol = 1e-5
        @test convert_units_full_age(w, h) |> typeof <: Quantity{Float64}
        @test convert_units_full_age(vs, h)[1] ≈ b rtol = 1e-5
        @test convert_units_full_age(vs, h) |> eltype <: Quantity{Float64}
        @test convert_units_full_age(ws, h)[1] ≈ b rtol = 1e-5
        @test convert_units_full_age(ws, h) |> eltype <: Quantity{Float64}
        @test convert_units_physical(vs, :age, h) == convert_units_physical_age(vs, h)
        @test convert_units_physical!(copy(vs), :age, h) == convert_units_physical_age(vs, h)
        @test convert_units_full(vs, :age, h) == convert_units_full_age(vs, h)

        vs2, ws2 = copy(vs), copy(ws)
        convert_units_physical_age!(vs2, h)
        @test vs2[1] ≈ a rtol = 1e-5
        convert_units_physical_age!(ws2, h)
        @test ws2[1] ≈ a rtol = 1e-5

        zs, mass = ones(11, 2), [42, 34]
        a = 24.0731
        @test convert_units_solar_metallicity(zs, mass)[1] ≈ a rtol = 1e-5

        p = GG.Particles(:stars, Dict())
        p.pos, p.vel, p.zs = rand(3, 2), rand(3, 2), rand(11, 2)
        p.iM, p.mass, p.age, p.temp = rand(2), rand(2), rand(2), rand(2)
        convert_units!(p, h, :physical)
        @test eltype(p.pos) === Float64
        @test eltype(p.age) === Float64
        @test eltype(p.zs) === Float64
        @test size(p.zs) == (2,)

        p.pos, p.vel, p.zs = rand(3, 2), rand(3, 2), rand(11, 2)
        p.iM, p.mass, p.age, p.temp = rand(2), rand(2), rand(2), rand(2)
        convert_units!(p, h, :full)
        @test eltype(p.pos) <: Unitful.Length
        @test eltype(p.vel) <: Unitful.Velocity
        @test eltype(p.age) <: Unitful.Time
        @test eltype(p.temp) <: Unitful.Temperature
        @test eltype(p.mass) <: Unitful.Mass
        @test eltype(p.zs) === Float64
        @test size(p.zs) == (2,)
    end

    @testset "Utils" begin
        boxsize, boxsize_half = 10, 5
        @test GG.shift_across_box_border(1, 2, boxsize, boxsize_half) == 1
        @test GG.shift_across_box_border(8, 2, boxsize, boxsize_half) == -2
        @test GG.shift_across_box_border(2, 8, boxsize, boxsize_half) == 12
    end

    @testset "Read" begin
        @test read_redshift(snapshot) ≈ 4 rtol = 1e-5
        @test read_header_particle_mass(snapshot, :dm, :physical) == 0
        @test read_header_particle_mass(snapshot, :dm) == 0u"Msun"
        @test read_header_particle_mass(snapshot, :gas) == 0u"Msun"
        @test read_galaxy_prop(g, "MSUB", :sim) ≈ 9.17776 rtol = 1e-5
        @test read_galaxy_prop(g, "MSUB", :physical) ≈ 1.27468f11 rtol = 1e-5
        @test read_galaxy_prop(g, "MSUB", :full) ≈ 1.27468f11u"Msun" rtol = 1e-5
        @test read_galaxy_prop(g, "SPOS") |> length == 3
        @test read_galaxy_prop(g, "SPOS")[1] ≈ -413.854u"kpc" rtol = 1e-5
        # @test read_galaxy_prop(g, "SVEL")[1] ≈ 83.1962u"km/s" rtol = 1e-5
        @test read_galaxy_prop(g, "SVEL")[1] ≈ 186.032u"km/s" rtol = 1e-5
        @test read_galaxy_prop(g, "GRNR") == 0
        info = "Reading property GRNR of halo HaloID(0, 1)"
        warn = "The quantity GRNR is returned in simulation units"
        @test_logs (:info, info) (:warn, warn) read_galaxy_prop(g, "GRNR"; verbose=true)
        @test GG.convert_units_subfind_prop(0.1, "SAGE", h) ≈ 1.06078Unitful.Gyr rtol = 1e-5
        @test GG.convert_units_subfind_prop(0.1, "SSFR", h) == 0.1u"Msun" / Unitful.yr

        @test read_galaxy_pos(g) == read_galaxy_pos(gr)
        @test read_galaxy_vel(g) == read_galaxy_vel(gr)
        @test is_main_halo(g)

        @test get_group(g).groupid == gr.groupid
        @test get_first_subhalo(gr).subid == g.subid
    end

    @testset "Read Halo" begin
        props = ((:dm, ["POS"]), (:gas, ["MASS", "POS", "VEL"]))
        read_halo!(g; use_keys=false, props, units=:full)

        @test typeof(g) <: Galaxy
        @test typeof(g.dm) <: GG.Particles
        @test typeof(g.gas) <: GG.Particles

        io = IOBuffer()
        show(io, "text/plain", g)
        @test String(take!(io)) ==
              "Galaxy at z=4.0\n isub 0\ndm: 7 Particles\n ID POS\ngas: 20 Particles\n ID MASS POS VEL\n"

        read_halo!(gr; use_keys=false, props, units=:full)

        @test typeof(gr) <: GalaxyGroup
        @test typeof(gr.dm) <: GG.Particles
        @test typeof(gr.gas) <: GG.Particles

        io = IOBuffer()
        show(io, "text/plain", gr)
        @test String(take!(io)) ==
              "GalaxyGroup at z=4.0\n igroup 0\ndm: 70 Particles\n ID POS\ngas: 76 Particles\n ID MASS POS VEL\n"
    end

    @testset "Transformation" begin
        a = [1 2 3 4; 5 6 7 8; 9 0 1 1]
        b = [9 0 1 1; 5 6 7 8; -1 -2 -3 -4]
        R = [0 0 1; 0 1 0; -1 0 0]
        @test rotate(a, R) == b
        rotate!(a, R)
        @test a == b

        p11 = 2.45564u"kpc"
        # v12 = 90.39079u"km/s"
        v12 = -12.44545u"km/s"
        @test g.dm.pos[1, 1] ≈ p11 rtol = 1e-5
        @test g.gas.vel[1, 2] ≈ v12 rtol = 1e-5
        gc = rotate(g, R)
        dmc = rotate(g.dm, R)
        @test gc.dm.pos[3, 1] ≈ -p11 rtol = 1e-5
        @test dmc.pos[3, 1] ≈ -p11 rtol = 1e-5
        @test gc.gas.vel[3, 2] ≈ -v12 rtol = 1e-5
        @test g.dm.pos != gc.dm.pos
        @test g.gas.vel != gc.gas.vel
        rotate!(g, R)
        @test g.dm.pos[3, 1] ≈ -p11 rtol = 1e-5
        @test g.gas.vel[3, 2] ≈ -v12 rtol = 1e-5
        rotate!(dmc, R)
        @test dmc.pos[3, 1] ≈ 3.05406u"kpc" rtol = 1e-5

        @test rotate_edgeon(gr.gas; algorithm=:unw_i_vol).pos != gr.gas.pos
        @test rotate_edgeon(gr, :gas; radius=20u"kpc").dm.pos != gr.dm.pos

        gr2 = deepcopy(gr)
        rotate_edgeon!(gr2.gas).pos
        @test gr2 != gr.gas.pos
        gr2 = deepcopy(gr)
        rotate_edgeon!(gr2, :gas).dm.pos
        @test gr2.dm.pos != gr.dm.pos

        @test rotate_edgeon(gr.gas; axis_ratios=true) isa Tuple
        @test rotate_edgeon!(gr2.gas; axis_ratios=true) isa Tuple
        @test rotate_edgeon(gr, :gas; axis_ratios=true) isa Tuple
        @test rotate_edgeon!(gr2, :gas; axis_ratios=true) isa Tuple
    end

    @testset "Shapes" begin
        let err = nothing
            try
                GG.get_algorithm_variables(:abc)
            catch err
            end
            @test err isa ErrorException
        end

        # 2d to 3d rotation matrix
        θ = π / 3
        R2d = [cos(θ) -sin(θ); sin(θ) cos(θ)]
        R3dx = GG.rotation_matrix_to_3d(R2d, [2, 3])
        R3dy = GG.rotation_matrix_to_3d(R2d, [3, 1])
        R3dz = GG.rotation_matrix_to_3d(R2d, [1, 2])
        @test R3dx == [1 0 0; 0 cos(θ) -sin(θ); 0 sin(θ) cos(θ)]
        @test R3dy == [cos(θ) 0 sin(θ); 0 1 0; -sin(θ) 0 cos(θ)]
        @test R3dz == [cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1]

        # parameters
        @test eccentricity(0.3) ≈ sqrt(1 - 0.3^2)
        @test ellipticity(0.3) ≈ 1 - 0.3
        @test triaxiality(0.7, 0.6) ≈ (1 - 0.7^2) / (1 - 0.6^2)

        # rotation matrix
        𝐈 = [2 -1 0; -1 2 -1; 0 -1 2]
        Q⁻¹, q, s = rotation_matrix_axis_ratios(𝐈)
        @test Q⁻¹[1, 1] ≈ Q⁻¹[1, 3] ≈ Q⁻¹[3, 1] ≈ Q⁻¹[3, 3] ≈ 0.5
        @test -Q⁻¹[1, 2] ≈ Q⁻¹[2, 1] ≈ -Q⁻¹[2, 3] ≈ Q⁻¹[3, 2] ≈ √0.5
        @test Q⁻¹[2, 2] ≈ 0 atol = 1e-14
        @test q ≈ 0.765366 rtol = 1e-5
        @test s ≈ 0.414213 rtol = 1e-5

        Q⁻¹, q, s, axes = rotation_matrix_axis_ratios(𝐈; return_axes=true)
        @test axes[1][1] ≈ axes[1][3] ≈ axes[3][1] ≈ axes[3][3] ≈ 0.5
        @test axes[2][1] ≈ -axes[2][3] ≈ -axes[1][2] ≈ axes[3][2] ≈ √0.5
        @test axes[2][2] ≈ 0 atol = 1e-14

        # inertia matrix: upper triangle
        𝐈 = zeros(Float32, 3, 3)
        𝐈_unw = [3103.08 -240.712 24.434; 0.0 19302.6 3955.77; 0.0 0.0 5671.47]
        𝐈_red = [16.1944 -2.31039 0.309565; 0.0 39.9156 3.32378; 0.0 0.0 19.89]
        GG.fill_𝐈_up!(𝐈, gr.gas.pos, gr.gas.mass, false, :unweighted)
        @test isapprox.(𝐈, 𝐈_unw) |> all
        GG.fill_𝐈_up!(𝐈, gr.gas.pos, gr.gas.mass, false, :reduced)
        @test isapprox.(𝐈, 𝐈_red) |> all
        GG.fill_𝐈_up!(𝐈, gr.gas.pos, gr.gas.mass, true, :unweighted)
        @test isapprox.(𝐈 ./ ustrip(gr.gas.mass[1]), 𝐈_unw) |> all # only works because of same gas mass of all
        GG.fill_𝐈_up!(𝐈, gr.gas.pos, gr.gas.mass, true, :reduced)
        @test isapprox.(𝐈 ./ ustrip(gr.gas.mass[1]), 𝐈_red) |> all # only works because of same gas mass of all

        # distances
        @test r²_ellipsoid(gr.gas.pos, 0.7, 0.5)[1] ≈ 573.1933u"kpc^2" rtol = 1e-5
        @test r²_ellipsoid(gr.gas.pos, 0.7, nothing)[1] ≈ 514.1716u"kpc^2" rtol = 1e-5
        @test r²_sphere(gr.gas.pos)[1] ≈ r²_ellipsoid(gr.gas.pos, 1, 1)[1] rtol = 1e-5
        @test r²_circle([3 2 3; 4 3 1])[1] == 25
        @test r²_sphere(gr.gas.pos[[1, 2], :])[1] ≈ r²_circle(gr.gas.pos[[1, 2], :])[1] rtol = 1e-5

        # ellipsoidal masking
        sph = Sphere(20u"kpc")
        @test GG.ellipsoidal_mask(gr.gas.pos, sph) |> count == 47
        ell = Ellipsoid(20u"kpc", 0.7, 0.5, :volume)
        @test GG.ellipsoidal_mask(gr.gas.pos, ell) |> count == 45
        @test GG.ellipsoidal_mask(gr.gas.pos, ell; return_r²_ellipsoid=true)[2][1] ≈ 573.1933u"kpc^2" rtol =
            1e-5
        ell = Ellipsoid(20u"kpc", 0.7, 0.5, :axis)
        @test GG.ellipsoidal_mask(gr.gas.pos, ell) |> count == 35

        # inertia matrices
        𝐈_unw = Symmetric(𝐈_unw, :U)
        𝐈_red = Symmetric(𝐈_red, :U)
        @test isapprox.(GG.inertia_matrix(gr.gas.pos, nothing; mass_weighted=false), 𝐈_unw) |> all
        @test isapprox.(GG.inertia_matrix(gr.gas.pos, gr.gas.mass) ./ ustrip(gr.gas.mass[1]), 𝐈_unw) |> all
        @test isapprox.(
            GG.inertia_matrix(gr.gas.pos, nothing; mass_weighted=false, inertia_matrix_type=:reduced),
            𝐈_red,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix(gr.gas.pos, gr.gas.mass; inertia_matrix_type=:reduced) ./
            ustrip(gr.gas.mass[1]),
            𝐈_red,
        ) |> all

        @test isapprox.(
            GG.inertia_matrix_iterative(gr.gas.pos, nothing, :unweighted; mass_weighted=false),
            𝐈_unw,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(gr.gas.pos, gr.gas.mass, :unweighted) ./ ustrip(gr.gas.mass[1]),
            𝐈_unw,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(gr.gas.pos, nothing, :reduced; mass_weighted=false),
            𝐈_red,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(gr.gas.pos, gr.gas.mass, :reduced) ./ ustrip(gr.gas.mass[1]),
            𝐈_red,
        ) |> all

        # iterative inertia matrices
        radius = 20u"kpc"
        𝐈_unwi = [1749.14 -1014.96 453.757; -1014.96 2522.14 -459.417; 453.757 -459.417 1034.89]
        𝐈_redi = [13.8342 -2.85969 1.23256; -2.85969 18.0241 -1.74392; 1.23256 -1.74392 14.1417]
        𝐈_redelli = [7.99452 -4.20033 1.454; -4.20033 12.72 -1.5627; 1.454 -1.5627 6.82152]
        𝐈_unwi_ax = [285.81 -420.516 164.467; -420.516 667.925 -272.656; 164.467 -272.656 145.668]
        𝐈_redi_ax = [13.1609 -2.25873 1.22684; -2.25873 16.9281 -2.09812; 1.22684 -2.09812 13.911]
        @test isapprox.(
            GG.inertia_matrix_iterative(gr.gas.pos, nothing, :unweighted; radius, mass_weighted=false),
            𝐈_unwi,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(gr.gas.pos, gr.gas.mass, :unweighted; radius) ./
            ustrip(gr.gas.mass[1]),
            𝐈_unwi,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(gr.gas.pos, nothing, :reduced; radius, mass_weighted=false),
            𝐈_redi,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(
                gr.gas.pos,
                nothing,
                :reduced;
                radius,
                mass_weighted=false,
                ellipsoidal_distance=true,
            ),
            𝐈_redelli,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(
                gr.gas.pos,
                nothing,
                :unweighted;
                radius,
                mass_weighted=false,
                constant_volume=false,
            ),
            𝐈_unwi_ax,
        ) |> all
        @test isapprox.(
            GG.inertia_matrix_iterative(
                gr.gas.pos,
                nothing,
                :reduced;
                radius,
                mass_weighted=false,
                constant_volume=false,
            ),
            𝐈_redi_ax,
        ) |> all

        # full rotation matrix
        R_unw = [
            0.013092664 -0.96556526 -0.25983113
            0.05690142 -0.25871283 0.96427697
            0.99829394 0.027409703 -0.05155479
        ]
        R_unw_r = [
            0.5127014 -0.8163861 0.2658031
            0.7872874 0.57054496 0.23378807
            0.34251392 -0.08939998 -0.9352497
        ]
        R_unw_r_i = [
            0.5581614 -0.78767705 0.2608086
            0.67917794 0.61428696 0.40170747
            0.4766269 0.047082156 -0.8778439
        ]
        R, q, s = rotation_matrix_edgeon(gr.gas; algorithm=:unw)
        @test isapprox.(R, R_unw; rtol=1e-5) |> all
        R, q, s = rotation_matrix_edgeon(gr.gas; algorithm=:unw, radius)
        @test isapprox.(R, R_unw_r; rtol=1e-5) |> all
        R, q, s = rotation_matrix_edgeon(gr.gas; algorithm=:unw_i_vol, radius)
        @test isapprox.(R, R_unw_r_i; rtol=1e-5) |> all

        # 2D rotation matrix
        R_unw = [0.009512025 0.0 0.9999547; 0.0 1.0 0.0; 0.9999547 0.0 -0.009512025]
        R_unw_r = [1.0 0.0 0.0; 0.0 0.9836683 -0.17999081; 0.0 0.17999081 0.98366827]
        R_unw_r_i = [0.43801472 -0.89896786 0.0; 0.898968 0.43801472 0.0; 0.0 0.0 1.0]
        R, q = rotation_matrix_2D(gr.gas; algorithm=:unw, perspective=:edgeon)
        @test isapprox.(R, R_unw; rtol=1e-5) |> all
        R, q = rotation_matrix_2D(gr.gas; algorithm=:unw, radius, perspective=:sideon)
        @test isapprox.(R, R_unw_r; rtol=1e-5) |> all
        R, q = rotation_matrix_2D(gr.gas; algorithm=:unw_i_vol, radius, perspective=:faceon)
        @test isapprox.(R, R_unw_r_i; rtol=1e-5) |> all
        @test rotation_matrix_2D(gr.gas; algorithm=:unw, matrix_2d=true)[1] |> size == (2, 2)

        @test_nowarn rotation_matrix_edgeon(gr.gas; algorithm=:red, radius)
        @test_nowarn rotation_matrix_edgeon(gr.gas; algorithm=:unw_i_ax, radius)
        @test_nowarn rotation_matrix_edgeon(gr.gas; algorithm=:red_i_vol, radius)
        @test_nowarn rotation_matrix_edgeon(gr.gas; algorithm=:red_i_ax, radius)
        @test_nowarn rotation_matrix_edgeon(gr.gas; algorithm=:red_ell_i_vol, radius)
        @test_nowarn rotation_matrix_edgeon(gr.gas; algorithm=:red_ell_i_ax, radius)

        let err = nothing
            try
                GG.rotation_matrix_2D(gr.gas; algorithm=:unw, matrix_2d=false, perspective=:abc)
            catch err
            end
            @test err isa ErrorException
        end

        # translation
        Δx = [1, 2, 3] .* u"kpc"
        gr2 = deepcopy(gr)
        x = gr2.gas.pos[:, 1]
        @test translate(gr2.gas, Δx).pos[:, 1] == x .+ Δx
        translate!(gr2.gas, Δx)
        @test gr2.gas.pos[:, 1] == x .+ Δx

        x = gr2.gas.pos[:, 1]
        @test translate(gr2, Δx).gas.pos[:, 1] == x .+ Δx
        translate!(gr2, Δx)
        @test gr2.gas.pos[:, 1] == x .+ Δx

        # center of mass
        gr2 = deepcopy(gr)
        a = [-0.036875404, -0.22094725, 0.5704074] .* u"kpc"
        rstart = 40u"kpc"
        @test isapprox.(center_of_mass_iterative(gr2, rstart, :gas; limit_fraction=0.01), a; rtol=1e-5) |> all
        @test isapprox.(
            translate_to_center_of_mass_iterative(gr2, rstart, :gas).gas.pos,
            gr2.gas.pos .- a;
            rtol=1e-4,
        ) |> all
        @test gr2.gas.pos == gr.gas.pos
        translate_to_center_of_mass_iterative!(gr2, rstart, :gas)
        @test isapprox.(gr2.gas.pos, gr.gas.pos .- a; rtol=1e-5) |> all
        gr2 = deepcopy(gr)
        @test isapprox.(
            translate_to_center_of_mass_iterative(gr2.gas, rstart).pos,
            gr2.gas.pos .- a;
            rtol=1e-4,
        ) |> all
        translate_to_center_of_mass_iterative!(gr2.gas, rstart)
        @test isapprox.(gr2.gas.pos, gr.gas.pos .- a; rtol=1e-5) |> all

        # center of velocity
        gr2 = deepcopy(gr)
        # a = [99.92409, -50.980724, -64.156525] .* u"km/s"
        a = [-2.912147, 3.523651, 13.41185] .* u"km/s"
        radius = 40u"kpc"
        @test isapprox.(center_of_velocity(gr2, radius, :gas; p=0.9), a; rtol=1e-5) |> all
        @test isapprox.(
            translate_to_center_of_velocity(gr2, radius, :gas).gas.vel,
            gr2.gas.vel .- a;
            rtol=1e-4,
        ) |> all
        @test gr2.gas.vel == gr.gas.vel
        translate_to_center_of_velocity!(gr2, radius, :gas)
        @test isapprox.(gr2.gas.vel, gr.gas.vel .- a; rtol=1e-4) |> all
        gr2 = deepcopy(gr)
        @test isapprox.(translate_to_center_of_velocity(gr2.gas, radius).vel, gr2.gas.vel .- a; rtol=1e-4) |>
              all
        translate_to_center_of_velocity!(gr2.gas, radius)
        @test isapprox.(gr2.gas.vel, gr.gas.vel .- a; rtol=1e-4) |> all
    end

    @testset "Global Properties" begin
        a = 14.7951u"kpc"
        b = 11.4378u"kpc"
        @test half_mass_radius(gr.gas.pos) ≈ a rtol = 1e-5
        @test half_mass_radius(gr.gas.pos, nothing) ≈ a rtol = 1e-5
        @test half_mass_radius(gr.gas.pos, gr.gas.mass) ≈ a rtol = 1e-5
        @test half_mass_radius(gr.gas.pos, gr.gas.mass; rmax=20u"kpc") ≈ b rtol = 1e-5
        @test half_mass_radius(gr.gas; rmax=20u"kpc") ≈ b rtol = 1e-5
        @test half_mass_radius(gr, :gas; rmax=20u"kpc") ≈ b rtol = 1e-5
        @test half_mass_radius(gr.dm; rmax=20u"kpc") ≈ 11.9408u"kpc" rtol = 1e-5

        c = 12.2065u"kpc"
        @test half_mass_radius_2D(gr.gas.pos) ≈ 8.31205u"kpc" rtol = 1e-5
        @test half_mass_radius_2D(gr.gas.pos, nothing; perspective=:sideon) ≈ c rtol = 1e-5
        @test half_mass_radius_2D(gr.gas.pos, gr.gas.mass; rmax=20u"kpc", perspective=:faceon) ≈ 8.32071u"kpc" rtol =
            1e-5
        @test half_mass_radius_2D(gr.gas; perspective=:sideon) ≈ c rtol = 1e-5
        @test half_mass_radius_2D(gr, :gas; perspective=:sideon) ≈ c rtol = 1e-5
        @test half_mass_radius_2D(gr.dm; rmax=20u"kpc") ≈ 9.79665u"kpc" rtol = 1e-5
    end

    @testset "Kinematic Properties" begin
        # total mass of particle component
        a = 1.65146e10u"Msun"
        @test GG.total_mass(gr.gas.pos, gr.gas.mass) ≈ a rtol = 1e-5
        @test GG.total_mass(gr.gas.pos, gr.gas.mass[1]) ≈ a rtol = 1e-5

        sph = Sphere(20u"kpc")

        # angular momentum
        # a = -8.0944f12u"Msun*kpc*km/s"
        a = -1.8541f11u"Msun*kpc*km/s"
        # b = -2.5619f12u"Msun*kpc*km/s"
        b = 1.10949f12u"Msun*kpc*km/s"
        @test angular_momentum(gr.gas.pos, gr.gas.vel, gr.gas.mass)[1] ≈ a rtol = 1e-4
        @test angular_momentum(gr.gas.pos, gr.gas.vel, gr.gas.mass, sph)[1] ≈ b rtol = 1e-4
        @test angular_momentum(gr.gas.pos, gr.gas.vel, gr.gas.mass[1], sph)[1] ≈ b rtol = 1e-4
        @test angular_momentum(gr.gas)[1] ≈ a rtol = 1e-4
        @test angular_momentum(gr.gas, sph)[1] ≈ b rtol = 1e-4

        # specific angular momentum
        # a = -490.137u"kpc*km/s"
        # b = -250.847u"kpc*km/s"
        a = -11.2274u"kpc*km/s"
        b = 108.6355u"kpc*km/s"
        @test specific_angular_momentum(gr.gas.pos, gr.gas.vel, gr.gas.mass)[1] ≈ a rtol = 1e-4
        @test specific_angular_momentum(gr.gas.pos, gr.gas.vel, gr.gas.mass, sph)[1] ≈ b rtol = 1e-4
        @test specific_angular_momentum(gr.gas.pos, gr.gas.vel, gr.gas.mass[1], sph)[1] ≈ b rtol = 1e-4
        @test specific_angular_momentum(gr.gas)[1] ≈ a rtol = 1e-4
        @test specific_angular_momentum(gr.gas; Mtot=sum(gr.gas.mass))[1] ≈ a rtol = 1e-4
        @test specific_angular_momentum(gr.gas, sph)[1] ≈ b rtol = 1e-4
        @test specific_angular_momentum(gr.gas, sph; Mtot=sum(gr.gas.mass))[1] ≈ b rtol = 1e-4

        # b-value
        # a = -3.9111
        # b = -3.9541
        a = -4.4449
        b = -4.2247
        @test b_value(gr.gas.pos, gr.gas.vel, gr.gas.mass) ≈ a rtol = 1e-4
        @test b_value(gr.gas.pos, gr.gas.vel, gr.gas.mass, sum(gr.gas.mass)) ≈ a rtol = 1e-4
        @test b_value(ustrip(gr.gas.pos), ustrip(gr.gas.vel), ustrip(gr.gas.mass)) ≈ a rtol = 1e-4
        @test b_value(gr.gas.pos, gr.gas.vel, gr.gas.mass, sph) ≈ b rtol = 1e-4
        @test b_value(gr.gas.pos, gr.gas.vel, gr.gas.mass[1], sph) ≈ b rtol = 1e-4
        @test b_value(gr.gas) ≈ a rtol = 1e-4
        @test b_value(gr.gas, sph) ≈ b rtol = 1e-4

        # circular velocity
        a = 103.693
        @test circular_velocity(2f10, 8) isa Float32
        @test circular_velocity(2e10, 8) isa Float64
        @test circular_velocity(2f10, 8.0) isa Float64
        @test circular_velocity(2, 8) isa Float64
        @test circular_velocity(2, 8.0) isa Float64
        @test circular_velocity(2u"Msun", 8u"kpc") isa Quantity{Float64}
        @test circular_velocity(2.0u"Msun", 8u"kpc") isa Quantity{Float64}
        @test circular_velocity(2.0f0u"Msun", 8u"kpc") isa Quantity{Float32}
        @test circular_velocity(2.0f0u"Msun", 8.0f0u"kpc") isa Quantity{Float32}
        @test circular_velocity(2.0e0u"Msun", 8.0f0u"kpc") isa Quantity{Float64}
        @test circular_velocity(2f10, 8) ≈ a rtol = 1e-5
        @test circular_velocity(2f10u"Msun", 8u"kpc") ≈ a * u"km/s" rtol = 1e-5

        # spin parameter
        r = 40u"kpc"
        # a = 0.333810
        # b = 0.310343
        # c = 0.356085
        a = 0.0976458
        b = 0.166437
        c = 0.170777
        @test spin_parameter(gr.gas.pos, gr.gas.vel, gr.gas.mass, sum(gr.gas.mass), r) ≈ a rtol = 1e-4
        @test spin_parameter(gr.gas.pos, gr.gas.vel, gr.gas.mass, sum(gr.gas.mass), sph) ≈ b rtol = 1e-4
        @test spin_parameter(gr.gas.pos, gr.gas.vel, gr.gas.mass[1], sum(gr.gas.mass), sph) ≈ b rtol = 1e-4
        @test spin_parameter(gr.gas, sum(gr.gas.mass), sph.radius) ≈ b rtol = 1e-4
        mask = gr.gas.pos[1, :] .< 3.5u"kpc"
        @test spin_parameter(gr.gas, sum(gr.gas.mass), sph.radius; mask) ≈ c rtol = 1e-4

        gr2 = deepcopy(gr)
        gr2.dm = gr2.gas
        gr2.stars = gr2.gas
        @test spin_parameter(gr2, 3sum(gr.gas.mass), r) ≈ 0.05637587 rtol = 1e-4
    end

    @info "deleting test data..."
    for i in 0:3
        rm("./sub_002.$i")
        rm("./snap_002.$i")
        rm("./snap_002.$i.key")
    end
    rm("./snap_002.key.index")
    @info "done!"
end


text:
# https://github.com/JuliaGraphics/Immerse.jl

supported_attrs(::ImmerseBackend) = supported_attrs(GadflyBackend())
supported_types(::ImmerseBackend) = supported_types(GadflyBackend())
supported_styles(::ImmerseBackend) = supported_styles(GadflyBackend())
supported_markers(::ImmerseBackend) = supported_markers(GadflyBackend())
supported_scales(::ImmerseBackend) = supported_scales(GadflyBackend())
is_subplot_supported(::ImmerseBackend) = true

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

function _initialize_backend(::ImmerseBackend; kw...)
  @eval begin
    import Immerse, Gadfly, Compose, Gtk
    export Immerse, Gadfly, Compose, Gtk
    include(joinpath(dirname(@__FILE__), "gadfly_shapes.jl"))
  end
end

function createImmerseFigure(plotattributes::KW)
  w,h = plotattributes[:size]
  figidx = Immerse.figure(; name = plotattributes[:window_title], width = w, height = h)
  Immerse.Figure(figidx)
end

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

# create a blank Gadfly.Plot object
# function _create_plot(pkg::ImmerseBackend, plotattributes::KW)
#   # create the underlying Gadfly.Plot object
#   gplt = createGadflyPlotObject(plotattributes)
#
#   # save both the Immerse.Figure and the Gadfly.Plot
#   Plot((nothing,gplt), pkg, 0, plotattributes, KW[])
# end
function _create_backend_figure(plt::Plot{ImmerseBackend})
    (nothing, createGadflyPlotObject(plt.attr))
end

# # plot one data series
# function _series_added(::ImmerseBackend, plt::Plot, plotattributes::KW)
#   addGadflySeries!(plt, plotattributes)
#   push!(plt.seriesargs, plotattributes)
#   plt
# end

function _series_added(plt::Plot{ImmerseBackend}, series::Series)
    addGadflySeries!(plt, series.plotattributes)
end

function _update_plot_object(plt::Plot{ImmerseBackend}, plotattributes::KW)
  updateGadflyGuides(plt, plotattributes)
  updateGadflyPlotTheme(plt, plotattributes)
end

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

function _add_annotations(plt::Plot{ImmerseBackend}, anns::AVec{Tuple{X,Y,V}}) where {X,Y,V}
  for ann in anns
    push!(getGadflyContext(plt).guides, createGadflyAnnotationObject(ann...))
  end
end

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

# accessors for x/y data

function getxy(plt::Plot{ImmerseBackend}, i::Integer)
  mapping = getGadflyMappings(plt, i)[1]
  mapping[:x], mapping[:y]
end

function setxy!(plt::Plot{ImmerseBackend}, xy::Tuple{X,Y}, i::Integer) where {X,Y}
  for mapping in getGadflyMappings(plt, i)
    mapping[:x], mapping[:y] = xy
  end
  plt
end

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

# function _create_subplot(subplt::Subplot{ImmerseBackend}, isbefore::Bool)
#   return false
#   # isbefore && return false
# end
#
# function showSubplotObject(subplt::Subplot{ImmerseBackend})
#   # create the Gtk window with vertical box vsep
#   plotattributes = getattr(subplt,1)
#   w,h = plotattributes[:size]
#   vsep = Gtk.GtkBoxLeaf(:v)
#   win = Gtk.GtkWindowLeaf(vsep, plotattributes[:window_title], w, h)
#
#   figindices = []
#   row = Gtk.GtkBoxLeaf(:h)
#   push!(vsep, row)
#   for (i,(r,c)) in enumerate(subplt.layout)
#     plt = subplt.plts[i]
#
#     # get the components... box is the main plot GtkBox, and canvas is the GtkCanvas where it's plotted
#     box, toolbar, canvas = Immerse.createPlotGuiComponents()
#
#     # add the plot's box to the row
#     push!(row, box)
#
#     # create the figure and store the index returned for destruction later
#     figidx = Immerse.figure(canvas)
#     push!(figindices, figidx)
#
#     fig = Immerse.figure(figidx)
#     plt.o = (fig, plt.o[2])
#
#     # add the row
#     if c == ncols(subplt.layout, r)
#       row = Gtk.GtkBoxLeaf(:h)
#       push!(vsep, row)
#     end
#
#   end
#
#   # destructor... clean up plots
#   Gtk.on_signal_destroy((x...) -> ([Immerse.dropfig(Immerse._display,i) for i in figindices]; subplt.o = nothing), win)
#
#   subplt.o = win
#   true
# end

function _remove_axis(plt::Plot{ImmerseBackend}, isx::Bool)
  gplt = getGadflyContext(plt)
  addOrReplace(gplt.guides, isx ? Gadfly.Guide.xticks : Gadfly.Guide.yticks; label=false)
  addOrReplace(gplt.guides, isx ? Gadfly.Guide.xlabel : Gadfly.Guide.ylabel, "")
end

function _expand_limits(lims, plt::Plot{ImmerseBackend}, isx::Bool)
  for l in getGadflyContext(plt).layers
    _expand_limits(lims, l.mapping[isx ? :x : :y])
  end
end

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

getGadflyContext(plt::Plot{ImmerseBackend}) = plt.o[2]
# getGadflyContext(subplt::Subplot{ImmerseBackend}) = buildGadflySubplotContext(subplt)

function Base.display(::PlotsDisplay, plt::Plot{ImmerseBackend})

  fig, gplt = plt.o
  if fig == nothing
    fig = createImmerseFigure(plt.attr)
    Gtk.on_signal_destroy((x...) -> (Immerse.dropfig(Immerse._display, fig.figno); plt.o = (nothing,gplt)), fig.canvas)
    plt.o = (fig, gplt)
  end

  Immerse.figure(fig.figno; displayfig = false)
  display(gplt)
end

# function Base.display(::PlotsDisplay, subplt::Subplot{ImmerseBackend})
#
#   # if we haven't created the window yet, do it
#   if subplt.o == nothing
#     showSubplotObject(subplt)
#   end
#
#   # display the plots by creating a fresh Immerse.Figure object from the GtkCanvas and Gadfly.Plot
#   for plt in subplt.plts
#     fig, gplt = plt.o
#     Immerse.figure(fig.figno; displayfig = false)
#     display(gplt)
#   end
#
#   # o is the window... show it
#   showall(subplt.o)
# end


text:
using Revise, Pkg, DataFrames, Finometrics, Random

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

Pkg.activate(pwd())
using Production

Random.seed!(1111)

if @isdefined(DEFINED)
  revise(Production)
end

const REDEFINE = true
#this is to avoid annoying redefinition warnings when possible
if REDEFINE || !(@isdefined DEFINED)
  const STATES_PER_VAR = 60

  const N_SIMS = 200
  const BURN_IN = 10^5

  const TOL = 1e-10
  const D_TOL = 1e-8
  const MAX_ITER = 10000

  const τ_SCENARIOS = [[0.0, 0.0], [0.4, 0.0], [0.0, -0.4], [0.2, -0.2]]
  const τ_SCENARIO_NAMES = [:baseline, :tax40, :subsidy40, :tax20subsidy20]
  const τ_VALS = collect(-0.8:0.05:0.8)

  #***************competitive parameters**************************
  const COM_τ = [0.0, 0.0]
  const COM_δ = 0.3
  const COM_Θ = (0.6, 0.4)
  const COM_α = 1.05
  const COM_β = 0.98
  const COM_γ = Float64(1/3)
  const COM_ζ = (1.0, 0.0)

  const COM_p = [0.81 0.68; 0.19 0.32]
  const COM_M = 2

  #some heuristics for the maximum values
  const COM_KMAX = (COM_Θ[1]/(COM_p[1,1]^(1/COM_α)*COM_δ))^(1/(1-COM_γ))
  const COM_YMAX = COM_KMAX^COM_γ*COM_Θ[1]/COM_p[1,2]^(1/COM_α)

  const COM_YHSTATES = collect(range(TOL, stop=COM_YMAX, length=STATES_PER_VAR))
  const COM_KSTATES = collect(range(TOL, stop=COM_KMAX, length=STATES_PER_VAR))

  #***************monopolistic parameters**************************
  const MON_τ = COM_τ
  const MON_δ = COM_δ
  const MON_Θ = COM_Θ
  const MON_α = COM_α
  const MON_β = COM_β
  const MON_γ = COM_γ
  const MON_ζ = (3.0, 1.0)

  const MON_p = COM_p
  const MON_M = COM_M

  #technically this sould be based on a sepearte
  const MON_YHSTATES = COM_YHSTATES
  const MON_KSTATES = COM_KSTATES

end
const DEFINED = true

function productionscript(firmtype::Symbol)::Nothing
  local F::Firm
  local Φ::Problem

  #the problem parameters
  local τ::Vector{Float64}
  local δ::Float64
  local Θ::NTuple{2,Float64}
  local α::Float64
  local β::Float64
  local γ::Float64
  local kstates::Vector{Float64}
  local yhstates::Vector{Float64}
  local ζ::NTuple{2,Float64}
  local p::Matrix{Float64}

  #assign values based on the scenario
  if firmtype ∈ (:competitive, :COM, :com)
    τ= COM_τ
    δ= COM_δ
    Θ= COM_Θ
    α= COM_α
    β= COM_β
    γ= COM_γ
    kstates= COM_KSTATES
    yhstates= COM_YHSTATES
    ζ= COM_ζ
    p= COM_p
  elseif firmtype ∈ (:monpolistic, :MON, :mon)
    τ = MON_τ
    δ = MON_δ
    Θ = MON_Θ
    α = MON_α
    β= MON_β
    γ= MON_γ
    kstates = MON_KSTATES
    yhstates = MON_YHSTATES
    ζ= MON_ζ
    p = MON_p
  else
    error("firmtype $firmtype not recognized")
  end

  println("\n#**********RUNNING $firmtype using $STATES_PER_VAR states/var***********#")
  println("\nMax K: $(maximum(kstates)) Max Yh: $(maximum(yhstates))")
  F = Firm(τ=τ, δ=δ, Θ=Θ, α=α, β=β, γ=γ, kstates=kstates, yhstates=yhstates, ζ=ζ, p=p)
  Φ = Problem(F, p, dtol = D_TOL, maxiter = MAX_ITER, verbose=false, findpolicy=false)
  Ψ::Simulation = Simulation(Φ, nsims= N_SIMS, burnin = BURN_IN, simulate = false)

  graphscript!(Ψ, τscenarios = τ_SCENARIOS, τscenarionames = τ_SCENARIO_NAMES, prefix="$firmtype",
    τvals = τ_VALS, scenariographs=true, taxanalysis=false)
  return nothing
end

@time productionscript(:COM)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "AABox" begin
    @testset "AABox{2}" begin
        for T in Floats
            P₁ = Point{2,T}(1, 0)
            P₂ = Point{2,T}(3, 2)

            # getproperty
            aab = AABox(P₁, P₂)
            @test xmin(aab) ≈ 1
            @test ymin(aab) ≈ 0
            @test xmax(aab) ≈ 3
            @test ymax(aab) ≈ 2

            # Δx, Δy
            @test Δx(aab) ≈ 2
            @test Δy(aab) ≈ 2

            # measure
            @test measure(aab) ≈ 4
        end
    end
    @testset "AABox{3}" begin
        for T in Floats
            # getproperty
            aab = AABox(Point{3,T}(1, 0, -1), Point{3,T}(3, 2, 1))
            @test xmin(aab) ≈ 1
            @test ymin(aab) ≈ 0
            @test zmin(aab) ≈ -1
            @test xmax(aab) ≈ 3
            @test ymax(aab) ≈ 2
            @test zmax(aab) ≈ 1

            # Δx, Δy, Δz
            @test Δx(aab) ≈ 2
            @test Δy(aab) ≈ 2
            @test Δz(aab) ≈ 2

            # measure
            @test measure(aab) ≈ 8
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
module CommunalHelperConnectedZipMover

using ..Ahorn, Maple
using Ahorn.CommunalHelper

const theme_options = String["Normal", "Moon", "Cliffside"]

@mapdef Entity "CommunalHelper/ConnectedZipMover" ConnectedZipMover(
    x::Integer,
    y::Integer,
    width::Integer=16,
    height::Integer=16,
    theme::String="Normal",
    permanent::Bool=false,
    waiting::Bool=false,
    ticking::Bool=false,
    nodes::Array{Tuple{Integer,Integer},1}=Tuple{Integer,Integer}[],
    customSkin::String="",
    colors::String="",
)

const placements = Ahorn.PlacementDict(
    "Connected Zip Mover ($(uppercasefirst(theme))) (Communal Helper)" => Ahorn.EntityPlacement(
        ConnectedZipMover,
        "rectangle",
        Dict{String,Any}(
            "theme" => theme,
        ),
        function (entity)
            entity.data["nodes"] = [(
                Int(entity.data["x"]) + Int(entity.data["width"]) + 8,
                Int(entity.data["y"]),
            )]
        end,
    ) for theme in theme_options
)

placements["Connected Zip Mover (Reskinnable) (Communal Helper)"] = Ahorn.EntityPlacement(
    ConnectedZipMover,
    "rectangle",
    Dict{String,Any}(
        "customBlockTexture" => "CommunalHelper/customConnectedBlock/customConnectedBlock",
    ),
)

Ahorn.editingOptions(entity::ConnectedZipMover) = Dict{String,Any}(
    "theme" => theme_options,
)

Ahorn.nodeLimits(entity::ConnectedZipMover) = 1, -1

Ahorn.minimumSize(entity::ConnectedZipMover) = 16, 16
Ahorn.resizable(entity::ConnectedZipMover) = true, true

function Ahorn.selection(entity::ConnectedZipMover)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    res = [Ahorn.Rectangle(x, y, width, height)]

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        push!(res, Ahorn.Rectangle(nx + floor(Int, width / 2) - 5, ny + floor(Int, height / 2) - 5, 10, 10))
    end

    return res
end

function getTextures(entity::ConnectedZipMover)
    theme = lowercase(get(entity, "theme", "normal"))
    isCliffside = theme == "cliffside"
    folder = isCliffside ? "CommunalHelper/connectedZipMover" : "zipmover"
    themePath = (theme == "normal") ? "" : string(theme, "/")

    customSkin = get(entity, "customSkin", "")
    if customSkin != ""
        return (
            "$(customSkin)/block",
            "$(customSkin)/light01",
            "$(customSkin)/cog",
            "$(customSkin)/innerCorners"
        )
    end

    return (
        "objects/$(folder)/$(themePath)block",
        "objects/$(folder)/$(themePath)light01",
        "objects/$(folder)/$(themePath)cog",
        "objects/$((isCliffside ? "" : "CommunalHelper/") * folder)/$(themePath)innerCorners"
    )
end

const ropeColor = (102, 57, 49) ./ 255

function renderZipMover(ctx::Ahorn.Cairo.CairoContext, entity::ConnectedZipMover, room::Maple.Room)
    rects = getExtensionRectangles(room)

    x, y = Ahorn.position(entity)
    px, py = x, y

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    block, light, cog, incorners = getTextures(entity)
    lightSprite = Ahorn.getSprite(light, "Gameplay")

    customBlockTexture = String(get(entity.data, "customBlockTexture", ""))
    hasCustomTexture = customBlockTexture != ""
    txOffset = 0

    if hasCustomTexture
        block = incorners = "objects/" * customBlockTexture
        txOffset = 24
    end

    tileWidth = div(width, 8)
    tileHeight = div(height, 8)

    rect = Ahorn.Rectangle(x, y, width, height)

    if !(rect in rects)
        push!(rects, rect)
    end

    # Iteration through all the nodes
    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)
        cx, cy = px + width / 2, py + height / 2
        cnx, cny = nx + width / 2, ny + height / 2

        length = sqrt((px - nx)^2 + (py - ny)^2)
        theta = atan(cny - cy, cnx - cx)

        Ahorn.Cairo.save(ctx)

        Ahorn.translate(ctx, cx, cy)
        Ahorn.rotate(ctx, theta)

        Ahorn.setSourceColor(ctx, ropeColor)
        Ahorn.set_antialias(ctx, 1)
        Ahorn.set_line_width(ctx, 1)

        # Offset for rounding errors
        Ahorn.move_to(ctx, 0, 4 + (theta <= 0))
        Ahorn.line_to(ctx, length, 4 + (theta <= 0))

        Ahorn.move_to(ctx, 0, -4 - (theta > 0))
        Ahorn.line_to(ctx, length, -4 - (theta > 0))

        Ahorn.stroke(ctx)

        Ahorn.Cairo.restore(ctx)

        Ahorn.drawSprite(ctx, cog, cnx, cny)

        px, py = nx, ny
    end

    Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (0.0, 0.0, 0.0, 1.0))

    for i in 1:tileWidth, j in 1:tileHeight
        drawX, drawY = (i - 1) * 8, (j - 1) * 8

        closedLeft = !notAdjacent(entity, drawX - 8, drawY, rects)
        closedRight = !notAdjacent(entity, drawX + 8, drawY, rects)
        closedUp = !notAdjacent(entity, drawX, drawY - 8, rects)
        closedDown = !notAdjacent(entity, drawX, drawY + 8, rects)
        completelyClosed = closedLeft && closedRight && closedUp && closedDown

        if completelyClosed
            Ahorn.drawRectangle(ctx, x + drawX, y + drawY, 8, 8, (0.0, 0.0, 0.0, 1.0))
            if notAdjacent(entity, drawX + 8, drawY - 8, rects)
                # up right
                Ahorn.drawImage(ctx, incorners, x + drawX, y + drawY, 8 + txOffset, 0, 8, 8)

            elseif notAdjacent(entity, drawX - 8, drawY - 8, rects)
                # up left
                Ahorn.drawImage(ctx, incorners, x + drawX, y + drawY, 0 + txOffset, 0, 8, 8)

            elseif notAdjacent(entity, drawX + 8, drawY + 8, rects)
                # down right
                Ahorn.drawImage(ctx, incorners, x + drawX, y + drawY, 8 + txOffset, 8, 8, 8)

            elseif notAdjacent(entity, drawX - 8, drawY + 8, rects)
                # down left
                Ahorn.drawImage(ctx, incorners, x + drawX, y + drawY, 0 + txOffset, 8, 8, 8)
            end
        else
            if closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawRectangle(ctx, x + drawX, y + drawY, 8, 8, (0.0, 0.0, 0.0, 1.0))
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 8, 0, 8, 8)

            elseif closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawRectangle(ctx, x + drawX, y + drawY, 8, 8, (0.0, 0.0, 0.0, 1.0))
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 8, 16, 8, 8)

            elseif closedLeft && !closedRight && closedUp && closedDown
                Ahorn.drawRectangle(ctx, x + drawX, y + drawY, 8, 8, (0.0, 0.0, 0.0, 1.0))
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 16, 8, 8, 8)

            elseif !closedLeft && closedRight && closedUp && closedDown
                Ahorn.drawRectangle(ctx, x + drawX, y + drawY, 8, 8, (0.0, 0.0, 0.0, 1.0))
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 0, 8, 8, 8)

            elseif closedLeft && !closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 16, 0, 8, 8)

            elseif !closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 0, 0, 8, 8)

            elseif !closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 0, 16, 8, 8)

            elseif closedLeft && !closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 16, 16, 8, 8)
            end
        end
    end

    Ahorn.drawImage(ctx, lightSprite, x + floor(Int, (width - lightSprite.width) / 2), y)
end

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ConnectedZipMover, room::Maple.Room) = renderZipMover(ctx, entity, room)

end