text:
function BesselJ(arg0::jdouble)
    return BesselJ((jdouble,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function rj_besl(::Type{BesselJ}, arg0::jdouble, arg1::jdouble, arg2::jint)
    return jcall(BesselJ, "rjBesl", BesselJ_BesselJResult, (jdouble, jdouble, jint), arg0, arg1, arg2)
end

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

function value(::Type{BesselJ}, arg0::jdouble, arg1::jdouble)
    return jcall(BesselJ, "value", jdouble, (jdouble, jdouble), arg0, arg1)
end

function value(obj::BesselJ, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using Test
using AutomotiveDrivingModels
using DataFrames
using Distributions

@testset "Vec" begin
    include("vec-tests/vec_runtests.jl")
end

@testset "AutomotiveDrivingModels" begin
    include("test_roadways.jl")
    include("test_agent_definitions.jl")
    include("test_scenes.jl")
    include("test_states.jl")
    include("test_collision_checkers.jl")
    include("test_actions.jl")
    include("test_features.jl")
    include("test_behaviors.jl")
    include("test_simulation.jl")
end

# Restore after AutoViz v0.9 is tagged
# @testset "doc tutorials" begin
#     @testset "basics" begin
#         include("../docs/lit/tutorials/straight_roadway.jl")
#     end
#     @testset "cameras" begin
#         include("../docs/lit/tutorials/stadium.jl")
#     end
#     @testset "overlays" begin
#         include("../docs/lit/tutorials/intersection.jl")
#     end
#     @testset "basics" begin
#         include("../docs/lit/tutorials/crosswalk.jl")
#     end
#     @testset "basics" begin
#         include("../docs/lit/tutorials/sidewalk.jl")
#     end
# end

ext: jl


text:
using Random: MersenneTwister

funcs = ( dof,
          dof_residual,
          confint,
          deviance,
          nulldeviance,
          loglikelihood,
          nullloglikelihood,
          dispersion,
          nobs,
          stderror,
          vcov,
          residuals,
          predict,
          response,
          weights,
          leverage,
          modelmatrix,
          dispersion,
          scale
        )

m1 = fit(LinearModel, form, data)

loss1 = RobustModels.L2Loss()
loss2 = RobustModels.TukeyLoss()
est1 = MEstimator(loss1)
est2 = MEstimator(loss2)

@testset "linear: Ridge M-estimator $(lossname)" for lossname in ("L2", "Huber", "Tukey")
    typeloss = getproperty(RobustModels, Symbol(lossname * "Loss"))
    l = typeloss()
    est = MEstimator(typeloss())

    # Formula, dense and sparse entry  and methods :cg and :chol
    @testset "(type, method): ($(typeof(A)),\t$(method))" for (A, b) in ((form, data), (X, y), (sX, y)), method in (:cg, :chol)
        aspace = if method==:cg; "  " else "    " end
        name  = "MEstimator($(typeloss)),\t"
        name *= if A==form; "formula" elseif A==X; "dense  " else "sparse " end
        name *= if method==:cg; ",  cg" else ",chol" end
        # use the dispersion from GLM to ensure that the loglikelihood is correct
        m2 = fit(RobustLinearModel, A, b, est; method=method, initial_scale=:L1)
        m3 = fit(RobustLinearModel, A, b, est; method=method, initial_scale=:L1, ridgeλ=1)
        m4 = fit(RobustLinearModel, A, b, est; method=method, initial_scale=:L1, ridgeλ=1, ridgeG=float([0 0; 0 1]))
        m5 = fit(RobustLinearModel, A, b, est; method=method, initial_scale=:L1, ridgeλ=0.1, ridgeG=[0, 1])
        m6 = fit(RobustLinearModel, A, b, est; method=method, initial_scale=:L1, ridgeλ=1, ridgeG=[0, 1], βprior=[0.0, 2.0])
        println("\n\t\u25CF Estimator: $(name)")
        println(" lm$(aspace)               : ", coef(m1))
        println("rlm($(method))             : ", coef(m2))
        println("ridge λ=1   rlm3($(method)): ", coef(m3))
        println("ridge λ=1   rlm4($(method)): ", coef(m4))
        println("ridge λ=0.1 rlm5($(method)): ", coef(m5))
        println("ridge λ=1 βprior=[0,2] rlm6($(method)): ", coef(m6))
        @test isapprox(coef(m3), coef(m4); rtol=1e-5)

        # Test printing the model
        @test_nowarn println(m3)

        # refit
        refit!(m5; ridgeλ=1, initial_scale=:L1)
        @test isapprox(m5.pred.λ, 1.0; rtol=1e-5)
        @test isapprox(coef(m3), coef(m5); rtol=1e-5)
    end
end

@testset "linear: Ridge L2 estimator methods" begin
    m2 = fit(RobustLinearModel, form, data, est1; method=:chol, initial_scale=:L1)
    m3 = fit(RobustLinearModel, form, data, est1; method=:chol, initial_scale=:L1, ridgeλ=1)

    @testset "method: $(f)" for f in funcs
        # make sure the interfaces for RobustLinearModel are well defined
        @test_nowarn f(m3)
    end

    # Check coefficients and dof change
    λs = vcat([0], 10 .^ range(-2, 1, length=4))

    L = length(coef(m3))
    βs = zeros(L, 5)
    dofs = zeros(5)
    for (i,λ) in enumerate(λs)
        refit!(m3; ridgeλ=λ, initial_scale=:L1)
        βs[:, i] = coef(m3)
        dofs[i] = dof(m2)
    end

    @test isapprox(dofs[1], dof(m2); rtol=1e-5)
    @test isapprox(βs[:, 1], coef(m2); rtol=1e-5)
    @test sort(dofs; rev=true)==dofs
    @test all(sort(βs[r, :]; rev=true)==βs[r, :] for r in 2:L)  # not the intercept
end


text:
# v0.6

function pascal(n::Int)
    r = ones(Int, n, n)
    for i in 2:n, j in 2:n
        r[i, j] = r[i-1, j] + r[i, j-1]
    end
    return r
end

function catalan_num(n::Int)
    p = pascal(n + 2)
    p[n+4:n+3:end-1] - diag(p, 2)
end

@show catalan_num(15)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
##
using Pkg, Revise, BenchmarkTools
Pkg.activate(".")
##
import Images
import NPZ: npzread
using GLMakie
using ZStackAnalysis
##
using Base.Threads
nthreads()
##
file_content = npzread("./data/Fed_X63_Z2_SIM.npz")
zstack = Preprocessing.convert_to_float32(file_content["data"])
scale = file_content["scale"]
##
size(zstack)
##

##
xy_region = (1201:1500, 201:500)
actin = interpolate_zstack(zstack[1, :, xy_region...], scale)
actin_adj = Images.adjust_histogram(actin, Images.Equalization(nbins = 256, minval = 0.0, maxval = 0.7))
##
intensity, orientation = LineFilterTransform(actin_adj[1:20, 1:20, 1:20], 10, 7, 20)
@btime OrientationFilterTransform(intensity, orientation, 10, 7, 20)
##
@btime line_segment()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale <: SwaggerModel
    labelSelectorPath::Any # spec type: Union{ Nothing, String } # spec name: labelSelectorPath
    specReplicasPath::Any # spec type: Union{ Nothing, String } # spec name: specReplicasPath
    statusReplicasPath::Any # spec type: Union{ Nothing, String } # spec name: statusReplicasPath

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale(;labelSelectorPath=nothing, specReplicasPath=nothing, statusReplicasPath=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale, Symbol("labelSelectorPath"), labelSelectorPath)
        setfield!(o, Symbol("labelSelectorPath"), labelSelectorPath)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale, Symbol("specReplicasPath"), specReplicasPath)
        setfield!(o, Symbol("specReplicasPath"), specReplicasPath)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale, Symbol("statusReplicasPath"), statusReplicasPath)
        setfield!(o, Symbol("statusReplicasPath"), statusReplicasPath)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale = Dict{Symbol,Symbol}(Symbol("labelSelectorPath")=>Symbol("labelSelectorPath"), Symbol("specReplicasPath")=>Symbol("specReplicasPath"), Symbol("statusReplicasPath")=>Symbol("statusReplicasPath"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale = Dict{Symbol,String}(Symbol("labelSelectorPath")=>"String", Symbol("specReplicasPath")=>"String", Symbol("statusReplicasPath")=>"String")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale)
    (getproperty(o, Symbol("specReplicasPath")) === nothing) && (return false)
    (getproperty(o, Symbol("statusReplicasPath")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresourceScale }, name::Symbol, val)
end


text:
# wls object is not actively maintained
"""
# User provided attributes
- `inputs`: the vector or array view of input points
- `outputs`: the vector or array view of output scalars
- `EPS::Float64`: ε of the method
- `weightFunc::Function`: weighting function of the method

# Automatically created attributes
If created by the `wls` function, the attributes in this section are created automatically.

- `vars::Vector{PolyVar{true}}`: variables of the polynomial,
- `b::Vector{Monomial{true}}`: basis of the polynomial,
- `matrix`: the result of `b * transpose(b)`.
"""
struct WlsObject
  inputs
  outputs
  EPS::Float64
  weightFunc::Function
  vars::Vector{PolyVar{true}}
  b::Vector{Monomial{true}}
  matrix
end

function WlsObject(inputs, outputs, EPS, weightFunc, vars, b)
  WlsObject(inputs, outputs, EPS, weightFunc, vars, b, b * transpose(b))
end

"""
    wls(inputs::Array{Float64, 2}, outputs::Vector{Float64}, maxDegree::Int, EPS::Float64, weightFunc::Function)

# Arguments
- `inputs`: a 2d array where each point is on a single row,
- `outputs`: an array of output scalars,
- `maxDegree::Int64`: the maximal degree of each polynomial term in the method,
- `EPS::Float64`: ε of the method,
- `weightFunc::Function`: weighting function of the method. It should be in form of `(distance, EPS) -> Float64`.
"""
function wls(inputs, outputs, maxDegree::Int, EPS::Float64, weightFunc::Function)
  inputDim = size(inputs, 2)
  vars::Vector{PolyVar{true}}, b::Vector{Monomial{true}} = generateMonomials(inputDim, maxDegree)
  return WlsObject(inputs, outputs, EPS, weightFunc, vars, b)
end


text:
## Learning a connect four agent for the JuliaCon 2020 demo

ENV["JULIA_CUDA_MEMORY_POOL"] = "split" # "binned" / "split"

using AlphaZero

include("../games/connect-four/main.jl")
using .ConnectFour: Game, Training

session = Session(
  Game,
  Training.Network{Game},
  Training.params,
  Training.netparams,
  benchmark=Training.benchmark,
  dir="sessions/connect-four",
  autosave=true,
  save_intermediate=false)

resume!(session)


text:
using SparseArrays
using LinearAlgebra

function setupFWI(m,filenamePrefix::String,plotting::Bool,
		workersFWI::Array{Int64,1}=workers(),maxBatchSize::Int64 = 48,
		Ainv::AbstractSolver = getJuliaSolver(), misfun::Function=SSDFun,useFilesForFields::Bool = false,useFreqOnlySplit::Bool = true)

file = matopen(string(filenamePrefix,"_PARAM.mat"));
n_cells = read(file,"n");
OmegaDomain = read(file,"domain");
Minv = getRegularMesh(OmegaDomain,n_cells);
gamma = read(file,"gamma");
omega = read(file,"omega");
waveCoef = read(file,"waveCoef");
if length(omega)==1 # matlab saves a 1 variable array as scalar.
	omega = [omega];
	waveCoef = [waveCoef];
end

boundsLow = read(file,"boundsLow");
boundsHigh = read(file,"boundsHigh");

mref =  read(file,"mref");
close(file);

### Read receivers and sources files
RCVfile = string(filenamePrefix,"_rcvMap.dat");
SRCfile = string(filenamePrefix,"_srcMap.dat");

srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);

Q = generateSrcRcvProjOperators(Minv.n.+1,srcNodeMap);
Q = Q.*1.0./(norm(Minv.h)^2);
println("We have ",size(Q,2)," sources");
P = generateSrcRcvProjOperators(Minv.n.+1,rcvNodeMap);

########################################################################################################
##### Set up remote workers ############################################################################
########################################################################################################

N = prod(Minv.n);

Iact = SparseMatrixCSC(1.0I, N, N);
mback   = zeros(Float64,N);
## Setting the sea constant:
mask = zeros(N);
mask[abs.(m[:] .- minimum(m)) .< 1e-2] .= 1.0;
mask[gamma[:] .>= 0.95*maximum(gamma)] .= 1.0;
# setup active cells
mback = vec(mref[:].*mask);
sback = velocityToSlowSquared(mback)[1];
sback[mask .== 0.0] .= 0.0;
Iact = Iact[:,mask .== 0.0];

boundsLow = Iact'*boundsLow;
boundsHigh = Iact'*boundsHigh;
mref = Iact'*mref[:];

####################################################################################################################
####################################################################################################################

println("Reading FWI data:");

batch = min(size(Q,2),maxBatchSize);
if useFreqOnlySplit
	(pForFWI,contDivFWI,SourcesSubIndFWI) = getFWIparamFreqOnlySplit(omega,waveCoef,vec(gamma),Q,P,Minv,Ainv,workersFWI,batch,useFilesForFields);
else
	(pForFWI,contDivFWI,SourcesSubIndFWI) = getFWIparam(omega,waveCoef,vec(gamma),Q,P,Minv,Ainv,workersFWI,batch,useFilesForFields);
end
# write data to remote workers
Wd   = Array{Array{ComplexF64,2}}(undef,length(pForFWI))
dobs = Array{Array{ComplexF64,2}}(undef,length(pForFWI))
for k = 1:length(omega)
	omRound = string(round((omega[k]/(2*pi))*100.0)/100.0);
	(DobsFWIwk,WdFWIwk) =  readDataFileToDataMat(string(filenamePrefix,"_freq",omRound,".dat"),srcNodeMap,rcvNodeMap);

	if useFreqOnlySplit
		Wd[k] = WdFWIwk
		dobs[k] = DobsFWIwk
	else
		for i = contDivFWI[k]:contDivFWI[k+1]-1
			I_i = SourcesSubIndFWI[i]; # subset of sources for ith worker.
			Wd[i] 	= WdFWIwk[:,I_i];
			dobs[i] = DobsFWIwk[:,I_i];
		end
	end
	DobsFWIwk = 0;
	WdFWIwk = 0;
end

pMisFWIRFs = getMisfitParam(pForFWI, Wd, dobs, misfun, Iact,sback);

########################################################################################################
##### Set up remote workers ############################################################################
########################################################################################################

SourcesSubInd    = SourcesSubIndFWI;
pMis 			 = pMisFWIRFs;
contDiv			 = contDivFWI;

return Q,P,pMis,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow,resultsFilename;
end


text:
using MendelImpute
using Test
using Random
using JLSO
using VCFTools
using CSV
using Statistics
using BenchmarkTools
using DataFrames
using VariantCallFormat
using BGEN

# change directory to where data is located
cd(normpath(MendelImpute.datadir()))

@testset "search breakpoints" begin
    X1 = Vector{Union{Float64, Missing}}([0,0,0,1,1,1])
    h1 = zeros(6)
    h2 = zeros(6)
    h3 = ones(6)
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X1, h1, h2, h3)
    @test bkpt_optim == 3
    @test err_optim == 0

    X2 = Vector{Union{Float64, Missing}}([0,0,1,1,1,1])
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X2, h1, h2, h3)
    @test bkpt_optim == 2
    @test err_optim == 0

    X3 = Vector{Union{Float64, Missing}}([missing,0,1,1,1,1])
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X3, h1, h2, h3)
    @test bkpt_optim == 2
    @test err_optim == 0

    X4 = Vector{Union{Float64, Missing}}([0,0,0,missing,1,2])
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X4, h1, h2, h3)
    @test bkpt_optim == 3
    @test err_optim == 1
end

@testset "intersect code" begin
    x = [1, 3, 4, 5, 7, 9]
    y = [2, 3, 5, 6]
    @test MendelImpute.intersect_size_sorted(x, y) == 2
    MendelImpute.intersect_lange!(x, y)
    @test all(x .== [3, 5])
    @test all(y .== [2, 3, 5, 6])

    x = [1, 3, 4, 7]
    y = [2, 3, 5, 6, 7, 10]
    @test MendelImpute.intersect_size_sorted(x, y) == 2
    MendelImpute.intersect_lange!(x, y)
    @test all(x .== [3, 7])
    @test all(y .== [2, 3, 5, 6, 7, 10])

    # allow repeats, although we don't have any in MendelImpute
    x = [3, 4, 7, 7, 7, 10]
    y = [2, 3, 5, 7, 7, 10]
    @test MendelImpute.intersect_size_sorted(x, y) == 4
    MendelImpute.intersect_lange!(x, y)
    @test all(x .== [3, 7, 7, 10])
    @test all(y .== [2, 3, 5, 7, 7, 10])
end

@testset "dynamic programming" begin
    T = Tuple{Int32, Int32}
    windows = 4
    haplotype_set = [T[] for i in 1:windows]
    push!(haplotype_set[1], (1, 4))
    push!(haplotype_set[1], (1, 5))
    push!(haplotype_set[1], (1, 6))
    push!(haplotype_set[1], (2, 4))
    push!(haplotype_set[1], (2, 5))
    push!(haplotype_set[1], (2, 6))
    push!(haplotype_set[1], (3, 4))
    push!(haplotype_set[1], (3, 5))
    push!(haplotype_set[1], (3, 6))

    push!(haplotype_set[2], (1, 5))
    push!(haplotype_set[2], (1, 7))
    push!(haplotype_set[2], (1, 8))
    push!(haplotype_set[2], (2, 5))
    push!(haplotype_set[2], (2, 7))
    push!(haplotype_set[2], (2, 8))
    push!(haplotype_set[2], (6, 5))
    push!(haplotype_set[2], (6, 7))
    push!(haplotype_set[2], (6, 8))

    push!(haplotype_set[3], (1, 1))
    push!(haplotype_set[3], (1, 5))
    push!(haplotype_set[3], (3, 1))
    push!(haplotype_set[3], (3, 5))

    push!(haplotype_set[4], (4, 4))
    push!(haplotype_set[4], (4, 6))
    push!(haplotype_set[4], (4, 8))
    push!(haplotype_set[4], (5, 4))
    push!(haplotype_set[4], (5, 6))
    push!(haplotype_set[4], (5, 8))
    push!(haplotype_set[4], (8, 4))
    push!(haplotype_set[4], (8, 6))
    push!(haplotype_set[4], (8, 8))

    sol_path, next_pair, subtree_err, best_err =
        MendelImpute.connect_happairs!(haplotype_set)

    @test best_err == 1.0
    @test sol_path == [(1, 5), (1, 5), (1, 5), (5, 4)]
    @test all(subtree_err[1] .== [2.0; 1.0; 2.0; 3.0; 2.0; 3.0; 5.0; 2.0; 3.0])
    @test all(subtree_err[2] .== [1.0; 2.0; 2.0; 2.0; 5.0; 5.0; 2.0; 5.0; 5.0])
    @test all(subtree_err[3] .== [4.0; 1.0; 4.0; 1.0])
    @test all(subtree_err[end] .== 0.0)
    @test all(next_pair[1] .== [1; 1; 1; 4; 1; 4; 1; 1; 7])
    @test all(next_pair[2] .== 2)
    @test all(next_pair[3] .== [1; 4; 1; 4])
    @test all(next_pair[4] .== 0)
end

@testset "compress" begin
    # perform compression
    reffile = "ref.excludeTarget.vcf.gz"
    tgtfile = "target.typedOnly.masked.vcf.gz"
    outfile = "ref.excludeTarget.jlso"
    @time compress_haplotypes(reffile, tgtfile, outfile)

    # read compressed data
    @time compressed_Hunique = JLSO.load(outfile)[:compressed_Hunique]

    # tests on summary statistics
    @test MendelImpute.nhaplotypes(compressed_Hunique) == 4800
    @test MendelImpute.nwindows(compressed_Hunique) == 16
    @test MendelImpute.max_dim(compressed_Hunique) == (313, 767)
    @test MendelImpute.get_window_widths(compressed_Hunique) == [313, 312,
        313, 312, 313, 312, 313, 312, 313, 312, 313, 312, 313, 312, 313, 312]
    @test MendelImpute.count_haplotypes_per_window(compressed_Hunique) == [736,
        752, 733, 742, 765, 650, 666, 745, 679, 695, 682, 767, 720, 665, 711,
        623]
    @test MendelImpute.avg_haplotypes_per_window(compressed_Hunique) == 708.1875

    # tests on conversion
    @test MendelImpute.unique_idx_to_complete_idx(100, 10, compressed_Hunique) == 127
    @test MendelImpute.unique_all_idx_to_complete_idx(100, 10, compressed_Hunique) == 112
    @test MendelImpute.unique_all_idx_to_unique_typed_idx(100, 10, compressed_Hunique) == 93
    @test MendelImpute.complete_idx_to_unique_all_idx(100, 10, compressed_Hunique) == 91
    @test MendelImpute.complete_idx_to_unique_typed_idx(100, 10, compressed_Hunique) == 85
end

@testset "impute & phasing" begin
    Random.seed!(2020)
    tgtfile = "target.typedOnly.masked.vcf.gz"
    reffile = "ref.excludeTarget.jlso"
    outfile = "imputed.vcf.gz"
    @time phase(tgtfile, reffile, outfile)

    # check error rate
    Ximputed = convert_gt(Float64, "imputed.vcf.gz")
    Xtrue = convert_gt(Float64, "target.full.vcf.gz")
    m, n = size(Xtrue)
    error_rate = sum(Xtrue .!= Ximputed) / m / n
    @test (m, n) == (100, 36063)
    @test error_rate ≈ 0.0007129190583146162

    # check per-sample error
    quality = CSV.read("imputed.sample.error", DataFrame)
    @test size(quality, 1) == 100
    @test count(isone, quality[!, :error]) == 65
    @test mean(quality[!, :error]) ≈ 0.9996911073712372
    @test std(quality[!, :error]) ≈ 0.0006256546642253605

    # check per-SNP error
    reader = VCF.Reader(openvcf(outfile, "r"))
    snpscores = Vector{Float64}(undef, nrecords(outfile))
    for (i, record) in enumerate(reader)
        snpscores[i] = parse(Float64, VCF.info(record)[1].second)
    end
    close(reader)
    @test length(snpscores) == 36063
    @test count(isone, snpscores) == 29934
    @test mean(snpscores) ≈ 0.9997004408951002
    @test std(snpscores) ≈ 0.0007984939339310984
end

@testset "PLINK and dosage inputs" begin
    # compute .vcf.gz (hard genotype) error rate
    tgtfile = "target.typedOnly.masked.vcf.gz"
    reffile = "ref.excludeTarget.jlso"
    outfile = "imputed.vcf.gz"
    @time phase(tgtfile, reffile, outfile);
    Ximputed = convert_gt(Float64, "imputed.vcf.gz")
    Xtrue = convert_gt(Float64, "target.full.vcf.gz")
    m, n = size(Xtrue)
    vcf_error = sum(Xtrue .!= Ximputed) / m / n

    # PLINK error rate
    tgtfile = "target.typedOnly.masked" # do no include bim/bed/fam in filename
    reffile = "ref.excludeTarget.jlso"
    outfile = "imputed.vcf.gz"
    @time phase(tgtfile, reffile, outfile)
    Ximputed = convert_gt(Float64, "imputed.vcf.gz")
    plink_error = sum(Xtrue .!= Ximputed) / m / n
    @test plink_error ≈ vcf_error

    # dosage error rate
    tgtfile = "target.typedOnly.dosages.masked.vcf.gz"
    reffile = "ref.excludeTarget.jlso"
    outfile = "imputed.vcf.gz"
    @time phase(tgtfile, reffile, outfile, dosage=true)
    Ximputed = convert_gt(Float64, "imputed.vcf.gz")
    dosage_error = sum(Xtrue .!= Ximputed) / m / n
    @test round(dosage_error, digits = 4) == round(vcf_error, digits=4)
end

@testset "decompress ultra-compressed files" begin
    # saving
    tgtfile = "target.typedOnly.masked.vcf.gz"
    reffile = "ref.excludeTarget.jlso"
    outfile = "imputed.jlso" # outfile ends in .jlso
    @time phase(tgtfile, reffile, outfile);

    # loading requires original reference panel in vcf format
    tgtfile = "imputed.jlso"
    reffile = "ref.excludeTarget.vcf.gz"
    X1, X2, phaseinfo, sampleID, H = convert_compressed(Float64,
        tgtfile, reffile);

    # test error rate
    Xtrue = convert_gt(Float64, "target.full.vcf.gz")
    Ximputed = (X1 + X2)' # transpose X1 and X2
    m, n = size(Ximputed)
    error_rate = sum(Xtrue .!= Ximputed) / m / n
    @test error_rate ≈ 0.0007129190583146162
end

@testset "0 allocations in phasing" begin
    reffile = "ref.excludeTarget.jlso"
    compressed_Hunique = JLSO.load(reffile)[:compressed_Hunique]

    # first person's optimal haplotype in each window (complete index)
    happair1_original = Int32[56, 624, 22, 20, 54, 221, 515, 172, 989, 19,
        171, 171, 570, 37, 570, 106]
    happair2_original = Int32[624, 989, 624, 624, 624, 515, 989, 989, 4115,
        171, 989, 202, 1681, 1681, 4115, 570]

    # preallocate vectors
    survivors1=Int32[]
    survivors2=Int32[]
    sizehint!(survivors1, 60000)
    sizehint!(survivors2, 60000)
    happair1 = copy(happair1_original)
    happair2 = copy(happair2_original)
    MendelImpute.phase_sample!(happair1, happair2,
        compressed_Hunique, survivors1, survivors2)

    # zero allocation and memory
    b = @benchmark MendelImpute.phase_sample!(happair1, happair2,
        $compressed_Hunique, $survivors1, $survivors2) setup=(happair1=
        copy($happair1_original);happair2=copy($happair2_original))

    @test b.allocs == 0
    @test b.memory == 0
end

@testset "0 allocation in haplopair" begin
    tgtfile = "target.typedOnly.masked.vcf.gz"
    reffile = "ref.excludeTarget.jlso"
    compressed_Hunique = JLSO.load(reffile)[:compressed_Hunique]

    # a bunch of setup code
    X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt = VCFTools.convert_gt(
        UInt8, tgtfile, trans=true, save_snp_info=true, msg =
        "Importing genotype file...");
    people = size(X, 2)
    tgt_snps = size(X, 1)
    ref_snps = length(compressed_Hunique.pos)
    threads = Threads.nthreads()
    max_width, max_d = MendelImpute.max_dim(compressed_Hunique)
    inv_sqrt_allele_var = nothing
    overlap = compressed_Hunique.overlap
    windows = MendelImpute.nwindows(compressed_Hunique)
    num_unique_haps = round(Int,
        MendelImpute.avg_haplotypes_per_window(compressed_Hunique))

    # preallocated arrays in phase
    ph = [MendelImpute.HaplotypeMosaicPair(ref_snps) for i in 1:people]
    haplotype1 = [zeros(Int32, windows) for i in 1:people]
    haplotype2 = [zeros(Int32, windows) for i in 1:people]

    # working arrys in happair
    happair1 = [ones(Int32, people)               for _ in 1:threads]
    happair2 = [ones(Int32, people)               for _ in 1:threads]
    hapscore = [zeros(Float32, people)            for _ in 1:threads]
    Xwork    = [zeros(Float32, max_width, people) for _ in 1:threads]
    Hwork    = [zeros(Float32, max_width, max_d)  for _ in 1:threads]

    # window 1
    w = 1
    Hw_aligned = compressed_Hunique.CW_typed[w].uniqueH
    Xrange = MendelImpute.extend_to_overlap_range(compressed_Hunique, w, overlap)
    Xw_aligned = view(X, Xrange, :)
    d = size(Hw_aligned, 2)
    id = Threads.threadid()

    # global search must allocate M, Hwork, N for each window
    # but there's no other allocation
    M = zeros(Float32, size(Hw_aligned, 2), size(Hw_aligned, 2))
    N = zeros(Float32, size(Xw_aligned, 2), size(Hw_aligned, 2))
    Hwork = convert(Matrix{Float32}, Hw_aligned)
    b = @benchmark MendelImpute.haplopair!($Xw_aligned, $Hw_aligned,
        happair1=$(happair1[id]), happair2=$(happair2[id]),
        hapscore=$(hapscore[id]), Xwork=$(Xwork[id]), M=$M, N=$N, Hwork=$Hwork)

    @test b.allocs == 0
    @test b.memory == 0
end

@testset "import BGEN" begin
    # compare BGEN genotypes to VCF genotypes
    @time G, Gsamples, Gchr, Gpos, GsnpID, Gref, Galt = MendelImpute.convert_gt(Float64,
        Bgen("target.typedOnly.masked.bgen"))
    @time Gtrue, sampleID, record_chr, record_pos, record_ids, record_ref,
        record_alt = convert_gt(Float64, "target.typedOnly.masked.vcf.gz",
        trans=true, save_snp_info=true)
    @test all(skipmissing(Gtrue .== G))
    @test all(Gsamples .== sampleID)
    @test all(Gchr .== record_chr)
    @test all(Gpos .== record_pos)
    @test all(GsnpID .== record_ids)
    @test all(Gref .== record_ref)
    @test all(Galt .== record_alt)

    # reference panel (BGEN vs VCF)
    @time H, Hsamples, Hchr, Hpos, HsnpID, Href, Halt =
        MendelImpute.convert_ht(Bgen("ref.excludeTarget.bgen"))
    @time Htrue, sampleID, record_chr, record_pos, record_ids, record_ref,
        record_alt = VCFTools.convert_ht(Bool, "ref.excludeTarget.vcf.gz", trans=true,
        save_snp_info=true)
    @test all(Htrue .== H)
    @test all(Hsamples .== sampleID)
    @test all(Hchr .== record_chr)
    @test all(Hpos .== record_pos)
    @test all(HsnpID .== record_ids)
    @test all(Href .== record_ref)
    @test all(Halt .== record_alt)
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using DataStructures, Distributions, CSV, DataFrames, Random, Parameters

## TO DO
# make move() failsafe
# input a threshold max number of cells
# increase compatibility with ODE model (e.g. polyharmonicSpline birth rates; use chemotaxis_form() fn's; & etc)
# modify birth_rate() to allow testing with different maxChrom
mutable struct karyotype
    state
    viable
end

function missegregate(k::karyotype, maxChrom,chromosome_selector)
    #select a chromosome randomly
    index = rand(chromosome_selector,1)[1]
    k_new=deepcopy(k)
    #generate daughters according to rule
    k.state[index]=k.state[index]+1
    k_new.state[index]=k_new.state[index]-1
    #check viability
    if k.state[index[1]]<1 || k.state[index[1]]>maxChrom
        k.viable=false
    end
    if k_new.state[index[1]]<1 || k_new.state[index[1]]>maxChrom
        k_new.viable=false
    end
    return k, k_new
end

## handy to set up some sort of forced trajectory of evolution for testing the model
## f(cn) = birthrate: this could easily be swapped with polyharmspline in future
function get_birthrate(cn)
    ## if [1,1,1] is an input state it should go through these states in order:
    if cn == [1,1,1]
        return 0.3
    end
    if cn == [1,1,2]
        return 0.3
    end

    if cn == [1,2,1]
        return 0.3
    end

    if cn == [2,1,1]
        return 0.3
    end

    if cn == [1,1,3]
        return 0.5
    end

    if cn == [1,2,3]
        return 0.6
    end

    if cn == [1,3,3]
        return 0.7
    end

    if cn == [2,3,3]
        return 0.8
    end

    return 0

end

#a cell object
mutable struct cell_allinfo
    x::Float32
    y::Float32
    state
    cell_allinfo(x,y,state) = new(x,y,state) # constructor
end

# a reduced cell object (cells often don't need to know their karyotype)
mutable struct cell
    x::Float32
    y::Float32
    cell(x,y) = new(x,y) # constructor
end

# migrate a cell
function move(c::cell,Np,dnorm, dE0,E0,χ,Ξ)
    ix = ceil(Int64,c.x)
    iy = ceil(Int64,c.y)
    new_x = c.x + rand(dnorm,1)[1] + χ*dE0[1,ix,iy]*log(E0[ix,iy]+Ξ)
    new_y = c.y + rand(dnorm,1)[1] + χ*dE0[2,ix,iy]*log(E0[ix,iy]+Ξ)

    ##implement boundary conditions (NOTE this can be broken if migration longer than Np])
    ## consider % operator to fix overstepping? or just do..while
    if new_x > Np[1]
        new_x=new_x-Np[1]
    end
    if new_y > Np[1]
        new_y=new_y-Np[1]
    end
    if new_x < 0
        new_x=new_x+Np[1]
    end
    if new_y < 0
        new_y=new_y+Np[1]
    end
    c.x = new_x
    c.y = new_y
    return c
end

#store all the cells for a given clone
mutable struct clone
    pop
    cn
    birthrate
end

# wrapper function for move()
function migrate(cl::clone,Np,dnorm,dE0,E0,χ,Ξ)
    cl.pop = map(x->move(x,Np,dnorm,dE0,E0,χ,Ξ),cl.pop)
    return(cl)
end

# loop thru all cells per clone, and perform divisions
function replicate(cl::clone, E0, ϕ, dt, uniform01,misrate,maxChrom,chromosome_selector,Np,dnorm,dE0,χ,Ξ)
    # for missegregations we will take a strategy of returning all the new viable karyotypes
    # and then allowing the stochasticCompartment object to deal with them
    missegs = MutableLinkedList{cell_allinfo}()
    for i in length(cl.pop):-1:1
        energy = E0[ceil(Int64,cl.pop[i].x),ceil(Int64,cl.pop[i].y)]
        #need to replace with call to energy function
        probBirth = (cl.birthrate-misrate)*dt*energy/(ϕ + energy)
        if rand(uniform01,1)[1] < probBirth  # division occurs
            if rand(uniform01,1)[1] < misrate #missegregation occurs

                k=karyotype(deepcopy(cl.cn),true)
                k1,k2 = missegregate(k,maxChrom,chromosome_selector)
                if k1.viable
                    c1=cell_allinfo(cl.pop[i].x,cl.pop[i].y,k1.state)
                    push!(missegs,c1)
                end
                if k2.viable
                    c2=cell_allinfo(cl.pop[i].x,cl.pop[i].y,k2.state)
                    push!(missegs,c2)
                end
                delete!(cl.pop,i)
            else
                new_c = cell(cl.pop[i].x,cl.pop[i].y)
                new_c=move(new_c,Np,dnorm,dE0,E0,χ,Ξ)
                push!(cl.pop,new_c)
            end
        end
    end
    return cl,missegs
end

# test if all cells per clone need to die, and kill them
function die(cl::clone,dt, deathRate, uniform01)
    # require iterating backwards when deleting or else
    # we will overrun the list

    for i in length(cl.pop):-1:1
        probDeath = deathRate*dt
        if rand(uniform01,1)[1] < probDeath  # division occurs
            delete!(cl.pop,i)
        end
    end
    return(cl)
end

# top level data structure for stochastic compartment, storing state of all clones
# and any other useful variables etc.
mutable struct stochasticCompartment
    popDict # structure to hold all clones (dictionary)
    #store parms used in stochastic model for convenience:
    Np
    ϕ
    χ
    Ξ
    maxChrom
    dt
    misrate
    deathRate
    # store the Distributions used in stochastic model:
    chromosome_selector
    uniform01
    dnorm
end

# perform migration for entire stochastic model
function migrate(s::stochasticCompartment,dE0,E0)
    for k in keys(s.popDict)
        s.popDict[k] = migrate(s.popDict[k],s.Np,s.dnorm,dE0,E0,s.χ,s.Ξ)
    end
    return(s)
end
# perform cell division for entire stochastic model
function replicate(s::stochasticCompartment, E0,dE0,sIndex)
    for k in keys(s.popDict)
        s.popDict[k],missegs = replicate(s.popDict[k],E0, s.ϕ, s.dt, s.uniform01,
        s.misrate,s.maxChrom,s.chromosome_selector,s.Np,s.dnorm,dE0,s.χ,s.Ξ)
        for m in missegs
            if (m.state in sIndex)==false # check if already in det. comp.
                c = cell(m.x,m.y)
                c=move(c,s.Np,s.dnorm,dE0,E0,s.χ,s.Ξ)
                if (m.state in keys(s.popDict)) == false
                    pop = MutableLinkedList{cell}()
                    push!(pop, c)
                    br=get_birthrate(m.state)
                    ci = clone(pop,m.state,br)
                    s.popDict[m.state]=ci
                else
                    push!(s.popDict[m.state].pop,c)
                end
            # else... ASSUME missegregations from stochastic compartment into det. compartment
            # are negligible (!?)
            end
        end
        if length(s.popDict[k].pop) < 1
            delete!(s.popDict,k)
        end
    end
    return(s)
end
# perform cell death for entire stochastic model
function die(s::stochasticCompartment)
    for k in keys(s.popDict)
        s.popDict[k] = die(s.popDict[k],s.dt,s.deathRate,s.uniform01)
        if length(s.popDict[k].pop) < 1
            delete!(s.popDict,k)
        end
    end
    return(s)
end

# return grid with locations of all cells in stochastic compartment
function spatial_summary(s::stochasticCompartment)
    u=zeros(s.Np...)
    for cl in values(s.popDict)
        for c in cl.pop
            u[ceil(Int64,c.x),ceil(Int64,c.x)] = u[ceil(Int64,c.x),ceil(Int64,c.x)]+1
        end
    end
    return u
end

# introduce new cells into stochastic model from deterministic compartment
function getDaughters(stoch::stochasticCompartment, E0, s0, sIndex,birthRates)
    for i in 1:(size(s0)[1])
        br = birthRates[i]
        for j in 1:(size(s0)[2])
            for k in 1:(size(s0)[3])
                pBirth = br*stoch.dt*s0[i,j,k]*E0[j,k]/(stoch.ϕ + E0[j,k]) # expected number of births per site
                pMis = pBirth*stoch.misrate # expected number of missegregations per site
                p = Poisson(pMis)
                Np=rand(p,1)[1]
                if Np>0
                    for n in 1:Np
                        k0=karyotype(deepcopy(sIndex[i]),true)
                        k1,k2=missegregate(k0, stoch.maxChrom, stoch.chromosome_selector)
                        for ki in [k1,k2]
                            if ki.viable
                                # create a new cell
                                c=cell((j-rand()),(k-rand()))
                                # check if cell karyotype exists in clones
                                if (ki.state in keys(stoch.popDict)) == false
                                #popDict[k.state]=1
                                    pop = MutableLinkedList{cell}()
                                    push!(pop, c)
                                    br=get_birthrate(ki.state) # will look this up later
                                    ci = clone(pop,ki.state,br)
                                    stoch.popDict[ki.state]=ci
                                else
                                    push!(stoch.popDict[ki.state].pop,c)
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return(stoch)
end

# get x,y & cn of each cell and return in a dataframe
function summarise_clone(ci::clone)
    df = repeat(ci.cn',outer=(length(ci.pop),1))
    xy = hcat(map(c->[c.x,c.y],ci.pop)...)'
    return DataFrame(hcat(df,xy),:auto)
end

## should take as arguments the energy array, the cell array, and stochasticCOmpartmentobject
## then return updated versions of each
function stochastic_step(s::stochasticCompartment, E0, s0, info)
    @unpack (sIndex,birthRates) = info
    dE0 = grad(E0,s.Np)
    #unclear what is best order to do these steps
    s=die(s)
    s=migrate(s,dE0,E0)
    s=replicate(s,E0,dE0,sIndex)
    s=getDaughters(s, E0, s0, sIndex,birthRates)
    u=spatial_summary(s)
    return s, u
end

## wrapper function to setup the stochastic model
## pass to the functions variables that remain constant
## throughout simulation etc. Should return a stochasticCompartment object
function setup_stochastic(params)

    @unpack (Np,ϕ, nChrom,maxChrom,dt,misrate,deathRate,Γ,χ,Ξ) = params
    chromosome_selector = DiscreteUniform(1,nChrom) # random integer 1:nchrom
    uniform01 = Uniform(0,1) # random continous uniform dist.
    dnorm = Normal(0,dt*Γ) #normal distribution for cell migration

    popDict = Dict{Array,clone}()
    s = stochasticCompartment(popDict,Np,ϕ,χ,Ξ, maxChrom,dt,misrate,deathRate,
    chromosome_selector,uniform01,dnorm)
    return(s)
end

# implementation of wrapped boundary coords
function wrap(index,Npi)
    if index>Npi
        return index-Npi
    end
    if index<1
        return index+Npi
    end
    return index
end

#get the gradient of a 2d array
function grad(R2,Np)
    dR2 = zeros(2,size(R2)...)
    for i in 1:size(R2)[1]
        for j in 1:size(R2)[2]
            #1st order central diff
            dR2[1,i,j]=(R2[wrap(i+1,Np[1]),j]-R2[wrap(i-1,Np[1]),j])/(2*Np[1])
            dR2[2,i,j]=(R2[i,wrap(j+1,Np[2])]-R2[i,wrap(j-1,Np[2])])/(2*Np[2])
        end
    end
    return dR2
end

max_stars_repo_name: MathOnco/ploidyEvolution

max_issues_repo_name: MathOnco/ploidyEvolution

max_forks_repo_name: MathOnco/ploidyEvolution


text:
function SharedVertices(convexexp1in2::AbstractArray{Float64, 2},
                        ordered_vertices1::Vector{Int},
                        ordered_vertices2::Vector{Int},
                        numof1in2::Int,
                        numof2in1::Int)

    reordered_vertices1 = copy(ordered_vertices1)
    reordered_vertices2 = copy(ordered_vertices2)
    Ncomm = 0
    M = numof1in2 * numof2in1

    if M > 0 # The simplices might share some vertex
        contained1in2 = collect(1:numof1in2)
        contained2in1 = collect(1:numof2in1)

        referenceexpansion = Matrix(1.0I, size(convexexp1in2, 1), size(convexexp1in2, 2))

        exp1 = convexexp1in2[:, reordered_vertices1[1:numof1in2]]
        exp2 = referenceexpansion[:, reordered_vertices2[1:numof2in1]]

        num1 = numof1in2
        num2 = numof2in1

        # Create matrices to compare the expansions
        exp1Xones_num2 = kron(exp1, ones(1, num2))
        ones_num1Xexp2 = kron(ones(1, num1), exp2)

        temp = maximum(abs.(exp1Xones_num2 - ones_num1Xexp2), dims = 1)
        x = transpose(heaviside0(-temp)) .* collect(1:M)
        coincidences = (LinearIndices(x))[findall(x->x!=0, x)]
        coincidences = coincidences
        Ncomm = length(coincidences)

        if Ncomm > 0
            # Indices in 1:numof1in2 corresponding to the vertices in simplex1 that are shared
            internal_indices_sharedvert_in_1 = ceil.(Int64, coincidences ./ num2)
            # Indices in 1:numof2in1 corresponding to the vertices in simplex2 that are shared
            internal_indices_sharedvert_in_2 = round.(Int64, coincidences .- num2 * (internal_indices_sharedvert_in_1 .- 1))

            # Common vertices in 1
            common_indices_in1 = contained1in2[internal_indices_sharedvert_in_1]
            Temp1 = reordered_vertices1[common_indices_in1]
            if Ncomm < numof1in2
                indices_not_common_in1 = complementary(common_indices_in1, numof1in2)
                Temp1 = [reordered_vertices1[common_indices_in1]; reordered_vertices1[indices_not_common_in1]]
            end
            reordered_vertices1[1:numof1in2] = Temp1

            # Common vertices in 2
            common_indices_in2 = contained2in1[internal_indices_sharedvert_in_2]
            Temp2 = reordered_vertices2[common_indices_in2]
            if Ncomm < numof2in1
                indices_not_common_in2 = complementary(common_indices_in2, numof2in1)
                Temp2 = [reordered_vertices2[common_indices_in2]; reordered_vertices2[indices_not_common_in2]]
            end
            reordered_vertices2[1:numof2in1] = Temp2

        end
    end

    return Ncomm, reordered_vertices1, reordered_vertices2
end


text:
module SpringCollab2020FlagTouchSwitch

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/FlagTouchSwitch" FlagTouchSwitch(x::Integer, y::Integer,
    flag::String="flag_touch_switch", icon::String="vanilla", persistent::Bool=false, inactiveColor::String="5FCDE4", activeColor::String="FFFFFF", finishColor::String="F141DF")

const bundledIcons = String["vanilla", "tall", "triangle", "circle"]

const placements = Ahorn.PlacementDict(
    "Flag Touch Switch (Spring Collab 2020)" => Ahorn.EntityPlacement(
        FlagTouchSwitch
    )
)

Ahorn.editingOptions(entity::FlagTouchSwitch) = Dict{String,Any}(
    "icon" => bundledIcons
)

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

    return  Ahorn.Rectangle(x - 7, y - 7, 14, 14)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagTouchSwitch, room::Maple.Room)
    Ahorn.drawSprite(ctx, "objects/touchswitch/container.png", 0, 0)

    icon = get(entity.data, "icon", "vanilla")

    iconPath = "objects/touchswitch/icon00.png"
    if icon != "vanilla"
        iconPath = "objects/SpringCollab2020/flagTouchSwitch/$(icon)/icon00.png"
    end

    Ahorn.drawSprite(ctx, iconPath, 0, 0)
end

end


text:
type Raster
    dataset::GDALDatasetH
    width::Cint
    height::Cint
    nband::Cint # number of raster bands
end

function Base.show(io::IO, raster::Raster)
    if raster.dataset == C_NULL
        print(io, "Closed Raster Dataset")
    else
        dptr = _getdatasetdriver(raster.dataset)
        driverstring = "$(drivershortname(dptr))/$(driverlongname(dptr))"
        println(io, "RasterIO.Raster (Driver: $driverstring)")
        print(io, "\nFile(s): ")
        for filename in getfilelist(raster.dataset)
            print(io, "$filename ")
        end
        print(io, "\nDataset (width x height): ")
        println(io, "$(raster.width) x $(raster.height) (pixels)")
        for i = 1:min(raster.nband, 5)
            band = getrasterband(raster.dataset, i)
            color = colorinterpname(getcolorinterp(band))
            xsize = _getrasterbandxsize(band)
            ysize = _getrasterbandysize(band)
            bandtype = dtype(_getrasterband(raster.dataset, i))
            println(io, "  Band $i ($color): $xsize x $ysize ($bandtype)")
        end
        raster.nband > 5 && println(io, "  ...")
        print(io, "Number of bands: $(raster.nband)")
    end
end

function Raster(dataset::GDALDatasetH)
    raster = Raster(dataset,
                    _getrasterxsize(dataset), # width
                    _getrasterysize(dataset), # height
                    _getrastercount(dataset)) # number of bands
    finalizer(raster, closeraster)
    raster
end

function closeraster(raster::Raster)
    if raster.dataset != C_NULL
        _close(raster.dataset)
        raster.dataset = C_NULL
    end
end


text:
Ns() = [1_000,]
# Ns() = [10,]
Ds() = [1,]

generate_x̄s(x, x̄s::Nothing) = [fill(x, N) for N in Ns()]
generate_x̄s(x, x̄s) = x̄s

function create_benchmarks(name, grads::Bool, f, x...)
    if grads
        @benchset name begin
            @bench "eval" $f($x...)
            __forward() = Zygote.pullback(f, x...)
            out, back = Zygote.pullback(f, x...)
            @bench "forward" $__forward()
            @bench "back" $back($out)
            # @bench "gradient" $__gradient()
        end
    else
        @bench name $f($x...)
    end
    return nothing
end

# Common benchmarks to run for mean functions.
function create_benchmarks(μ::MeanFunction; grads=true, x=5.0, x̄s=nothing)
    create_benchmarks("μ(x)", grads, μ, x)

    for x̄ in generate_x̄s(x, x̄s)
        N = length(x̄)
        create_benchmarks("map ($(length(x̄)))", grads, x->map(μ, x), x̄)
    end
end

# Common benchmarks for CrossKernel.
function create_benchmarks(k::CrossKernel; x=5.0, x′=4.0, x̄s=nothing, x̄′s=nothing, grads=true)
    create_benchmarks("k(x, x′)", grads, k, x, x′)

    for (x̄, x̄′) in zip(generate_x̄s(x, x̄s), generate_x̄s(x′, x̄′s))
        N = length(x̄)
        @benchset "$N" begin
            create_benchmarks("map k(x̄, x̄′) $N", grads, (x̄, x̄′)->map(k, x̄, x̄′), x̄, x̄′)
            create_benchmarks("pw k(x̄, x̄′) $N", grads, (x̄, x̄′)->pairwise(k, x̄, x̄′), x̄, x̄′)
        end
    end
end

# Common benchmarks to run for Kernels.
function create_benchmarks(k::Kernel; x=5.0, x′=4.0, x̄s=nothing, x̄′s=nothing, grads=true)
    create_benchmarks("k(x)", grads, k, x)
    create_benchmarks("k(x, x′)", grads, k, x, x′)

    for (x̄, x̄′) in zip(generate_x̄s(x, x̄s), generate_x̄s(x′, x̄′s))
        N = length(x̄)
        @benchset "$N" begin
            create_benchmarks("map k(x̄) $N", grads, x̄->map(k, x̄), x̄)
            create_benchmarks("map k(x̄, x̄′) $N", grads, (x̄, x̄′)->map(k, x̄, x̄′), x̄, x̄′)
            create_benchmarks("pw k(x̄) $N", grads, x̄->pairwise(k, x̄), x̄)
            create_benchmarks("pw k(x̄, x̄′) $N", grads, (x̄, x̄′)->pairwise(k, x̄, x̄′), x̄, x̄′)
        end
    end
end

function pretty_print(d::BenchmarkGroup, pre=1)
    pretty_print(d.data, pre)
end

function pretty_print(d::Dict, pre=1)
    for (k, v) in d
        if v isa BenchmarkGroup
            s = "$(repr(k))"
            println(join(fill(" ", pre)) * s)
            pretty_print(v, pre+1+4)
        else
            println(join(fill(" ", pre)) * "$(repr(k)) => $(repr(v))")
        end
    end
    nothing
end


text:
module BenefitCap
#
# Apply the benefit cap as described in CPAG 19/20
# ch 52. Somewhat impressionistic and could do
# with another go.
#
using ScottishTaxBenefitModel

using .ModelHousehold:
    BenefitUnit,
    Household,
    Person

using .STBParameters:
    BenefitCapSys

using .Results:
    BenefitUnitResult

using .Definitions

using .Intermediate:
    MTIntermediate

using .STBIncomes

export apply_benefit_cap!

"""
Apply a benefit cap to a benefit unit.
"""
function apply_benefit_cap!(
    benefit_unit_result :: BenefitUnitResult,
    region           :: Standard_Region,
    benefit_unit     :: BenefitUnit,
    intermed         :: MTIntermediate,
    caps             :: BenefitCapSys,
    route            :: LegacyOrUC )
    # println( "apply_benefit_cap entered; route = $route")
    bu = benefit_unit # shortcut
    bur = benefit_unit_result # shortcut

    if route == legacy_bens
        for pid in bu.adults
            if bur.pers[pid].income[WORKING_TAX_CREDIT] > 0
                # cpag 21/2 p 1182
                # println("wtc bailing out")
                return
            end
        end
    else
        gross_earnings = 0.0
        for pid in bu.adults
            gross_earnings += isum( bur.pers[pid].income, [WAGES,SELF_EMPLOYMENT_INCOME] )
        end
        if gross_earnings >= caps.uc_incomes_limit
            # println("gross earn; bailing out ")
            return
        end
    end

    if intermed.someone_pension_age ||
        intermed.someone_is_carer ||
        (intermed.num_severely_disabled_adults > 0) ||
        caps.abolished
        # println("pension age bailing out")
        return
    end

    cap = intermed.num_people == 1 ?
        caps.outside_london_single :
        caps.outside_london_couple
    if region == London
        cap = intermed.num_people == 1 ?
            caps.inside_london_single :
            caps.inside_london_couple
    end
    # println("got cap as $cap")
    totbens = 0.0
    included = UC_CAP_BENEFITS
    target_ben = UNIVERSAL_CREDIT
    min_amount = bur.uc.childcare_costs
    if route == legacy_bens
        included = LEGACY_CAP_BENEFITS
        target_ben = HOUSING_BENEFIT
        min_amount = 0.5
    end
    recip_pers :: BigInt = -1
    recip_ben = 0.0

    for pid in bu.adults
        totbens += isum( bur.pers[pid].income, included )
        if bur.pers[pid].income[target_ben] > 0
            recip_pers = pid
            recip_ben = bur.pers[pid].income[target_ben]
        end
    end
    if recip_ben == 0.0
        # println("uc/hb0; returning ")
        return
    end
    excess = totbens - cap
    # println("totbens=$totbens cap=$cap excess=$excess")
    if excess > min_amount
        rd = max( min_amount, recip_ben - excess )
        bur.bencap.reduction = recip_ben - rd
        bur.pers[recip_pers].income[target_ben] = rd
    end
    bur.bencap.cap_benefits = totbens
    bur.bencap.cap = cap
end # cap_benefits

end # module


text:
type SubtreeCrossover <: Crossover
  stage::AbstractString
  rate::Float
  max_tries::Int
  max_depth::Int
  rep::KozaTreeRepresentation

  SubtreeCrossover(stage::AbstractString, rate::Float, max_tries::Int, max_depth::Int, rep::KozaTreeRepresentation) =
    new(stage, rate, max_tries, max_depth, rep)
end

"""
Provides a definition for a subtree crossover operator.
"""
type SubtreeCrossoverDefinition <: CrossoverDefinition
  """
  The name of the stage that this operator works on.
  """
  stage::AbstractString

  """
  The crossover rate for this operator, i.e. the probability that a given pair
  of trees will be subject to crossover.
  """
  rate::Float

  """
  The maximum number of tries that the operator when attempting to make a
  viable crossover before giving up.
  """
  max_tries::Int

  SubtreeCrossoverDefinition(stage::AbstractString, r::Float, mt::Int) =
    new(stage, r, mt)
end

"""
Composes a given sub-tree crossover operator from its definition.
"""
function compose!(c::SubtreeCrossoverDefinition, sp::Species)
  c.stage = c.stage == "" ? genotype(sp).label : c.stage
  r = sp.stages[c.stage].representation
  SubtreeCrossover(c.stage, c.rate, c.max_tries, r.max_depth, r)
end

"""
TODO: Write koza.subtree_crossover documentation.
"""
subtree_crossover() = subtree_crossover("", 0.7, 1)
subtree_crossover(stage::AbstractString) =
  subtree_crssover(stage, 0.7, 1)
subtree_crossover(rate::Float) = subtree_crossover("", rate, 1)
subtree_crossover(limit::Int) = subtree_crossover("", 0.7, limit)
subtree_crossover(stage::AbstractString, rate::Float) =
  subtree_crossover(stage, rate, 1)
subtree_crossover(stage::AbstractString, limit::Int) =
  subtree_crossover(stage, 0.7, limit)
subtree_crossover(rate::Float, limit::Int) = subtree_crossover("", rate, limit)
subtree_crossover(stage::AbstractString, rate::Float, limit::Int) =
  SubtreeCrossoverDefinition(stage, rate, limit)
function subtree_crossover(def::Function)
  c = subtree_crossover()
  def(c)
  c
end

num_inputs(o::SubtreeCrossover) = 2
num_outputs(o::SubtreeCrossover) = 2

"""
Performs a sub-tree crossover on a pair of provided Koza trees (belonging to
the same species).
"""
function operate!{T <: KozaTree}(
  o::SubtreeCrossover,
  outputs::Vector{IndividualStage{T}},
  x::T,
  y::T
)
  # Enforce the crossover rate
  rand() > o.rate && return

  # Choose a node at random from the first tree.
  nx = sample(x)

  # Continue to sample nodes from the 2nd tree until either a viable
  # crossover candidate has been found, or till a maximum number of
  # attempts have passed.
  for i in 1:o.max_tries
    ny = sample(y)
    if  (depth(ny) + height(nx) - 1) <= o.max_depth &&
        (depth(nx) + height(ny) - 1) <= o.max_depth
      swap!(ny, nx)
      break
    end
  end
end


text:
using ITensors,
      Random

function makeRandomMPS(sites;
                       chi::Int=4)::MPS
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(chi, "Link,l=$n") for n=1:N-1]
  for n=1:N
    s = sites[n]
    if n == 1
      v[n] = ITensor(l[n], s)
    elseif n == N
      v[n] = ITensor(l[n-1], s)
    else
      v[n] = ITensor(l[n-1], l[n], s)
    end
    randn!(v[n])
    normalize!(v[n])
  end
  return MPS(N,v,0,N+1)
end

function makeRandomMPO(sites;
                       chi::Int=4)::MPO
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(chi, "Link,l=$n") for n=1:N-1]
  for n=1:N
    s = sites[n]
    if n == 1
      v[n] = ITensor(l[n],s,s')
    elseif n == N
      v[n] = ITensor(l[n-1],s,s')
    else
      v[n] = ITensor(l[n-1],s,s',l[n])
    end
    randn!(v[n])
    normalize!(v[n])
  end
  return MPO(N,v,0,N+1)
end


text:
using Docile
make_term(term, typ, parent) =
    [:(immutable $typ <: $parent end),
     :(const $(esc(term))  = $typ())]

@doc """
`@terms` allows you to create singleton types,
create instances and put them in constants.

E.g.

    abstract A
    @terms A begin
        x => X
        y => Y
        z => Z
    end

results in

    abstract A
    immutable X <: A end
    const x = X()

    immutable Y <: A end
    const y = Y()

    immutable Z <: A end
    const z = Z()

""" ->
macro terms(parent, terms)
    args = filter(x -> x.head != :line, terms.args)
    Expr(:block,
        reduce(vcat, [make_term(arg.args[1], arg.args[2], parent)
            for arg in args])...)
end

function argskwargs(exps)
    kwargs = filter(exp -> exp.args[1] in [:kwarg, :typedkwarg], exps)
    args = filter(exp -> !(exp.args[1] in [:kwarg, :typedkwarg]), exps)
    args, kwargs
end
function typebody(exp::Expr)
    if exp.head == :call
        if exp.args[2].head == :(::)
            return exp.args[2]
        elseif exp.args[2].head == :(kw)
            return exp.args[2].args[1]
        end
    end
    error("Invalid API definition")
end

function typebody(args::AbstractArray)
    map(typebody, args)
end

hasdefault(exp) = exp.args[2].head == :kw
argtype(exp) = exp.args[1]
function getvar(exp)
    decl = exp.args[2]
    hasdefault(exp) ?
        (argtype(exp) === :typedarg ? decl.args[1] : decl.args[1].args[1]) :
        (argtype(exp) === :typedarg ? decl : decl.args[1])
end
striptype(exp::Symbol) = exp
striptype(exp) = exp.head == :(::) ? exp.args[1] : exp

replaceparam(x::Symbol, params) = x
function replaceparam(x, params)
    if x.head === :(::)
        typ = get(params, x.args[2], x.args[2])
        return Expr(x.head, x.args[1], esc(typ))
    else
        return Expr(x.head, x.args[1], esc(typ))
    end
end

function states(arg, params)
    var = replaceparam(getvar(arg), params)

    if hasdefault(arg)
        val = arg.args[2].args[2]

        return Any[(false, var, striptype(var)), (false, nothing, esc(val))]
    else
        if argtype(arg) === :curry
            return Any[(false, striptype(var), striptype(var)),
                       (true, striptype(var), striptype(var))]
        else
            return Any[(false, var, striptype(var))]
        end
    end
end

immutable ApiMethod
    args
    lambda_args
    constructor_args
end

function prependarg(state, method)
    curried, var, val = state
    lambda_args = method.lambda_args
    args = method.args
    constructor_args = method.constructor_args

    if curried
        lambda_args = vcat(var, lambda_args)
    else
        if var != nothing
            args = vcat(var, args)
        end
    end
    constructor_args = vcat(val, constructor_args)
    ApiMethod(args, lambda_args, constructor_args)
end

function teeprint(x, fn=println)
    fn(x)
    x
end

function makeapimethods(arglist, params)
    if length(arglist) == 0
        return Any[ApiMethod(Any[], Any[], Any[])]
    else
        this = arglist[1]
        rest = arglist[2:end]
        return [prependarg(state, method)
            for state in states(this, params), method in makeapimethods(rest, params)]
    end
end

function kwize(argdef)
    @assert argdef.args[1] in [:typedkwarg, :kwarg]

    if argdef.args[1] === :typedkwarg
        return argdef.args[2]
    else
        var = striptype(argdef.args[2].args[1])
        return Expr(:kw, var, argdef.args[2].args[2])
    end
end

function methodexpr(fn, typ, kws, kwnames, m)
    m = isempty(m.lambda_args) ?
        :(foo($(m.args...)) = $(esc(typ))($(m.constructor_args...), $(kwnames...))) :
        :(foo($(m.args...)) = ($(m.lambda_args...)) -> $(esc(typ))($(m.constructor_args...), $(kwnames...)))

    if !isempty(kws)
        m.args[1].args = vcat(fn, Expr(:parameters, kws...), m.args[1].args[2:end])
    end
    m.args[1].args[1] = esc(fn)
    m
end

typexpr(typ::Symbol) = typ
typexpr(typ) =
    typ.head === :comparison ?
        Expr(:(<:), typ.args[1], typ.args[3]) : typ

typename(typ::Symbol) = typ
typename(typ) =
    typ.head in [:comparison, :curly] ?
        typename(typ.args[1]) : typ

paramdict(typ::Symbol) = Dict()
paramdict(typ) =
    typ.head === :curly ?
    [param.args[1] => param.args[2] for param in typ.args[2:end]] :
        (typ.head === :(<:) ?
            paramdict(typ.args[1]) : Dict())
@doc """
`@api` is used to create new tile types and associated constructors

    @api border => Bordered{T <: Side} <: Tile begin
        arg(side::T)
        curry(tile::Tile)
        kwarg(color::ColorValue=color("black"))
    end

 Should result in:

     immutable Bordered{T <: Side} <: Tile
         tile::Tile
         color::ColorValue
     end
     border(side, tiles; color::ColorValue=color("black")) = Bordered(side, tiles, color)
     border(side; kwargs...) = tiles -> border(tiles; kwargs...)

""" ->
macro api(names, body)

    fn, typ = names.args
    fields = body.args

    argdefs = filter(f -> f.head != :line, fields)

    typedef = Expr(:type, false, esc(typexpr(typ)),
        Expr(:block, typebody(argdefs)...))

    dict = paramdict(typexpr(typ))

    args, kwargs = argskwargs(argdefs)
    methoddefs = makeapimethods(args, dict)
    kws = map(kwize, kwargs)
    kwnames = map(x -> striptype(x.args[1]), kws)
    ms = map(m -> methodexpr(fn, typename(typ), kws, kwnames, m), methoddefs)

    Expr(:block, typedef, ms...)
end


text:
include("core.jl")

function load_espresso()
    for n in Base.names(Espresso, true) @eval import Espresso: $n end
    for n in Base.names(XDiff, true) @eval import XDiff: $n end
end

logistic(x) = 1 ./ (1 + exp.(-x))
@scalardiff logistic(x::Number) 1 (logistic(x) .* (1 .- logistic(x)))

function find_bad(g)
    for i=1:length(g.tape)
        println("Evaliating $(i)th node $(g[i])")
        evaluate!(g, g[i])
    end
end

function predict(W, b, x)
    return W * x .+ b
end

function predict_cost(W, b, x, y)
    sum(predict(W, x, b) .- y)
end

function main_1642()
    W = rand(30, 64); b = rand(30); x = rand(64, 10)
    inputs = [:W => W, :b => b, :x => x]
    ctx = Dict(:codegen => EinCodeGen())
    ex = :(W * x .+ b)

    dex = xdiff(ex; ctx=ctx, inputs...)
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    NoProximal

Empty functor object for running a quasi-gradient algorithm without a prox step.

"""
struct NoProximal <: AbstractProximal end

function initialize_prox!(::AbstractQuasiGradient, ::NoProximal)
    return nothing
end

function restore_proximal_master!(::AbstractQuasiGradient, ::NoProximal)
    return nothing
end

function prox!(::AbstractQuasiGradient, ::NoProximal, x::AbstractVector, ∇f::AbstractVector, γ::AbstractFloat)
    x .= x - γ*∇f
    return nothing
end

# API
# ------------------------------------------------------------
"""
    NoProx

Factory object for [`NoProximal`](@ref). Passed by default to `prox` in `QuasiGradient.Optimizer`.

"""
struct NoProx <: AbstractProx end

function (::NoProx)(::StageDecompositionStructure, ::AbstractVector, ::Type{T}) where T <: AbstractFloat
    return NoProximal()
end

function str(::NoProx)
    return ""
end


text:
"""
    test_modification_set_function_single_variable(
        model::MOI.ModelLike,
        config::Config,
    )

Test that modifying the function of a `VariableIndex`-in-`LessThan` constraint
throws a [`SettingVariableIndexNotAllowed`](@ref) error.
"""
function test_modification_set_function_single_variable(
    model::MOI.ModelLike,
    ::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x, y
    maxobjective: 1.0x + 1.0y
    x <= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    c = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(x.value)
    # We test this after the creation of every `VariableIndex` constraint
    # to ensure a good coverage of corner cases.
    @test c.value == x.value
    err = MOI.SettingVariableIndexNotAllowed()
    @test_throws err MOI.set(model, MOI.ConstraintFunction(), c, y)
    return
end

"""
    test_modification_set_singlevariable_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test set modification VariableIndex-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_set_singlevariable_lessthan(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    maxobjective: 1.0x
    x <= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(x.value)
    @test c.value == x.value
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, -1.0)],
    )
    MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(2.0))
    @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.LessThan(2.0)
    _test_model_solution(
        model,
        config;
        objective_value = 2.0,
        variable_primal = [(x, 2.0)],
        constraint_primal = [(c, 2.0)],
        constraint_dual = [(c, -1.0)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_set_singlevariable_lessthan),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [2.0]),
            MOI.FEASIBLE_POINT,
        ),
    )
    return
end

"""
    test_modification_transform_singlevariable_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test set transformation of a VariableIndex-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_transform_singlevariable_lessthan(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    maxobjective: 1.0x
    x <= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(x.value)
    @test c.value == x.value
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, -1.0)],
    )
    c2 = MOI.transform(model, c, MOI.GreaterThan(2.0))
    @test !MOI.is_valid(model, c)
    @test MOI.is_valid(model, c2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = 2.0,
        variable_primal = [(x, 2.0)],
        constraint_primal = [(c2, 2.0)],
        constraint_dual = [(c2, 1.0)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_transform_singlevariable_lessthan),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [2.0]),
            MOI.FEASIBLE_POINT,
        ),
    )
    return
end

"""
    test_modification_set_scalaraffine_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying set of ScalarAffineFunction-in-LessThan constraint. If
`config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_set_scalaraffine_lessthan(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    maxobjective: 1.0x
    c: 1.0x <= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, -1.0)],
    )
    MOI.set(model, MOI.ConstraintSet(), c, MOI.LessThan(2.0))
    @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.LessThan(2.0)
    _test_model_solution(
        model,
        config;
        objective_value = 2.0,
        variable_primal = [(x, 2.0)],
        constraint_primal = [(c, 2.0)],
        constraint_dual = [(c, -1.0)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_set_scalaraffine_lessthan),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [2.0]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
    )
    return
end

"""
    test_modification_coef_scalaraffine_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying a variable coefficient in a ScalarAffineFunction-in-LessThan
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_coef_scalaraffine_lessthan(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    maxobjective: 1.0x
    c: 1.0x <= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, -1.0)],
    )
    MOI.modify(model, c, MOI.ScalarCoefficientChange(x, 2.0))
    _test_model_solution(
        model,
        config;
        objective_value = 0.5,
        variable_primal = [(x, 0.5)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, -0.5)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_coef_scalaraffine_lessthan),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.5]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-0.5],
        ),
    )
    return
end

"""
    test_modification_func_scalaraffine_lessthan(
        model::MOI.ModelLike,
        config::Config,
    )

Test setting the function in a ScalarAffineFunction-in-LessThan
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_func_scalaraffine_lessthan(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    maxobjective: 1.0x
    c: 1.0x <= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, -1.0)],
    )
    MOI.set(
        model,
        MOI.ConstraintFunction(),
        c,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 0.0),
    )
    foo = MOI.get(model, MOI.ConstraintFunction(), c)
    @test foo ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 0.0)
    _test_model_solution(
        model,
        config;
        objective_value = 0.5,
        variable_primal = [(x, 0.5)],
        # constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, -0.5)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_func_scalaraffine_lessthan),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.5]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-0.5],
        ),
    )
    return
end

"""
    test_modification_func_vectoraffine_nonneg(
        model::MOI.ModelLike,
        config::Config,
    )

Test setting the function in a VectorAffineFunction-in-Nonnegatives
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_func_vectoraffine_nonneg(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x, y
    minobjective: 1.0x + 2.0y
    c: [1.0x, 2.0y] in Nonnegatives(2)
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    c = MOI.get(model, MOI.ConstraintIndex, "c")
    _test_model_solution(
        model,
        config;
        objective_value = 0.0,
        variable_primal = [(x, 0.0), (y, 0.0)],
        constraint_primal = [(c, [0.0, 0.0])],
    )
    MOI.set(
        model,
        MOI.ConstraintFunction(),
        c,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, y)),
            ],
            [-1.0, -1.5],
        ),
    )
    foo = MOI.get(model, MOI.ConstraintFunction(), c)
    @test foo ≈ MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, y)),
        ],
        [-1.0, -1.5],
    )
    _test_model_solution(
        model,
        config;
        objective_value = 2.5,
        variable_primal = [(x, 1.0), (y, 0.75)],
        constraint_primal = [(c, [0.0, 0.0])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_func_vectoraffine_nonneg),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.0, 0.0]),
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 0.75]),
        ),
    )
    return
end

"""
    test_modification_const_vectoraffine_nonpos(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying the constant term in a VectorAffineFunction-in-Nonpositives
constraint. If `config.solve=true` confirm that it solves correctly, and if
`config.duals=true`, check that the duals are computed correctly.
"""
function test_modification_const_vectoraffine_nonpos(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x, y
    maxobjective: 1.0x + 2.0y
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, y)),
            ],
            [0.0, 0.0],
        ),
        MOI.Nonpositives(2),
    )
    _test_model_solution(
        model,
        config;
        objective_value = 0.0,
        variable_primal = [(x, 0.0), (y, 0.0)],
        constraint_primal = [(c, [0.0, 0.0])],
    )
    MOI.modify(model, c, MOI.VectorConstantChange([-1.0, -1.5]))
    _test_model_solution(
        model,
        config;
        objective_value = 2.5,
        variable_primal = [(x, 1.0), (y, 0.75)],
        constraint_primal = [(c, [0.0, 0.0])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_const_vectoraffine_nonpos),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.0, 0.0]),
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 0.75]),
        ),
    )
    return
end

"""
    test_modification_multirow_vectoraffine_nonpos(
        model::MOI.ModelLike,
        config::Config,
    )

Test modifying the variable coefficients in a
VectorAffineFunction-in-Nonpositives constraint. If `config.solve=true` confirm
that it solves correctly.
"""
function test_modification_multirow_vectoraffine_nonpos(
    model::MOI.ModelLike,
    config::Config,
)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    c = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x)),
            ],
            [-1.0, -1.0],
        ),
        MOI.Nonpositives(2),
    )
    _test_model_solution(
        model,
        config;
        objective_value = 0.5,
        variable_primal = [(x, 0.5)],
        constraint_primal = [(c, [-0.5, 0.0])],
    )
    MOI.modify(model, c, MOI.MultirowChange(x, [(1, 4.0), (2, 3.0)]))
    _test_model_solution(
        model,
        config;
        objective_value = 0.25,
        variable_primal = [(x, 0.25)],
        constraint_primal = [(c, [0.0, -0.25])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_multirow_vectoraffine_nonpos),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [0.5])),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.25]),
        ),
    )
    return
end

"""
    test_modification_const_scalar_objective(
        model::MOI.ModelLike,
        config::Config,
    )

Test the constant of a scalaraffine objective. If `config.solve=true` confirm
that it solves correctly.
"""
function test_modification_const_scalar_objective(
    model::MOI.ModelLike,
    config::Config,
)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 2.0)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.add_constraint(model, x, MOI.LessThan(1.0))
    _test_model_solution(
        model,
        config;
        objective_value = 3.0,
        variable_primal = [(x, 1.0)],
    )
    MOI.modify(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarConstantChange(3.0),
    )
    _test_model_solution(
        model,
        config;
        objective_value = 4.0,
        variable_primal = [(x, 1.0)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_const_scalar_objective),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [1.0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [1.0])),
    )
    return
end

"""
    test_modification_coef_scalar_objective(
        model::MOI.ModelLike,
        config::Config
    )

Test modifying a variable coefficient in a scalaraffine objective. If
`config.solve=true` confirm that it solves correctly.
"""
function test_modification_coef_scalar_objective(
    model::MOI.ModelLike,
    config::Config,
)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.add_constraint(model, x, MOI.LessThan(1.0))
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
    )
    MOI.modify(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarCoefficientChange(x, 3.0),
    )
    _test_model_solution(
        model,
        config;
        objective_value = 3.0,
        variable_primal = [(x, 1.0)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_coef_scalar_objective),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [1.0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [1.0])),
    )
    return
end

"""
    test_modification_delete_variable_with_single_variable_obj(
        model::MOI.ModelLike,
        config::Config,
    )

Test deleting a variable that is the objective function.
"""
function test_modification_delete_variable_with_single_variable_obj(
    model::MOI.ModelLike,
    config::Config,
)
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    c = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.delete(model, y)
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, 1.0)],
    )
    return
end

function setup_test(
    ::typeof(test_modification_delete_variable_with_single_variable_obj),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
            (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0],
        ),
    )
    return
end

"""
    test_modification_delete_variables_in_a_batch(
        model::MOI.ModelLike,
        config::Config,
    )

Test deleting many variables in a batch (i.e. using the delete method which
takes a vector of variable references). If `config.solve=true` confirm that it
solves correctly.
"""
function test_modification_delete_variables_in_a_batch(
    model::MOI.ModelLike,
    config::Config,
)
    MOIU.loadfromstring!(
        model,
        """
    variables: x, y, z
    minobjective: 1.0 * x + 2.0 * y + 3.0 * z
    x >= 1.0
    y >= 1.0
    z >= 1.0
""",
    )
    x, y, z = MOI.get(model, MOI.ListOfVariableIndices())
    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test MOI.is_valid(model, z)
    if _supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 6.0, config)
    end
    MOI.delete(model, [x, z])
    @test !MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test !MOI.is_valid(model, z)
    if _supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 2.0, config)
    end
    return
end

function setup_test(
    ::typeof(test_modification_delete_variables_in_a_batch),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 1.0, 1.0]),
            MOI.FEASIBLE_POINT,
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
        ),
    )
    return
end

"""
    test_modification_affine_deletion_edge_cases(
        model::MOI.ModelLike,
        config::Config,
    )

Test various edge cases relating to deleting affine constraints. This requires
    + ScalarAffineFunction-in-LessThan; and
    + VectorAffineFunction-in-Nonpositives.
"""
function test_modification_affine_deletion_edge_cases(
    model::MOI.ModelLike,
    config::Config,
)
    x = MOI.add_variable(model)
    # helpers. The function 1.0x + 0.0
    saf = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    vaf = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
        [0.0],
    )
    vaf2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
        [-2.0],
    )
    # max x
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        saf,
    )
    # test adding a VectorAffineFunction -in- LessThan
    c1 = MOI.add_constraint(model, vaf, MOI.Nonpositives(1))
    _test_model_solution(
        model,
        config;
        objective_value = 0.0,
        constraint_primal = [(c1, [0.0])],
    )
    # test adding a ScalarAffineFunction -in- LessThan
    c2 = MOI.add_constraint(model, saf, MOI.LessThan(1.0))
    _test_model_solution(
        model,
        config;
        objective_value = 0.0,
        constraint_primal = [(c1, [0.0]), (c2, 0.0)],
    )
    # now delete the VectorAffineFunction
    MOI.delete(model, c1)
    @test_throws MOI.InvalidIndex{typeof(c1)} MOI.delete(model, c1)
    try
        MOI.delete(model, c1)
    catch err
        @test err.index == c1
    end
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        constraint_primal = [(c2, 1.0)],
    )
    # add a different VectorAffineFunction constraint
    c3 = MOI.add_constraint(model, vaf2, MOI.Nonpositives(1))
    _test_model_solution(
        model,
        config;
        objective_value = 1.0,
        constraint_primal = [(c2, 1.0), (c3, [-1.0])],
    )
    # delete the ScalarAffineFunction
    MOI.delete(model, c2)
    _test_model_solution(
        model,
        config;
        objective_value = 2.0,
        constraint_primal = [(c3, [0.0])],
    )
    return
end

function setup_test(
    ::typeof(test_modification_affine_deletion_edge_cases),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [0.0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [0.0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [1.0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [1.0])),
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [2.0])),
    )
    return
end


text:
module Halite

using Memento

include("init.jl")
include("constants.jl")
include("commands.jl")
include("common.jl")
include("entity.jl")
include("game_map.jl")
include("networking.jl")
include("positionals.jl")

"""
    read_input()

Reads input from stdin, shutting down logging and exiting if an EOFError occurs
"""
function read_input()
    try
        return readline()
    catch err
        log(err)
        rethrow(err)
    end
end

export Direction, logging

end #module Halite

max_stars_repo_name: Affos/Halite-III

max_issues_repo_name: Affos/Halite-III

max_forks_repo_name: Affos/Halite-III


text: # TextAnalysis.jl does not need wrapping at this point (some things do not work, can be easily wrapped in pipes)


text:
@testset "importance sampling" begin

    @gen function model()
        x = @trace(normal(0, 1), :x)
        @trace(normal(x, 1), :y)
    end

    @gen function proposal()
        @trace(normal(0, 2), :x)
    end

    y = 2.
    observations = choicemap()
    set_value!(observations, :y, y)

    n = 4

    (traces, log_weights, lml_est) = importance_sampling(model, (), observations, n)
    @test length(traces) == n
    @test length(log_weights) == n
    @test isapprox(logsumexp(log_weights), 0., atol=1e-14)
    @test !isnan(lml_est)
    for trace in traces
        @test get_choices(trace)[:y] == y
    end

    (traces, log_weights, lml_est) = importance_sampling(model, (), observations, proposal, (), n)
    @test length(traces) == n
    @test length(log_weights) == n
    @test isapprox(logsumexp(log_weights), 0., atol=1e-14)
    @test !isnan(lml_est)
    for trace in traces
        @test get_choices(trace)[:y] == y
    end

    (trace, lml_est) = importance_resampling(model, (), observations, n)
    @test isapprox(logsumexp(log_weights), 0., atol=1e-14)
    @test !isnan(lml_est)
    @test get_choices(trace)[:y] == y

    (trace, lml_est) = importance_resampling(model, (), observations, proposal, (), n)
    @test isapprox(logsumexp(log_weights), 0., atol=1e-14)
    @test !isnan(lml_est)
    @test get_choices(trace)[:y] == y
end


text:
const AbstractStringOrSymbol = Union{AbstractString, Symbol}

# Verifies if the object has keys from vector
function verifykeys(
    d::Union{JSON3.Object, Params, Query},
    ks::Union{Tuple{Vararg{T}}, AbstractVector{T}}
) where {T<:AbstractStringOrSymbol}
    return all(k -> haskey(d, k), ks)
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Create arrays with DGSEM-specific structure to store the mean values and set them all to 0
function initialize_mean_values(mesh::TreeMesh{2}, equations::AbstractCompressibleEulerEquations{2},
                                dg::DGSEM, cache)
  uEltype = eltype(cache.elements)
  v_mean = zeros(uEltype, (ndims(equations), nnodes(dg), nnodes(dg), nelements(cache.elements)))
  c_mean = zeros(uEltype, (nnodes(dg), nnodes(dg), nelements(cache.elements)))
  rho_mean = zeros(uEltype, size(c_mean))
  vorticity_mean = zeros(uEltype, size(c_mean))

  return (; v_mean, c_mean, rho_mean, vorticity_mean)
end

# Create cache which holds the vorticity for the previous time step. This is needed due to the
# trapezoidal rule
function create_cache(::Type{AveragingCallback}, mesh::TreeMesh{2},
                      equations::AbstractCompressibleEulerEquations{2}, dg::DGSEM, cache)
  # Cache vorticity from previous time step
  uEltype = eltype(cache.elements)
  vorticity_prev = zeros(uEltype, (nnodes(dg), nnodes(dg), nelements(cache.elements)))
  return (; vorticity_prev)
end

# Calculate vorticity for the initial solution and store it in the cache
function initialize_cache!(averaging_callback_cache, u,
                           mesh::TreeMesh{2}, equations::AbstractCompressibleEulerEquations{2},
                           dg::DGSEM, cache)
  @unpack vorticity_prev = averaging_callback_cache

  # Calculate vorticity for initial solution
  calc_vorticity!(vorticity_prev, u, mesh, equations, dg, cache)

  return nothing
end

# Update mean values using the trapezoidal rule
function calc_mean_values!(mean_values, averaging_callback_cache, u, u_prev, integration_constant,
                           mesh::TreeMesh{2}, equations::AbstractCompressibleEulerEquations{2},
                           dg::DGSEM, cache)
  @unpack v_mean, c_mean, rho_mean, vorticity_mean = mean_values
  @unpack vorticity_prev = averaging_callback_cache

  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      vorticity = calc_vorticity_node(u, mesh, equations, dg, cache, i, j, element)
      vorticity_prev_node = vorticity_prev[i, j, element]
      vorticity_prev[i, j, element] = vorticity # Cache current vorticity for the next time step

      u_node_prim = cons2prim(get_node_vars(u, equations, dg, i, j, element), equations)
      u_prev_node_prim = cons2prim(get_node_vars(u_prev, equations, dg, i, j, element), equations)

      rho,      v1,      v2,      p      = u_node_prim
      rho_prev, v1_prev, v2_prev, p_prev = u_prev_node_prim

      c = sqrt(equations.gamma * p / rho)
      c_prev = sqrt(equations.gamma * p_prev / rho_prev)

      # Calculate the contribution to the mean values using the trapezoidal rule
      vorticity_mean[i, j, element] += integration_constant * (vorticity_prev_node + vorticity)
      v_mean[1, i, j, element]      += integration_constant * (v1_prev + v1)
      v_mean[2, i, j, element]      += integration_constant * (v2_prev + v2)
      c_mean[i, j, element]         += integration_constant * (c_prev + c)
      rho_mean[i, j, element]       += integration_constant * (rho_prev + rho)
    end
  end

  return nothing
end

end # @muladd