text:
using ExponentialAction
using FiniteDifferences
using ForwardDiff
using ReverseDiff
using Zygote

@testset "expv automatic differentiation" begin
    @testset "ForwardDiff" begin
        t = rand()
        A = randn(10, 10)
        B = randn(10)
        tjac_fd, Ajac_fd, Bjac_fd = FiniteDifferences.jacobian(
            central_fdm(5, 1), expv, t, A, B
        )
        tjac_ad = ForwardDiff.derivative(t -> expv(t, A, B), t)
        @test tjac_ad ≈ tjac_fd
        Ajac_ad = ForwardDiff.jacobian(A -> expv(t, A, B), A)
        @test Ajac_fd ≈ Ajac_fd
        Bjac_ad = ForwardDiff.jacobian(B -> expv(t, A, B), B)
        @test Bjac_fd ≈ Bjac_fd
    end

    @testset "ReverseDiff" begin
        t = rand()
        A = randn(10, 10)
        B = randn(10)
        tjac_fd, Ajac_fd, Bjac_fd = FiniteDifferences.jacobian(
            central_fdm(5, 1), expv, t, A, B
        )
        tjac_ad = vec(ReverseDiff.jacobian(t -> expv(first(t), A, B), [t]))
        @test tjac_ad ≈ tjac_fd
        Ajac_ad = ReverseDiff.jacobian(A -> expv(t, A, B), A)
        @test Ajac_fd ≈ Ajac_fd
        Bjac_ad = ReverseDiff.jacobian(B -> expv(t, A, B), B)
        @test Bjac_fd ≈ Bjac_fd
    end

    @testset "Zygote" begin
        @testset for T in (Float64, ComplexF64)
            t = rand(T)
            A = randn(T, 10, 10)
            B = randn(T, 10)
            Ȳ = randn(T, 10)
            t̄_fd, Ā_fd, B̄_fd = j′vp(central_fdm(5, 1), expv, Ȳ, t, A, B)
            Y, back = Zygote.pullback(expv, t, A, B)
            @test Y ≈ expv(t, A, B)
            t̄_ad, Ā_ad, B̄_ad = back(Ȳ)
            @test t̄_ad ≈ t̄_fd
            @test Ā_ad ≈ Ā_fd
            @test B̄_ad ≈ B̄_fd
        end
    end
end


text:
using Test
include("../src/Pinecone.jl")

@testset "PineconeIndex Type Tests" begin
   TESTINDEX = "test"
   pineconeindex = Pinecone.PineconeIndex(TESTINDEX)
   @test typeof(pineconeindex) == Pinecone.PineconeIndex
   @test pineconeindex.indexname == TESTINDEX
end;

@testset "PineconeContext Type Tests" begin
   KEY = "asdf-asdfadfs-adsfafds"
   ENV = "us.east.gcp"
   PROJECT = "testproject"
   pineconecontext = Pinecone.PineconeContext(KEY, ENV, PROJECT)
   @test typeof(pineconecontext) == Pinecone.PineconeContext
   @test pineconecontext.apikey == KEY
   @test pineconecontext.cloudenvironment == ENV
   @test pineconecontext.projectname == PROJECT
end

@testset "PineconeVector Type Tests" begin
   VID = "vectorid"
   VALUES = [0.1, 0.2, 0.3, 0.4]
   META = Dict{String, Any}("foo"=>"bar", "key"=>"value", "numeric"=>1234)
   pineconevector = Pinecone.PineconeVector(VID, VALUES, META)
   @test typeof(pineconevector) == Pinecone.PineconeVector
   @test pineconevector.id == VID
   @test length(pineconevector.values) == 4
   @test typeof(pineconevector.values) == Vector{Float64}
   @test length(pineconevector.metadata) == 3
   @test typeof(pineconevector.metadata) == Dict{String,Any}
end

ext: jl


text:
predict(x, Q) = transpose(x) * Q * x

loss(x, y, Q) = (predict(x, Q) .-y).^2

function make_df(X, Q, penalty, binarization, domains)
	df = DataFrame()
	for (i,x) in enumerate(X)
		if i == 1
			df = DataFrame(transpose(x), :auto)
		else
			push!(df, transpose(x))
		end
	end

	dim = length(df[1,:])

    if binarization == :none
        df[!,:penalty] = map(r -> penalty(Vector(r)), eachrow(df))
        df[!,:predict] = map(r -> predict(Vector(r), Q), eachrow(df[:, 1:dim]))
    else
        df[!,:penalty] = map(
            r -> penalty(binarize(Vector(r), domains; binarization)),
            eachrow(df)
        )
        df[!,:predict] = map(
            r -> predict(binarize(Vector(r), domains; binarization), Q),
            eachrow(df[:, 1:dim])
        )
    end

	min_false = minimum(
        filter(:penalty => >(minimum(df[:,:penalty])), df)[:,:predict];
        init = typemax(Int)
    )
    df[!,:shifted] = df[:,:predict] .- min_false
    df[!,:accurate] = df[:, :penalty] .* df[:,:shifted] .≥ 0.

	return df
end

function oversample(X, f)
    X_true = Vector{eltype(X)}()
    X_false = Vector{eltype(X)}()

    μ = minimum(f, X)

    foreach(x -> push!(f(x) == μ ? X_true : X_false, x), X)

    b = length(X_true) > length(X_false)
    Y = reverse(b ? X_true : X_false)
    it = Iterators.cycle(b ? X_false : X_true)

    Z = Vector{eltype(X)}()
    l = length(Y)
    for (i, x) in enumerate(it)
        push!(Z, x, Y[i])
        i == l && break
    end

    return Z
end

function preliminaries(X, domains, binarization)
    if binarization==:none
        n = length(first(X))
        return X, zeros(n,n)
    else
        Y = map(x -> collect(binarize(x, domains; binarization)), X)
        n = length(first(Y))
        return Y, zeros(n,n)
    end
end

function preliminaries(X, _)
    n = length(first(X))
    return X, zeros(n,n)
end

function train!(Q, X, penalty, η, precision, X_test, oversampling, binarization, domains)
    θ = params(Q)
    for x in (oversampling ? oversample(X, penalty) : X)
        grads = gradient(() -> loss(x, penalty(x), Q), θ)
        Q .-= η * grads[Q]
    end

    Q[:,:] = round.(precision*Q)

    df = make_df(X_test, Q, penalty, binarization, domains)
    return pretty_table(describe(df[!, [:penalty, :predict, :shifted, :accurate]]))
end

function train(
    X,
    domains::Vector{D},
    penalty;
    η = .001,
    precision = 5,
    X_test = X,
    binarization = :none,
    oversampling = false,
) where {D <: DiscreteDomain}
    Y, Q = preliminaries(X, domains, binarization)
    return train!(Q, Y, penalty, η, precision, X_test, oversampling, binarization, domains)
end

function train(
    X, penalty;
    η = .001, precision = 5, X_test = X, binarization = :one_hot, oversampling = false,
)
    return train(
        X, to_domains(X), penalty;
        η, precision, X_test, binarization, oversampling,
    )
end

function train(
    X, domain_size::Int, penalty;
    η = .001, precision = 5, X_test = X, binarization = :one_hot, oversampling = false,
)
    return train(
        X, to_domains(X, domain_size), penalty;
        η, precision, X_test, binarization, oversampling
    )
end

function train(
    X, domain_sizes::Vector{Int}, penalty;
    η = .001, precision = 5, X_test = X, binarization = :one_hot, oversampling = false,
)
    return train(
        X, to_domains(X, domain_sizes), penalty;
        η, precision, X_test, binarization, oversampling
    )
end

function train(
    X, d::D, penalty;
    η = .001, precision = 5, X_test = X, binarization = :one_hot, oversampling = false,
) where {D <: DiscreteDomain}
    return train(
        X, to_domains(X, d), penalty;
        η, precision, X_test, binarization, oversampling
    )
end


text:
a₀ = [1.0 2.0; 2.0 3.0]
a₁ = [1.0 0.0; 0.0 1.0]
a₂ = [4.0 5.0; 5.0 6.0]
a = [a₀,a₁,a₂]

b₀ = [1.0 5.0; 5.0 5.0]
b₁ = [2.0 3.0; 3.0 4.0]
b₂ = [5.0 1.0; 1.0 6.0]
b = [b₀,b₁,b₂]

c₀ = [1.0 3.0; 3.0 1.0]
c₁ = [2.0 0.0; 0.0 3.0]
c₂ = [4.0 5.0; 5.0 1.0]
c = [c₀,c₁,c₂]

n = 2

using Convex
using SCS

x = Variable(n)

A = a[1] + x[1]*a[2] + x[2]*a[3]
B = b[1] + x[1]*b[2] + x[2]*b[3]
C = c[1] + x[1]*c[2] + x[2]*c[3]

for λ = 5.0:-0.1:0.7
    constr1 = (λ*B-A ⪰ 0)
    constr2 = (B ⪰ 0)
    constr3 = (C ⪰ 0)
    problem = satisfy([constr1,constr2,constr3])
    solve!(problem,SCS.Optimizer(verbose=false))
    @show λ
    @show evaluate(x)
end

using LinearAlgebra

evaluate(A)


text:
    ####### SIMULATION OF MARGINAL ###########################################

    ### +++ INITIALIZATION +++

    ### import packages
    using Random, LinearAlgebra, RDatasets, MixedModels, Distributions, CategoricalArrays

    ### define functions
    vcov2 = function(x, xbar, σ, τ, n, γ)
        m = length(xbar)
        a = γ / τ^2 * m # 1-1 entry
        d = (sum(x.^2) - n * γ * sum(xbar.^2) ) / σ^2 # 2-2 entry
        b = γ / τ^2 / n * sum(x) # 2-1 entry
        [d -b; -b a] / (a * d - b^2) # return inverse (X^t V^-1 X)^-1
    end

    beta = function(vcov1, x, xbar, y, ybar, σ, τ, n, γ)
        c1 = γ / τ^2 * sum(ybar) # = 1^t V^-1 y
        c2 = (x' * y - n * γ * sum(ybar .* xbar))/σ^2 # = x^t V^-1 y
        vcov1 * [c1 c2]'  # (X^tV^-1X)^-1 X^tV^-1y
    end

    ### fix parameters
    p = 2
    nsim = 10000
    Δ = [0:0.1:2;]
    Random.seed!(123)

    ### +++ PER SPECIFICATION SETTING +++

    τ0 = 8^0.5
    σ0 = 4^0.5
    m = 5
    n = 5
    REdistr = "normal"

    ### variable parameters
    N = m * n

    β0 = [-4.9, 0.03]
    x = 10 * randn(N)
    xbar = vec(sum(reshape(x, n, :), dims=1) / n)

    ### conditional evaluation
    if REdistr == "normal"
        v = τ0 * randn(m)
    elseif REdistr == "tdist"
        v = rand(TDist(2), m)
    elseif REdistr == "chi2"
        v = rand(Chisq(τ0^2 / 2), m) .- τ0^2 / 2 - m
    end
    u = β0[1] .+ β0[2] * x + repeat(v, inner = n)
    μ0 = β0[1] .+ β0[2] .* xbar + v

    ### obtain conditional variance
    σ1 = zeros(10000, 1)
    τ1 = zeros(10000, 1)
    for i in 1:10000
        data = DataFrame(x = x, v = repeat(categorical(v), inner = n), y = u + σ0 * randn(N))
        model = fit(MixedModel, @formula(y ~ x + (1 | v)), data)
        τ1[i] = VarCorr(model).σρ.v.σ[1]
        σ1[i] = VarCorr(model).s
    end
    σ1 = mean(σ1) # should be (close to) σ0, but theory for conditional convergence is missing
    τ1 = mean(τ1) # should be (close to) τ0, but theory for conditional convergence is missing

    ### +++ KNOWN DELTA +++

    ###     MARGINAL COVARIANCE MATRIX

    ### helpful parameters
    γ0 = τ1^2 / (τ1^2 + σ1^2 / n)
    vcov0 = vcov2(x, xbar, σ1, τ1, n, γ0) # = (X^tV^-1X)^-1
    d0 = (1 - γ0) .* hcat(fill(1, m), xbar) # d' = l' - b'X
    dd0 = d0' * d0

    ### assemble covariance matrix parts
    g1 = γ0 * σ1^2 / n

    ### assemble inverse of covariance matrix with Woodbury
    Σ0diag = g1
    Σ0i = (I - d0 * inv(Σ0diag * inv(vcov0) + dd0) * d0' ) / Σ0diag

    ###     CONDITIONAL COVARIANCE MATRIX
    l1 = γ0^2 / n * σ1^2

    ### assemble inverse of covariance matrix with Woodbury
    cΣ0diag = l1
    cΣ0i = (I - d0 * inv(cΣ0diag * σ1^2 *
                inv(vcov0 * (2 * γ0 * I + n^2 * dd0 * vcov0)) + dd0) * d0' ) / cΣ0diag

    cΣ0sq = eigvecs(cΣ0i)' * diagm(eigvals(cΣ0i).^0.5) * eigvecs(cΣ0i)

    ###     NONCENTRALITY PARAMETER
    AZv = (γ0 - 1) .* v + n / τ1^2 .* d0 * vcov0 * d0' * v
    λ0 = AZv' * cΣ0i * AZv

    ###     QUANTILES
    qc0 = cquantile(NoncentralChisq(m, λ0), 0.05) # quantile of alpha = 0.05
    qm = cquantile(Chisq(m), 0.05) # quantile of alpha = 0.05

    ### simulation
    Coverage = zeros(5, length(Δ))
    λ = zeros(nsim, 2)
    Volume = ones(nsim, 5)

    for i in 1:nsim

        ### +++ EMPIRICAL +++

        ### variable data
        y = u + σ0 * randn(N) # conditional
        ybar = vec(sum(reshape(y, n, :), dims=1) / n)
        data = DataFrame(x = x, v = repeat(categorical(v), inner = n), y = y)

        ### fit model
        model = fit(MixedModel, @formula(y ~ x + (1 | v)), data)

        ### make sure the variance parameters are positive
        τ = maximum([VarCorr(model).σρ.v.σ[1], 0.001])
        σ = maximum([VarCorr(model).s, 0.001])

        ###     MARGINAL

        ### helpful parameters
        γ = τ^2 / (τ^2 + σ^2 / n)
        vcov1 = vcov2(x, xbar, σ, τ, n, γ) # = (X^tV^-1X)^-1
        d = (1-γ) .* hcat(fill(1, m), xbar) # d' = l' - b'X
        dd = d' * d
        a = n * τ^2 + σ^2

        ### assemble diagonal matrix parts
        g1 = γ * σ^2 / n
        g3 = 2 * σ^4 / τ^2 * γ / N/ (n - 1)

        ### assemble inverse of covariance matrix with Woodbury
        Σdiag = g1 + 2 * g3
        #Σ = Σdiag * I + d * vcov1 * d'
        Σi = (I - d * inv(Σdiag * inv(vcov1) + dd) * d' ) / Σdiag

        ###     CONDITIONAL

        ### assemble diagonal matrix parts
        l1 = γ^2 / n * σ^2

        l3_1 = 2 * σ^2 * τ^2 / N/ a^3
        l3_2 = 8 / N/ (n - 1) / a^4 / n^2 * τ^2 * σ^2 * (
            + (a^2 - σ^2) / n * (
                + σ^6 / a^2 * n
                - (σ^4  / a^2) * (τ + a / 2 * n)
                - (τ^2 - a / 2 + σ^2 * (1 + n^2)) * (n - 1))
            + σ^10 / a^2 * (n^3 - (n - 1)^2 + 1 / n)
            + a * σ^4 * (n - 1 + 1 / n / 2)  +
            - σ^6 * (2 * n - 1 + 1 / n) +
            + σ^8 / a * (1 + (n-1) * (n^2 + 1) / n) / 2)
        l3_3 = - 8 / N/ n * τ^2 * σ^4 / a^4 * (n - 1) * (a^2 - σ^4)
        l3_4 = - 8 / N/ n * τ^2 * σ^4 / a^3 / (n - 1) * (
            + a * ((n - 1) * n - 1)
            + σ^2 * (n + 1))
        l3 = 2 * (l3_1 + l3_2 + l3_3 + l3_4)

        l4 = 2 * σ^4 / N / a^2 * ( 1 / n + ((n - 1) * τ^2 +
                (σ^2 + τ^2)^2) / ((n - 1) * a^2))

        l5 = 2 * n / (n - 1) / N * σ^4 * (σ^2 - 2 * n * τ^2) / a^2

        ### assemble inverse of covariance matrix with Woodbury
        cΣdiag = maximum([l1 + l3 + l4 - l5, 1/100000])
        cΣi = (I - d * inv(cΣdiag * σ^2 *
                    inv(vcov1 * (2 * γ * I + n^2 * dd * vcov1)) + dd) * d' ) / cΣdiag

        ###     COMPUTE MEAN
        β1 = beta(vcov0, x, xbar, y, ybar, σ1, τ1, n, γ0)
        μ1 = (1 - γ0) * (β1[1] .+ β1[2] .* xbar) + γ0 * ybar

        β2 = beta(vcov1, x, xbar, y, ybar, σ, τ, n, γ) # = fixef(model)
        μ2 = (1 - γ) * (β2[1] .+ β2[2] .* xbar) + γ * ybar

        ###     NONCENTRALITY PARAMETER
        ay = d0 * β1 - (1 - γ0) * ybar
        λ[i, 1] = maximum([(ay' * cΣ0i * ay)[1] - tr(cΣ0i) * γ0^2 / n^3 / τ1^4 * σ1^2, 0])

        ay = d * β2 - (1 - γ) * ybar
        λ[i, 2] = maximum([(ay' * cΣi * ay)[1] - tr(cΣi) * σ^2 / a^2 / n, 0])

        ###     QUANTILES
        qc1 = cquantile(NoncentralChisq(m, λ[i, 1]), 0.05) # quantile of alpha = 0.05
        qc2 = cquantile(NoncentralChisq(m, λ[i, 2]), 0.05) # quantile of alpha = 0.05

        ###     COMPUTE STATISTICS
        t1 = (μ1 - μ0) .+ Δ'
        t2 = (μ2 - μ0) .+ Δ'
        Coverage[1,:] += (sum(t1' .* (t1' * Σ0i), dims = 2) .< qm) / nsim
        Coverage[2,:] += (sum(t2' .* (t2' * Σi), dims = 2) .< qm) / nsim
        tmp1 = sum(t1' .* (t1' * cΣ0i), dims = 2)
        Coverage[3,:] += (tmp1 .< qc0) / nsim
        Coverage[4,:] += (tmp1 .< qc1) / nsim
        Coverage[5,:] += (sum(t2' .* (t2' * cΣi), dims = 2) .< qc2) / nsim

        ###     COMPUTE RELATIVE VOLUME
        # (via approximations as the determinant of large matrices is numerically unstable)
        Volume[i, 3] = (Σ0diag / cΣ0diag * qm / qc0)^(-m / 2)
        Volume[i, 4] = (Σ0diag / cΣ0diag * qm / qc1)^(-m / 2)
        Volume[i, 5] = (Σdiag / cΣdiag * qm / qc2)^(-m / 2)

    end


text:
import DiskArrays
const AllowedXY = Union{Integer,Colon,AbstractRange}
const AllowedBand = Union{Integer,Colon,AbstractArray}

"""
    RasterDataset(dataset::AbstractDataset)

This data structure is returned by the [`ArchGDAL.readraster`](@ref) function and
is a wrapper for a GDAL dataset. This wrapper is to signal the
user that the dataset should be treated as a 3D AbstractArray
where the first two dimensions correspond to longitude and latitude
and the third dimension corresponds to different raster bands.

As it is a wrapper around a GDAL Dataset, it supports the usual
raster methods for a GDAL Dataset such as `getgeotransform`,
`nraster`, `getband`, `getproj`, `width`, and `height`. As it
is also a subtype of `AbstractDiskArray{T,3}`, it supports the
following additional methods: `readblock!`, `writeblock!`,
`eachchunk`, `haschunks`, etc.
This satisfies the DiskArray interface, allowing us to
be able to index into it like we would an array.

Constructing a RasterDataset will error if the raster bands do not
have all the same size and a common element data type.
"""
struct RasterDataset{T,DS} <: AbstractDiskArray{T,3}
    ds::DS
    size::Tuple{Int,Int,Int}
end

function RasterDataset(ds::AbstractDataset)
    iszero(nraster(ds)) && throw(ArgumentError("The Dataset does not contain any raster bands"))
    s = _common_size(ds)
    RasterDataset{_dataset_type(ds), typeof(ds)}(ds, s)
end

# Forward a few functions
# Here we try to include all functions that are relevant
# for raster-like datasets.
for f in (:getgeotransform, :nraster, :getband, :getproj,
    :width, :height, :destroy, :getdriver, :filelist, :listcapability,
    :ngcp, :copy, :write, :testcapability, :setproj!, :buildoverviews!)
    eval(:($(f)(x::RasterDataset, args...; kwargs...) = $(f)(x.ds, args...; kwargs...)))
end

# Here we need to special-case, because source and dest might be rasters
copywholeraster(x::RasterDataset,y::AbstractDataset;kwargs...) = copywholeraster(x.ds,y;kwargs...)
copywholeraster(x::RasterDataset,y::RasterDataset;kwargs...) = copywholeraster(x.ds,y.ds;kwargs...)
copywholeraster(x::AbstractDataset,y::RasterDataset;kwargs...) = copywholeraster(x.ds,y.ds;kwargs...)

"""
    _dataset_type(ds::AbstractDataset)

Tries to determine a common dataset type for all the bands
in a raster dataset.
"""
function _dataset_type(ds::AbstractDataset)
    alldatatypes = map(1:nraster(ds)) do i
        b = getband(ds, i)
        pixeltype(b)
    end
    reduce(promote_type, alldatatypes)
end

"""
    _common_size(ds::AbstractDataset)

Determines the size of the raster bands in a dataset and errors
if the sizes are not unique.
"""
function _common_size(ds::AbstractDataset)
    nr = nraster(ds)
    allsizes = map(1:nr) do i
        b = getband(ds, i)
        size(b)
    end
    s = unique(allsizes)
    length(s) == 1 || throw(DimensionMismatch("Can not coerce bands to single dataset, different sizes found"))
    Int.((s[1]..., nr))
end

getband(ds::RasterDataset, i::Integer) = getband(ds.ds, i)
unsafe_readraster(args...; kwargs...)  = RasterDataset(unsafe_read(args...; kwargs...))

"""
    readraster(s::String; kwargs...)

Opens a GDAL raster dataset. The difference to `ArchGDAL.read` is that
this function returns a `RasterDataset`, which is a subtype of `AbstractDiskArray{T,3}`,
so that users can operate on the array using direct indexing.
"""
readraster(s::String; kwargs...) = RasterDataset(read(s; kwargs...))

function DiskArrays.eachchunk(ds::RasterDataset)
    subchunks = DiskArrays.eachchunk(getband(ds, 1))
    DiskArrays.GridChunks(ds,(subchunks.chunksize..., 1))
end

DiskArrays.haschunks(::RasterDataset) = DiskArrays.Chunked()
DiskArrays.haschunks(::AbstractRasterBand) = DiskArrays.Chunked()

Base.size(band::AbstractRasterBand) = (width(band), height(band))

function DiskArrays.eachchunk(band::AbstractRasterBand)
    wI = windows(band)
    cs = (wI.blockiter.xbsize, wI.blockiter.ybsize)
    DiskArrays.GridChunks(band, cs)
end

function DiskArrays.readblock!(band::AbstractRasterBand, buffer, x::AbstractUnitRange, y::AbstractUnitRange)
    xoffset, yoffset = first(x)-1, first(y)-1
    xsize, ysize = length(x), length(y)
    read!(band, buffer, xoffset, yoffset, xsize, ysize)
end

function DiskArrays.writeblock!(band::AbstractRasterBand, value, x::AbstractUnitRange, y::AbstractUnitRange)
    xoffset, yoffset = first(x)-1, first(y)-1
    xsize, ysize = length(x), length(y)
    write!(band, value, xoffset, yoffset, xsize, ysize)
end

# AbstractDataset indexing

Base.size(dataset::RasterDataset) = dataset.size

function DiskArrays.readblock!(dataset::RasterDataset, buffer, x::AbstractUnitRange, y::AbstractUnitRange, z::AbstractUnitRange)
    buffer2 = Array(buffer)
    DiskArrays.readblock!(dataset::RasterDataset, buffer2, x, y, z)
    buffer .= buffer2
end

function DiskArrays.readblock!(dataset::RasterDataset, buffer::Array, x::AbstractUnitRange, y::AbstractUnitRange, z::AbstractUnitRange)
    xoffset, yoffset = first.((x, y)) .- 1
    xsize, ysize= length.((x, y))
    indices  = [Cint(i) for i in z]
    read!(dataset.ds, buffer, indices, xoffset, yoffset, xsize, ysize)
end

function DiskArrays.writeblock!(dataset::RasterDataset, value, x::AbstractUnitRange, y::AbstractUnitRange, bands::AbstractUnitRange)
    xoffset, yoffset = first.((x, y)) .- 1
    xsize, ysize= length.((x, y))
    indices  = [Cint(i) for i in bands]
    write!(dataset.ds, value, indices, xoffset, yoffset, xsize, ysize)
end

Base.Array(dataset::RasterDataset) = dataset[:,:,:]


text:
__precompile__()

"Dictionary with text strings for printing out mortgage terms, and term values in years, where meaningful."
mortgageTermDescriptors = Dict("1/4"=>("3 months",1//4),
	"1/2"=>("6 months",1//2),
	"1"=>("1 year",1),
	"a"=>("until paid off","a"))

"Dictionary with text strings for printouts of payment frequency."
frequencyDescriptors = Dict("m"=>"monthly", "b"=>"bimonthly", "w"=>"weekly", "t"=>"every 2 weeks")

"Dictionary with text strings for printouts of compounding."
compoundingDescriptors = Dict("s"=>"simple interest (no compounding, but accrued daily)",
	"m"=>"monthly (American) compounding",
	"c"=>"Canadian, or semi-annual, compounding")

"Dictionary for durations of payment intervals corresponding to payment frequencies."
timeValues = Dict("m"=>1//12, "b"=>1//24, "w"=>7, "t"=>14)


text:
mutable struct BiDict{K,V}
    key_to_value::Dict{K,V}
    value_to_key::Dict{V,K}
end

function Base.haskey(d::BiDict{K,V}, key::K) where {K,V}
    return haskey(d.key_to_value, key)
end

function Base.haskey(d::BiDict{K,V}, value::V) where {K,V}
    return haskey(d.value_to_key, value)
end

function Base.getindex(d::BiDict{K,V}, key::K) where {K,V}
    d.key_to_value[key]
end

function Base.getindex(d::BiDict{K,V}, value::V) where {K,V}
    d.value_to_key[value]
end

function Base.setindex!(d::BiDict{K,V}, key::K, value::V) where {K,V}
    if haskey(d.value_to_key, value)
        old_key = d.value_to_key[value]
        throw(Base.KeyError("$value is already associated with $old_key"))
    end
    d.key_to_value[key] = value
    d.value_to_key[value] = key
end

function Base.setindex!(d::BiDict{K,V}, value::V, key::K) where {K,V}
    if haskey(d.key_to_value, key)
        old_value = d.key_to_value[value]
        throw(Base.KeyError("$key is already associated with $old_value"))
    end
    d.key_to_value[key] = value
    d.value_to_key[value] = key
end

function Base.delete!(d::BiDict{K,V}, key::K) where {K,V}
    value = d.key_to_value[key]
    delete!(d.key_to_value, key)
    delete!(d.value_to_key, value)
end

function Base.delete!(d::BiDict{K,V}, value::V) where {K,V}
    key = d.value_to_key[value]
    delete!(d.key_to_value, key)
    delete!(d.value_to_key, value)
end

function Base.length(d::BiDict{K,V}) where {K,V}
    length(d.value_to_key)
end

BiDict{K,V}() where {K,V} = BiDict(Dict{K,V}(), Dict{V,K}())


text:
using ArrayViews, Test, LinearAlgebra

println("    -- testing linalg")

avparent = rand(7, 8)
bvparent = rand(8, 5)

_ab = avparent * bvparent

@test _ab == aview(avparent, :, :) * bvparent
@test _ab == unsafe_aview(avparent, :, :) * bvparent

@test _ab == avparent * aview(bvparent, :, :)
@test _ab == avparent * unsafe_aview(bvparent, :, :)

@test _ab == aview(avparent, :, :) * aview(bvparent, :, :)
@test _ab == aview(avparent, :, :) * unsafe_aview(bvparent, :, :)

for j = 1:size(bvparent,2)
    @test _ab[:,j] ≈ aview(avparent,:,:) * aview(bvparent,:,j)
    @test _ab[:,j] ≈ aview(avparent,:,:) * unsafe_aview(bvparent,:,j)
end

# @test avparent[:, 2:2:6] * bvparent[1:3, 1:2:5] == aview(avparent, :, 2:2:6) * aview(bvparent, 1:3, 1:2:5)
# @test avparent[:, 2:2:6] * bvparent[1:3, 1:2:5] == unsafe_aview(avparent, :, 2:2:6) * unsafe_aview(bvparent, 1:3, 1:2:5)

cvparent = rand(2, 3)
dvparent = rand(3, 2)

@test cvparent[1:2, 1:2:3] * dvparent[1:2:3, 1:2] == aview(cvparent, 1:2, 1:2:3) * aview(dvparent, 1:2:3, 1:2)
@test cvparent[1:2, 1:2:3] * dvparent[1:2:3, 1:2] == unsafe_aview(cvparent, 1:2, 1:2:3) * unsafe_aview(dvparent, 1:2:3, 1:2)

evparent = rand(3, 3)
fvparent = rand(3, 3)

@test evparent[1:2:3, 1:2:3] * fvparent[1:2:3, 1:2:3] == aview(evparent, 1:2:3, 1:2:3) * aview(fvparent, 1:2:3, 1:2:3)
@test evparent[1:2:3, 1:2:3] * fvparent[1:2:3, 1:2:3] == unsafe_aview(evparent, 1:2:3, 1:2:3) * unsafe_aview(fvparent, 1:2:3, 1:2:3)

@test evparent[1:2:3, 1:2:3] * fvparent[1:2, 1:2:3] == aview(evparent, 1:2:3, 1:2:3) * aview(fvparent, 1:2, 1:2:3)
@test evparent[1:2:3, 1:2:3] * fvparent[1:2, 1:2:3] == unsafe_aview(evparent, 1:2:3, 1:2:3) * unsafe_aview(fvparent, 1:2, 1:2:3)


text:
module BioXMLSciMLTest

using ModelingToolkit, OrdinaryDiffEq, CSV, DataFrames, Distributed

const cols = [:filename, :to_system, :to_problem, :to_solve, :states, :parameters, :error]
const N = length(cols)

"info on how far the lowering got"
function sciml(fn, f=ODESystem; dir="out")
    row = Vector{Any}(fill(missing, N))
    sciml!(row, fn, f) #; dir=dir, write=true)
    row
end

"""this writes a file per row so that I can pmap.
want to generalize this to use both sbml and cellml

takes a row length(cols) a filename to write the csv to
and a function `f: filename -> ODESystem`.

For CellML this would be `x -> ODESystem(CellModel(x), (0, 1.))`

At some point replace to use FileTrees and a pmapreduce or something.

Currently its really dumb. im not trying to make this perfect

tries to lower it to sys, then prob, then sol.
saves the # of states and parameters
stores the error in the last row.
"""
function sciml!(row, fn, f=ODESystem; dir="out", write=true)
    d, file  = splitdir(fn)
    name, ext = splitext(file)
    row[1] = fn
    try
        sys = f(fn)
        row[2] = true
        row[5] = length(states(sys))
        row[6] = length(parameters(sys))
        prob = ODEProblem(sys, Pair[], (0., 1.))
        row[3] = true
        sol = solve(prob, Tsit5())
        row[4] = true
    catch e
        row[end] = e
    end
    replace!(row, nothing => missing)
    write && CSV.write(joinpath(dir, name * ".csv"), DataFrame(cols .=> row))
    row
end

function scimls_pmap!(mat, fns, f=ODESystem; dir="out")
    # mat[:, :] = permutedims(reduce(hcat, @showprogress(pmap(test_cellml!, eachrow(mat), fns))))
    mat[:, :] = permutedims(reduce(hcat, pmap(sciml!, eachrow(mat), fns)))
    nothing
end

function scimls!(mat, fns, f=ODESystem; dir="out")
    # mat[:, :] = permutedims(reduce(hcat, @showprogress(pmap(test_cellml!, eachrow(mat), fns))))
    mat[:, :] = permutedims(mapreduce(hcat, sciml!, eachrow(mat), fns))
    nothing
end

"creates a df from the created row files,"
function results_df(dir="out/suite/")
    vcat(CSV.read.(readdir(dir;join=true), DataFrame)...)
end

"
THIS IS THE IMPORTANT FXN IN THE PACKAGE

f: filename -> ODESystem.

idrk how performant it is allocing the mat up front"
function files_to_sciml(fns, f=ODESystem; dir="out", pmap=false)
    mkpath(dir)
    n = length(cols)
    mat = Array{Any,2}(nothing, length(fns), n)
    pmap ? scimls_pmap!(mat, fns, f; dir=dir) : scimls!(mat, fns, f; dir=dir)
    nothing
end

export files_to_sciml, results_df, sciml, sciml!

end


text:
"""
    bbo_floor_int_searchrange(L::Int)

Generate a BlackBoxOptim SearchRange for optimizing integers on a
float scale when we are going to use floor to get back an Int.
"""
bbo_floor_int_searchrange(L::Int) = (0.499999999-L, L+0.499999999)

@inline floortoints(λ::Vector{N}) where {N<:Number} = floor.(Int, λ)
@inline floortoints(λ::Vector{Int}) = λ

"""
    nnzeroasint(λ::Vector{N})

Count number of nonzero elements of a vector when floored to integers.
"""
nnzeroasint(λ::Vector{N}) where {N<:Number} = sum(floortoints(λ) .!= 0)

"""
    nonzeroindices(λ::Vector{N}, sorted = false)

Return the indices of the nonzero elements of a vector when they are floored to integers.
"""
function nonzeroindices(λ::Vector{N}, sorted = true) where {N<:Number}
    idxs = findall(v -> v != 0, λ)
    return sorted ? idxs[sortperm(λ[idxs], rev=true)] : idxs
end

"""
    sign_predict(λ::Vector, X)

Predict class, indicated with +1 or -1, i.e. sign-based, by multiplying
a vector `λ` interpreted as integer coefficients and a design matrix `X`.
"""
sign_predict(λ::Vector{Int}, X::AbstractMatrix) = Int[(v<0 ? -1 : 1) for v in (X * λ)]
sign_predict(λ::Vector{N}, X::AbstractMatrix) where {N<:Number} =
    sign_predict(floortoints(λ), X)

"""
    sign_accuracy(yhat, y)

Accuracy with which yhat predicts the (signed) class labels in y.
"""
@inline function sign_accuracy(yhat::Vector{N1}, y::Vector{N2}) where {N1<:Number,N2<:Number}
    @assert length(yhat) == length(y)
    sum((yhat .* y) .> 0)/float(length(yhat))
end

"""
    sign_accuracy(λ, X, y)

Accuracy with which a coefficient vector `λ` predicts `y` given design matrix `X`.
"""
sign_accuracy(λ::Vector{N}, X::AbstractMatrix{N1}, y::Vector{N2}) where {N<:Number,N1<:Number,N2<:Number} =
    sign_accuracy(sign_predict(λ, X), y)

"""
    sign_misclassification_rate(yhat, y)

Misclassification rate for the prediction yhat of the (signed) class labels in y.
"""
@inline sign_misclassification_rate(yhat, y) = 1.0 - sign_accuracy(yhat, y)

"""
    aspct(v::Float64)

Convert to a percentage value.
"""
@inline aspct(v::Float64) = round(100.0*v, digits=3)

as_sign_labels(y::Vector{N}) where {N<:Number} = Int[(c>0 ? 1 : -1) for c in y]
function as_sign_labels(y::AbstractVector)
    uniquevalues = sort(unique(y))
    @assert 1 <= length(uniquevalues) <= 2
    Int[(c==uniquevalues[1] ? 1 : -1) for c in y]
end

function uptobutexcluding(P::Int, excl::AbstractVector)
    setdiff(collect(1:P), excl)
end


text:
function timezone(location::Tuple{Float64,Float64}, timestamp::DateTime=now())
    lat, lng = location
    timestamp = convert(BigInt, trunc(datetime2unix(timestamp))*big(1))
    url = "https://maps.googleapis.com/maps/api/timezone/json?location=$lat,$lng&timestamp=$timestamp&key=$(ENV["GOOGLE_MAPS_KEY"])"
    response = HTTP.get(url)
    return(JSON.parse(String(response.body)))
end


text:
using CUDA.CUDNN: handle, @workspace, cudnnTensorDescriptor, cudnnFilterDescriptor,
                  cudnnDataType, math_mode, CUDNN_DEFAULT_REORDER, CUDNN_CROSS_CORRELATION,
                  CUDNN_NOT_PROPAGATE_NAN, CUDNN_TENSOR_NCHW, dim4

cudnnversion() = CUDA.CUDNN.version()

function nnlibPadding(dims)
    pd = NNlib.padding(dims)
    if !all(pd[1:2:end] .== pd[2:2:end])
        @warn "cuDNN does not support asymmetric padding; defaulting to symmetric choice" maxlog=1
    end
    return pd[1:2:end]
end


text:
using GeometricIntegrators.Solutions

import GeometricProblems.Diagnostics
import GeometricProblems.Diagnostics: compute_invariant, compute_one_form,
       compute_invariant_error, compute_error_drift

export compute_energy, compute_toroidal_momentum,
       compute_energy_error, compute_toroidal_momentum_error,
       compute_momentum_error, compute_one_form, compute_error_drift

function convert_coordinates_RZphi_to_xyz(q::DataSeries{DT,2}) where {DT}
    x = SDataSeries(DT, q.nd, q.nt, q.ni)
    for k in 1:q.ni
        for j in 0:q.nt
            for i in 1:q.nd
                x[1,j,k] = q[1,j,k] * cos(q[3,j,k])
                x[2,j,k] = q[1,j,k] * sin(q[3,j,k])
                x[3,j,k] = q[2,j,k]
                x[4,j,k] = q[4,j,k]
            end
        end
    end
    return x
end

compute_energy(t::TimeSeries, q::DataSeries, params::NamedTuple) = compute_invariant(t, q, (t,q) -> hamiltonian(t, q, params))
compute_energy(sol::Solution, params::NamedTuple) = compute_energy(sol.t, sol.q, params)

compute_energy_error(t::TimeSeries, q::DataSeries, params::NamedTuple) = compute_invariant_error(t, q, (t,q) -> hamiltonian(t, q, params))
compute_energy_error(sol::Solution, params::NamedTuple) = compute_energy_error(sol.t, sol.q, params)

compute_toroidal_momentum(t::TimeSeries, q::DataSeries) = compute_invariant(t, q, toroidal_momentum)
compute_toroidal_momentum(sol::Solution) = compute_toroidal_momentum(sol.t, sol.q)

compute_toroidal_momentum_error(t::TimeSeries, q::DataSeries) = compute_invariant_error(t, q, toroidal_momentum)
compute_toroidal_momentum_error(sol::Solution) = compute_toroidal_momentum_error(sol.t, sol.q)

compute_momentum_error(t::TimeSeries, q::DataSeries, p::DataSeries) = GeometricProblems.Diagnostics.compute_momentum_error(t, q, p, ϑ)
compute_momentum_error(sol::Solution) = compute_momentum_error(sol.t, sol.q, sol.p)


text:
using Diversity.API
using EcoBase
using EcoBase: AbstractAssemblage, AbstractThings, AbstractPlaces
using RecipesBase
using DataFrames

# Now satisfy the EcoBase interface
import EcoBase: nplaces, placenames
nplaces(part::AbstractPartition) = countsubcommunities(part)
placenames(part::AbstractPartition) = getsubcommunitynames(part)

import EcoBase: coordinates
# Where there's no spatial information, make it up!
function coordinates(part::AbstractPartition{Nothing})
    dimx = round(Int, sqrt(countsubcommunities(part)), RoundUp)
    coords = Matrix{Float64}(undef, countsubcommunities(part), 2)
    for i in Base.OneTo(countsubcommunities(part))
        coords[i, 1] = (i - 1) % dimx + 1
        coords[i, 2] = (i - 1) ÷ dimx + 1
    end
    return coords
end

import EcoBase: nthings, thingnames
nthings(types::AbstractTypes) = counttypes(types)
thingnames(types::AbstractTypes) = gettypenames(types)

import EcoBase: occurrences, places, things
occurrences(mc::AbstractMetacommunity) = getabundance(mc)
places(mc::AbstractMetacommunity) = getpartition(mc)
things(mc::AbstractMetacommunity) = gettypes(mc)

# And use the EcoBase interface to provide a basic diversity interface
import Diversity.API: _getpartition
_getpartition(p::AbstractAssemblage) = places(p)

import Diversity.API: _gettypes
_gettypes(p::AbstractAssemblage) = things(p)

import Diversity.API: _getaddedoutput
_getaddedoutput(::AbstractThings) = nothing

import Diversity.API: _addedoutputcols
_addedoutputcols(::AbstractThings) = Dict{Symbol, Type}()

import Diversity.API: _getscale
_getscale(::AbstractAssemblage) = 1

import Diversity.API: _getabundance
_getabundance(a::AbstractAssemblage, raw::Bool) =
    raw ? occurrences(a) : occurrences(a) / sum(occurrences(a))

import Diversity.API: _getsubcommunitynames
_getsubcommunitynames(p::AbstractPlaces) = placenames(p)

import Diversity.API: _countsubcommunities
_countsubcommunities(p::AbstractPlaces) = length(_getsubcommunitynames(p))

import Diversity.API: _gettypenames
_gettypenames(p::AbstractThings, ::Bool) = thingnames(p)

import Diversity.API: _counttypes
_counttypes(t::AbstractThings, raw::Bool) = length(_gettypenames(t, raw))

import Diversity.API: _calcsimilarity
_calcsimilarity(t::AbstractThings, ::Real) =
    Matrix(1.0I, counttypes(t), counttypes(t))

import Diversity.API: _getweight
function _getweight(a::AbstractAssemblage)
    ab = _getabundance(a, false)
    w = sum(ab, dims=1)
    return reshape(w, length(w))
end

import Diversity.API: _calcabundance
function _calcabundance(::AbstractThings, a::AbstractArray)
    return a, one(eltype(a))
end

import Diversity.API: _calcordinariness
function _calcordinariness(t::AbstractThings, a::AbstractArray, ::Real)
    abundance, scale = _calcabundance(t, a)
    return _calcsimilarity(t, scale) * abundance
end

import Diversity.API: _getordinariness!
function _getordinariness!(a::AbstractAssemblage)
    return _calcordinariness(_gettypes(a), _getabundance(a, false),
                             _getscale(a))
end

import Diversity.API: _getmetaordinariness!
function _getmetaordinariness!(a::AbstractAssemblage)
    ord = sum(_getordinariness!(a), dims=2)
    return reshape(ord, length(ord))
end

import Diversity.API: _getdiversityname
_getdiversityname(::AbstractThings) = "species"

import Diversity.API: _hassimilarity
_hassimilarity(::AbstractThings) = false

RecipesBase.@recipe function f(var::DataFrame, asm::AbstractAssemblage)
    var[:diversity], getcoords(places(asm))
end


text:
import Base: rand,
             split

export Zonotope,
       scale,
       reduce_order

"""
    Zonotope{N<:Real} <: AbstractZonotope{N}

Type that represents a zonotope.

### Fields

- `center`     -- center of the zonotope
- `generators` -- matrix; each column is a generator of the zonotope

### Notes

Mathematically, a zonotope is defined as the set

```math
Z = \\left\\{ x ∈ \\mathbb{R}^n : x = c + ∑_{i=1}^p ξ_i g_i,~~ ξ_i \\in [-1, 1]~~ ∀ i = 1,…, p \\right\\},
```
where ``c \\in \\mathbb{R}^n`` is its *center* and ``\\{g_i\\}_{i=1}^p``,
``g_i \\in \\mathbb{R}^n``, is the set of *generators*.
This characterization defines a zonotope as the finite Minkowski sum of line
segments.
Zonotopes can be equivalently described as the image of a unit infinity-norm
ball in ``\\mathbb{R}^n`` by an affine transformation.

Zonotopes can be constructed in two different ways: either passing the generators as a matrix, where
each column represents a generator, or passing a list of vectors where each vector represents a generator.
Below we illustrate both ways.

The optional argument `remove_zero_generators` controls whether we remove zero
columns from the `generators` matrix. This option is active by default.

### Examples

A two-dimensional zonotope with given center and set of generators:

```jldoctest zonotope_label
julia> Z = Zonotope([1.0, 0.0], [0.1 0.0; 0.0 0.1])
Zonotope{Float64}([1.0, 0.0], [0.1 0.0; 0.0 0.1])

julia> dim(Z)
2
```
Here, each column of the second input corresponds to a generator.

We can collect its vertices using `vertices_list`:

```jldoctest zonotope_label
julia> vertices_list(Z)
4-element Array{Array{Float64,1},1}:
 [1.1, 0.1]
 [0.9, 0.1]
 [0.9, -0.1]
 [1.1, -0.1]
```

The support vector along a given direction can be computed using `σ`
(resp. the support function can be computed using `ρ`):

```jldoctest zonotope_label
julia> σ([1., 1.], Z)
2-element Array{Float64,1}:
 1.1
 0.1
```

Zonotopes admit an alternative constructor that receives a list of
vectors, each vector representing a generator:

```jldoctest
julia> Z = Zonotope(ones(2), [[1., 0.], [0., 1.], [1., 1.]])
Zonotope{Float64}([1.0, 1.0], [1.0 0.0 1.0; 0.0 1.0 1.0])

julia> Z.generators
2×3 Array{Float64,2}:
 1.0  0.0  1.0
 0.0  1.0  1.0
```
"""
struct Zonotope{N<:Real} <: AbstractZonotope{N}
    center::AbstractVector{N}
    generators::AbstractMatrix{N}

    function Zonotope(center::AbstractVector{N}, generators::AbstractMatrix{N};
                      remove_zero_generators::Bool=true) where {N<:Real}
        @assert length(center) == size(generators, 1) "the dimension of the " *
            "center ($(length(center))) and the generators " *
            "($(size(generators, 1))) need to match"
        if remove_zero_generators
            generators = delete_zero_columns!(generators)
        end
        new{N}(center, generators)
    end
end

isoperationtype(::Type{<:Zonotope}) = false

# constructor from center and list of generators
Zonotope(center::AbstractVector{N}, generators_list::AbstractVector{VN};
         remove_zero_generators::Bool=true
        ) where {N<:Real, VN<:AbstractVector{N}} =
    Zonotope(center, hcat(generators_list...);
             remove_zero_generators=remove_zero_generators)

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(Z::Zonotope{N})::Vector{N} where {N<:Real}

Return the center of a zonotope.

### Input

- `Z` -- zonotope

### Output

The center of the zonotope.
"""
function center(Z::Zonotope{N})::Vector{N} where {N<:Real}
    return Z.center
end

# --- AbstractZonotope interface functions ---

"""
   genmat(Z::Zonotope)

Return the generator matrix of a zonotope.

### Input

- `Z` -- zonotope

### Output

A matrix where each column represents one generator of the zonotope `Z`.
"""
function genmat(Z::Zonotope)
    return Z.generators
end

"""
    generators(Z::Zonotope)

Return an iterator over the generators of a zonotope.

### Input

- `Z` -- zonotope

### Output

An iterator over the generators of `Z`.
"""
function generators(Z::Zonotope)
    return generators_fallback(Z)
end

"""
    ngens(Z::Zonotope)::Int

Return the number of generators of a zonotope.

### Input

- `Z` -- zonotope

### Output

Integer representing the number of generators.
"""
ngens(Z::Zonotope)::Int = size(Z.generators, 2)

# --- LazySet interface functions ---

"""
    rand(::Type{Zonotope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
        )::Zonotope{N}

Create a random zonotope.

### Input

- `Zonotope`       -- type for dispatch
- `N`              -- (optional, default: `Float64`) numeric type
- `dim`            -- (optional, default: 2) dimension
- `rng`            -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`           -- (optional, default: `nothing`) seed for reseeding
- `num_generators` -- (optional, default: `-1`) number of generators of the
                      zonotope (see comment below)

### Output

A random zonotope.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

The number of generators can be controlled with the argument `num_generators`.
For a negative value we choose a random number in the range `dim:2*dim` (except
if `dim == 1`, in which case we only create a single generator).
"""
function rand(::Type{Zonotope};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing,
              num_generators::Int=-1
             )::Zonotope{N}
    rng = reseed(rng, seed)
    center = randn(rng, N, dim)
    if num_generators < 0
        num_generators = (dim == 1) ? 1 : rand(dim:2*dim)
    end
    generators = randn(rng, N, dim, num_generators)
    return Zonotope(center, generators)
end

# --- Zonotope functions ---

"""
    scale(α::Real, Z::Zonotope)

Concrete scaling of a zonotope.

### Input

- `α` -- scalar
- `Z` -- zonotope

### Output

The zonotope obtained by applying the numerical scale to the center and
generators of ``Z``.
"""
function scale(α::Real, Z::Zonotope)
    c = α .* Z.center
    gi = α .* Z.generators
    return Zonotope(c, gi)
end

"""
    reduce_order(Z::Zonotope, r::Union{Integer, Rational})::Zonotope

Reduce the order of a zonotope by overapproximating with a zonotope with less
generators.

### Input

- `Z` -- zonotope
- `r` -- desired order

### Output

A new zonotope with less generators, if possible.

### Algorithm

See `overapproximate(Z::Zonotope{N}, ::Type{<:Zonotope}, r::Union{Integer, Rational}) where {N<:Real}` for details.
"""
function reduce_order(Z::Zonotope{N}, r::Union{Integer, Rational})::Zonotope{N} where {N<:Real}
    return overapproximate(Z, Zonotope, r)
end

"""
    split(Z::Zonotope, j::Int)

Return two zonotopes obtained by splitting the given zonotope.

### Input

- `Z` -- zonotope
- `j` -- index of the generator to be split

### Output

The zonotope obtained by splitting `Z` into two zonotopes such that
their union is `Z` and their intersection is possibly non-empty.

### Algorithm

This function implements [Prop. 3, 1], that we state next. The zonotope
``Z = ⟨c, g^{(1, …, p)}⟩`` is split into:

```math
Z₁ = ⟨c - \\frac{1}{2}g^{(j)}, (g^{(1, …,j-1)}, \\frac{1}{2}g^{(j)}, g^{(j+1, …, p)})⟩ \\\\
Z₂ = ⟨c + \\frac{1}{2}g^{(j)}, (g^{(1, …,j-1)}, \\frac{1}{2}g^{(j)}, g^{(j+1, …, p)})⟩,
```
such that ``Z₁ ∪ Z₂ = Z`` and ``Z₁ ∩ Z₂ = Z^*``, where

```math
Z^* = ⟨c, (g^{(1,…,j-1)}, g^{(j+1,…, p)})⟩.
```

[1] *Althoff, M., Stursberg, O., & Buss, M. (2008). Reachability analysis of
nonlinear systems with uncertain parameters using conservative linearization.
In Proc. of the 47th IEEE Conference on Decision and Control.*
"""
function split(Z::Zonotope, j::Int)
    @assert 1 <= j <= ngens(Z) "cannot split a zonotope with $(ngens(Z)) generators along index $j"
    c, G = Z.center, Z.generators
    Gj = G[:, j]
    Gj_half = Gj / 2

    c₁ = c - Gj_half
    c₂ = c + Gj_half

    G₁ = copy(G)
    G₁[:, j] = Gj_half
    G₂ = copy(G₁)

    Z₁ = Zonotope(c₁, G₁)
    Z₂ = Zonotope(c₂, G₂)
    return Z₁, Z₂
end


text:
using BenchmarkTools
using Glimpse
using Glimpse.Parameters
using Glimpse: Vec3, ComponentData, SystemData, System, Diorama, Entity, AbstractComponent,
               Singleton
const Gl = Glimpse

SUITE = BenchmarkGroup()

@with_kw struct Spring <: ComponentData
    center::Point3{Float64} = zero(Point3{Float64})
    k::Float64 = 0.01
    damping::Float64 = 0.000
end

struct Spatial <: ComponentData
    p::Vec3{Float64}
    v::Vec3{Float64}
end

struct Oscillator{S<:SystemData} <: System
    data::S
end

Oscillator(m::Diorama) = Oscillator(SystemData(m, (Spatial, Spring), ()))

function Gl.update_indices!(sys::Oscillator)
    return sys.data.indices = [Gl.valid_entities(sys, Spatial, Spring)]
end

function dio_setup()
    dio = Diorama(:test, Entity[], AbstractComponent[], Singleton[], System[];
                  interactive = false)
    Gl.add_component!(dio, Spring)
    Gl.add_component!(dio, Spatial)
    Gl.add_system!(dio, Oscillator(dio))
    return dio
end

SUITE["Overseer"] = BenchmarkGroup()
SUITE["Overseer"]["setup"] = @benchmarkable dio_setup()

function create_fill_entities(dio)
    for i in 1:100
        if i % 2 == 0
            e = Gl.add_entity!(dio;
                               separate = [Spatial(Point3(30.0, 1.0, 1.0),
                                                   Vec3(1.0, 1.0, 1.0)), Spring()])
        else
            e = Gl.add_entity!(dio;
                               separate = [Spatial(Point3(30.0, 1.0, 1.0),
                                                   Vec3(1.0, 1.0, 1.0))])
        end
    end
end
function fill_entities(dio)
    for i in 1:100
        if i % 2 == 0
            Gl.set_entity_component!(dio, i,
                                     Spatial(Point3(30.0, 1.0, 1.0), Vec3(1.0, 1.0, 1.0)),
                                     Spring())
        else
            Gl.set_entity_component!(dio, i,
                                     Spatial(Point3(30.0, 1.0, 1.0), Vec3(1.0, 1.0, 1.0)))
        end
    end
end

SUITE["Overseer"]["create and fill entities"] = @benchmarkable create_fill_entities(dio) setup = (dio = dio_setup())

SUITE["Overseer"]["fill entities"] = @benchmarkable fill_entities(dio) setup = (dio = dio_setup(); create_fill_entities(dio))

SUITE["Overseer"]["update system indices"] = @benchmarkable Gl.update_system_indices!(dio) setup = (dio = dio_setup(); create_fill_entities(dio))

function Gl.update(sys::Oscillator)
    spat, spring = Gl.component(sys, Spatial), Gl.component(sys, Spring)
    for e in Gl.indices(sys)[1]
        e_spat = spat[e]
        spr = spring[e]
        v_prev = e_spat.v
        new_v = v_prev - (e_spat.p - spr.center) * spr.k - v_prev * spr.damping
        new_p = e_spat.p + v_prev * 1.0
        Gl.overwrite!(spat, Spatial(new_p, new_v), e)
    end
end

SUITE["Overseer"]["update oscillator"] = @benchmarkable Gl.update(dio.systems[1]) setup = (dio = dio_setup(); create_fill_entities(dio); Gl.update_system_indices!(dio))


text:
using MathOptInterface, LinearAlgebra, Statistics, Plots, QuantEcon, Interpolations, NLsolve, Optim, Random, IterTools, JuMP, Ipopt, MathOptInterface;

u(x) = sum(log(x[i]) for i in 1:numSectors)
w(x) = sum(5*log(x[i]) for i in 1:numSectors)

β = 0.96
α = 0.4
f(x) = x^α

numSectors = 2

numPoints1D = 2

grid = Vector{Vector{Float64}}(undef,(numPoints1D)^numSectors)

gridMax = 2
gridMin = 1

iter=1
for p in product(LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D))
    grid[iter] = collect(p)
    global iter += 1
end

wVal = w.(grid)

function T(wVal, grid, β, f ; compute_policy = false)

    wVals = zeros(numPoints1D,numPoints1D);
    for j in numPoints1D
        for i in numPoints1D
            wVals[i,j] = wVal[i+(i-1)*j]
        end
    end
    function wFunc(x,y)
        return interpolate((LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D)), wVals, Gridded(Linear()))(x,y)
    end
    Tw = zeros(length(grid))
    σ = similar(grid)
    for n in 1:length(grid)
        y = grid[n]
    end

    objectives = ( c -> sum(log(c[i]) for i in 1:numSectors) + β*wFunc(f(k[1]),f(k[2])))

    results = maximize.(objectives, 1e-10, grid) # solver result for each grid point
    Tw = Optim.maximum.(results)
    if compute_policy
        σ = Optim.maximizer.(results)
        return Tw, σ
    end

    return Tw

        #=
        modTrial = Model(Ipopt.Optimizer);
        @variable(modTrial,  c[1:numSectors])
        @variable(modTrial, k[1:numSectors], start=1.000001)
        for i in 1:numSectors
            @constraint(modTrial, 0.99*gridMin <= c[i] <= 0.999*y[i])
	    @constraint(modTrial, k[i] == y[i] - c[i])
        end
	register(modTrial, :wFunc, 2, wFunc; autodiff = true)
	register(modTrial, :f, 1, f; autodiff = true)
	@NLobjective(modTrial, Max, sum(log(c[i]) for i in 1:numSectors) + β*wFunc(f(k[1]),f(k[2])))
        optimize!(modTrial)
        Tw[n] = JuMP.objective_value(modTrial)
        if compute_policy
            σ[n] = value.(c)
            return Tw, σ
        end
    end

    return Tw
    =#
end

wVal = T(wVal, grid, β, f)
using LinearAlgebra, Statistics
using Plots, QuantEcon, Interpolations, NLsolve, Optim, Random

function T(w, grid, β, u, f; compute_policy = false)
    w_func = w
    # objective for each grid point
    objectives(c) = u(c) + β * mean(w_func.(f(c)))
    global bobjectives = objectives
    #println("here's the objectives")
    #println(objectives)
    results = maximize(objectives, 1e-10, grid) # solver result for each grid point
    global besults = results
    #println(results)

    Tw = Optim.maximum.(results)
    if compute_policy
        σ = Optim.maximizer.(results)
        return Tw, σ
    end

    return Tw
end
α = 0.4
β = 0.96
μ = 0
s = 0.1

c1 = log(1 - α * β) / (1 - β)
c2 = (μ + α * log(α * β)) / (1 - α)
c3 = 1 / (1 - β)
c4 = 1 / (1 - α * β)

# Utility
u(c) = log(c)

∂u∂c(c) = 1 / c

# Deterministic part of production function
f(k) = k^α

f′(k) = α * k^(α - 1)

# True optimal policy
c_star(y) = (1 - α * β) * y

# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log(y)
v(y) = c1 + c2 * (c3 - c4) + c4 * y
using Random
Random.seed!(42) # For reproducible results.

grid_max = 4         # Largest grid point
grid_size = 200      # Number of grid points
shock_size = 250     # Number of shock draws in Monte Carlo integral

grid = 4#range(1e-5,  grid_max, length = grid_size)
shocks = 0 #exp.(μ .+ s * randn(shock_size))
w = T(v(grid), grid, β, log, k -> k^α)
#=plt = plot(ylim = (-35,-24))
plot!(plt, grid_y, w, linewidth = 2, alpha = 0.6, label = "T(v)")
plot!(plt, v_star, grid_y, linewidth = 2, alpha=0.6, label = "v_star")
plot!(plt, legend = :bottomright)=#

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using Dates
using LibSerialPort
using SnakeIterator

include("utils.jl")

portname = LibSerialPort.get_port_list()[1]
baudrate = 115_200 # Creality Ender 3D uses a non standard baud rate

origin = Point(0, 0, 55, "origin")
plate = Point(125, 152, 29, "plate top left")
waste = Point(65, plate.y, 55, "waste container")

# let 3D printer initialize, it reboots when we setup the serial connection
sp = LibSerialPort.open(portname, baudrate)
sleep(7)

# move head up and calibrate the origin

write(sp, "G1 Z55\n")
write(sp, "G28\n")

# raise head and wait, this is useful for priming the lines
moveto(sp, origin; wait=true)

wells = create_wells(plate, well_spacing = 9.0, well_depth = 11.0)

groups = [
    wells[:, 1:6],
    wells[:, 7:12]
    ]

for group in groups
    println("Fractionating first tube")

    moveto(sp, plate)
    wait_for_key("start flow now, press any key")

    waste_start = now()
    moveto(sp, waste; wait = true)
    waste_stop = now()
    println("$(waste_stop - waste_start)")

    start_pos = first(group)
    start_pos = Point(start_pos.x, start_pos.y, "Start pos")
    moveto(sp, start_pos)
    sleep(2)

    prev_pos = start_pos

    # visit each well in a group in a snaking pattern
    for well in snake(group)
        println(well.name)
        moveto(sp, well; safe_height=plate.z)
        sleep(13)
        prev_pos = well
    end

    moveto(sp, Point(prev_pos.x, prev_pos.y + 9, waste.z*1.0, "outside"))
    moveto(sp, waste; wait=true)
end

max_stars_repo_name: tlnagy/fractionation

max_issues_repo_name: tlnagy/fractionation

max_forks_repo_name: tlnagy/fractionation


text:
using TimeseriesSurrogates.StatsBase
using TimeseriesSurrogates.DSP

"""
    surroplot(x, s; kwargs...)

Plot a timeseries `x` along with its surrogate realization `s`, and compare the
periodogram and histogram of the two time series.

## Keyword arguments
- `cx` and `cs`: Colors of the original and the surrogate time series, respectively.
- `nbins`: The number of bins for the histograms.
- `resolution`: A tuple giving the resolution of the figure.
"""
function surroplot(x, s;
        cx = "#1B1B1B", cs = ("#2DB9C5", 0.9), resolution = (500, 600),
        nbins = 50,
    )

    t = 1:length(x)
    fig = Makie.Figure(resolution = resolution)

    # Time series
    ax1, _ = Makie.lines(fig[1,1], t, x; color = cx)
    Makie.lines!(ax1, t, s; color = cs)
    # Autocorrelation
    acx = autocor(x)
    ax2, _ = Makie.lines(fig[2,1], 0:length(acx)-1, acx; color = cx)
    Makie.lines!(ax2, 0:length(acx)-1, autocor(s); color = cs)

    # Binned multitaper periodograms
    p, psurr = DSP.mt_pgram(x), DSP.mt_pgram(s)
    ax3 = Makie.Axis(fig[3,1]; yscale = log10)
    Makie.lines!(ax3, p.freq, p.power; color = cx)
    Makie.lines!(ax3, psurr.freq, psurr.power; color = cs)

    # Histograms
    ax4 = Makie.Axis(fig[4,1])
    Makie.hist!(ax4, x; label = "Original", bins = nbins, color = (cx, 0.5))
    Makie.hist!(ax4, s; label = "Surrogate", bins = nbins, color = (cs, 0.5))
    Makie.axislegend(ax4)

    ax1.xlabel = "time step"
    ax1.ylabel = "value"
    ax2.xlabel = "lag"
    ax2.ylabel = "autocor"
    ax3.xlabel = "binned freq."
    ax3.ylabel = "power"
    ax4.xlabel = "binned value"
    ax4.ylabel = "histogram"
    return fig
end
export surroplot


text:
import WhittleLikelihoodInference: MaternAcv, MaternAcv1D
@testset "MaternAcv" begin
    θ₁ = [2.0, 3/2, 0.5]
    ω₁ = 0.7
    θ₂=[2.0,0.5,0.7,
        3/2,6.3,4.7,
        0.5,2.0,1.2]
    ω₂ = 1.1
    @testset "Function sdf" begin
        @test sdf(MaternAcv1D(θ₁), ω₁)         == sdf(Matern1D(θ₁), ω₁)
        @test sdf(MaternAcv{2,3}(θ₂), ω₂)      == sdf(Matern2D(θ₂), ω₂)
    end
    @testset "Gradient sdf" begin
        @test_throws ErrorException grad_sdf(MaternAcv1D(θ₁), ω₁)
    end
    @testset "Hessian sdf" begin
        @test_throws ErrorException hess_sdf(MaternAcv1D(θ₁), ω₁)
    end
    @testset "Error handling" begin
        @test_throws ArgumentError MaternAcv1D(-1.0,1.0,1.0)
        @test_throws ArgumentError MaternAcv1D(1.0,-1.0,1.0)
        @test_throws ArgumentError MaternAcv1D(1.0,1.0,-1.0)
        @test_throws ArgumentError MaternAcv1D(ones(2))
        @test_throws ArgumentError MaternAcv1D(ones(4))

        @test_throws ArgumentError MaternAcv{2,3}([[1.0,1.2]; fill(0.5,7)])
        @test_throws ArgumentError MaternAcv{2,3}([-1; fill(0.5,8)])
        @test_throws ArgumentError MaternAcv{2,3}(fill(0.5,8))
        @test_throws ErrorException MaternAcv{2,5}(fill(0.5,9))
    end
    @testset "parameter bounds" begin
        @test_throws ErrorException lowerbounds(Matern2D)
        @test_throws ErrorException upperbounds(Matern2D)
        @test lowerbounds(Matern1D) == zeros(3)
        @test upperbounds(Matern1D) == fill(Inf,3)
    end
end


text:
function onlineplot!(scene, layout, axis_dict, stats::Series, iter, data, variable, i)
    for (j, stat) in enumerate(stats.stats)
        axis_dict[(variable, stat)] = layout[i, j] = Axis(scene, title = "$(name(stat))")
        onlineplot!(axis_dict[(variable, stat)], stat, iter, data[variable], data[:iter], i, j)
        tight_ticklabel_spacing!(axis_dict[(variable, stat)])
    end
end

function onlineplot!(scene, layout, axis_dict, stats::AbstractVector, iter, data, variable, i)
    for (j, stat) in enumerate(stats)
        axis_dict[(variable, stat)] = layout[i, j] = Axis(scene, title = "$(name(stat))")
        onlineplot!(axis_dict[(variable, stat)], stat, iter, data[variable], data[:iter], i, j)
        tight_ticklabel_spacing!(axis_dict[(variable, stat)])
    end
end

function onlineplot!(axis, stat::Symbol, args...)
    onlineplot!(axis, Val(stat), args...)
end

onlineplot!(axis, stat::Val{:mean}, args...) = onlineplot!(axis, Mean(), args...)

onlineplot!(axis, stat::Val{:var}, args...) = onlineplot!(axis, Variance(), args...)

onlineplot!(axis, stat::Val{:autocov}, args...) = onlineplot!(axis, AutoCov(20), args...)

onlineplot!(axis, stat::Val{:hist}, args...) = onlineplot!(axis, KHist(50, Float32), args...)

# Generic fallback for OnlineStat objects
function onlineplot!(axis, stat::T, iter, data, iterations, i, j) where {T<:OnlineStat}
    window = data.b
    @eval TStat = $(nameof(T))
    stat = Node(TStat(Float32))
    on(iter) do i
        stat[] = fit!(stat[], last(value(data)))
    end
    statvals = Node(MovingWindow(window, Float32))
    on(stat) do s
        statvals[] = fit!(statvals[], Float32(value(s)))
    end
    statpoints = lift(statvals; init = Point2f0.([0], [0])) do v
        Point2f0.(value(iterations), value(v))
    end
    lines!(axis, statpoints, color = std_colors[i], linewidth = 3.0)
end

function onlineplot!(axis, stat::Val{:trace}, iter, data, iterations, i, j)
    trace = lift(iter; init = [Point2f0(0, 0f0)]) do i
        Point2f0.(value(iterations), value(data))
    end
    lines!(axis, trace, color = std_colors[i]; linewidth = 3.0)
end

function onlineplot!(axis, stat::KHist, iter, data, iterations, i, j)
    nbins = stat.k
    stat = Node(KHist(nbins, Float32))
    on(iter) do i
        stat[] = fit!(stat[], last(value(data)))
    end
    hist_vals = lift(stat; init = Point2f0.(range(0, 1, length = nbins), zeros(Float32, nbins))) do h
        edges, weights = OnlineStats.xy(h)
        weights = nobs(h) > 1 ? weights / OnlineStats.area(h) : weights
        return Point2f0.(edges, weights)
    end
    barplot!(axis, hist_vals, color = std_colors[i])
end

function expand_extrema(xs)
    xmin, xmax = xs
    diffx = xmax - xmin
    xmin = xmin - 0.1 * abs(diffx)
    xmax = xmax + 0.1 * abs(diffx)
    return (xmin, xmax)
end

function onlineplot!(axis, stat::Val{:kde}, iter, data, iterations, i, j)
    interpkde = Node(InterpKDE(kde([1f0])))
    on(iter) do i
        interpkde[] = InterpKDE(kde(value(data)))
    end
    xs = lift(iter; init = range(0.0, 2.0, length = 200)) do i
        range(expand_extrema(extrema(value(data)))..., length = 200)
    end
    kde_pdf = lift(xs) do xs
        pdf.(Ref(interpkde[]), xs)
    end
    lines!(axis, xs, kde_pdf, color = std_colors[i], linewidth = 3.0)
end

function onlineplot!(axis, stat::Val{:histkde}, iter, data, iterations, i, j)
    onlineplot!(axis, KHist(50), iter, data, iterations, i, j)
    onlineplot!(axis, Val(:kde), iter, data, iterations, i, j)
end

function onlineplot!(axis, stat::AutoCov, iter, data, iterations, i, j)
    b = length(stat.cross)
    stat = Node(AutoCov(b, Float32))
    on(iter) do i
        stat[] = fit!(stat[], last(value(data)))
    end
    statvals = lift(stat; init = zeros(Float32, b + 1)) do s
        value(s)
    end
    scatter!(axis, Point2f0.([0.0, 0.0], [-0.0, 1.0]), markersize = 0.0, color = RGBA(0.0, 0.0, 0.0, 0.0)) # Invisible points to keep limits fixed
    lines!(axis, 0:b, statvals, color = std_colors[i], linewidth = 3.0)
end


text:
using HybridArrays
using StaticArrays
using Test

@test length(Test.detect_ambiguities(HybridArrays)) == 0

@testset "Inner Constructors" begin
    @test HybridArray{Tuple{2}, Int, 1, 1, Vector{Int}}((3, 4)).data == [3, 4]
    @test HybridArray{Tuple{2}, Int, 1}([3, 4]).data == [3, 4]
    @test HybridArray{Tuple{2, 2}, Int, 2}(collect(3:6)).data == collect(3:6)
    @test size(HybridArray{Tuple{4, 5}, Int, 2}(undef).data) == (4, 5)
    @test size(HybridArray{Tuple{4, 5}, Int}(undef).data) == (4, 5)

    # Bad input
    @test_throws Exception SArray{Tuple{1},Int,1}([2 3])

    # Bad parameters
    @test_throws Exception HybridArray{Tuple{1},Int,2}(undef)
    @test_throws Exception SArray{Tuple{3, 4},Int,1}(undef)

    # Parameter/input size mismatch
    @test_throws Exception HybridArray{Tuple{1},Int,2}([2; 3])
    @test_throws Exception HybridArray{Tuple{1},Int,2}((2, 3))
end

@testset "Outer Constructors" begin
    # From Array
    @test @inferred(HybridArray{Tuple{2},Float64,1}([1,2]))::HybridArray{Tuple{2},Float64,1,1} == [1.0, 2.0]
    @test @inferred(HybridArray{Tuple{2},Float64}([1,2]))::HybridArray{Tuple{2},Float64,1,1} == [1.0, 2.0]
    @test @inferred(HybridArray{Tuple{2}}([1,2]))::HybridArray{Tuple{2},Int,1,1} == [1,2]
    @test @inferred(HybridArray{Tuple{2,2}}([1 2;3 4]))::HybridArray{Tuple{2,2},Int,2,2} == [1 2; 3 4]

    # Uninitialized
    @test @inferred(HybridArray{Tuple{2,2},Int,2}(undef)) isa HybridArray{Tuple{2,2},Int,2,2}
    @test @inferred(HybridArray{Tuple{2,2},Int}(undef)) isa HybridArray{Tuple{2,2},Int,2,2}

    # From Tuple
    @test @inferred(HybridArray{Tuple{2},Float64,1,1,Vector{Float64}}((1,2)))::HybridArray{Tuple{2},Float64,1,1} == [1.0, 2.0]
    @test @inferred(HybridArray{Tuple{2},Float64}((1,2)))::HybridArray{Tuple{2},Float64,1,1} == [1.0, 2.0]
    @test @inferred(HybridArray{Tuple{2}}((1,2)))::HybridArray{Tuple{2},Int,1,1} == [1,2]
    @test @inferred(HybridArray{Tuple{2,2}}((1,2,3,4)))::HybridArray{Tuple{2,2},Int,2,2} == [1 3; 2 4]
end

@testset "HybridVector and HybridMatrix" begin
    @test @inferred(HybridVector{2}([1,2]))::HybridArray{Tuple{2},Int,1,1} == [1,2]
    @test @inferred(HybridVector{2}((1,2)))::HybridArray{Tuple{2},Int,1,1} == [1,2]
    # Reshaping
    @test @inferred(HybridVector{2}([1 2]))::HybridArray{Tuple{2},Int,1,2} == [1,2]
    # Back to Vector
    @test Vector(HybridVector{2}((1,2))) == [1,2]
    @test convert(Vector, HybridVector{2}((1,2))) == [1,2]

    @test @inferred(HybridMatrix{2,2}([1 2; 3 4]))::HybridArray{Tuple{2,2},Int,2,2} == [1 2; 3 4]
    # Reshaping
    @test @inferred(HybridMatrix{2,2}((1,2,3,4)))::HybridArray{Tuple{2,2},Int,2,2} == [1 3; 2 4]
    # Back to Matrix
    @test Matrix(HybridMatrix{2,2}([1 2;3 4])) == [1 2; 3 4]
    @test convert(Matrix, HybridMatrix{2,2}([1 2;3 4])) == [1 2; 3 4]
end

@testset "setindex" begin
    sa = HybridArray{Tuple{2}, Int, 1}([3, 4])
    sa[1] = 2
    @test sa.data == [2, 4]
end

@testset "aliasing" begin
    a1 = rand(4)
    a2 = copy(a1)
    sa1 = HybridVector{4}(a1)
    sa2 = HybridVector{4}(a2)
    @test Base.mightalias(a1, sa1)
    @test Base.mightalias(sa1, HybridVector{4}(a1))
    @test !Base.mightalias(a2, sa1)
    @test !Base.mightalias(sa1, HybridVector{4}(a2))
    @test Base.mightalias(sa1, view(sa1, 1:2))
    @test Base.mightalias(a1, view(sa1, 1:2))
    @test Base.mightalias(sa1, view(a1, 1:2))
end

@testset "back to Array" begin
    @test Array(HybridArray{Tuple{2}, Int, 1}([3, 4])) == [3, 4]
    @test Array{Int}(HybridArray{Tuple{2}, Int, 1}([3, 4])) == [3, 4]
    @test Array{Int, 1}(HybridArray{Tuple{2}, Int, 1}([3, 4])) == [3, 4]
    @test Vector(HybridArray{Tuple{4}, Int, 1}(collect(3:6))) == collect(3:6)
    @test convert(Vector, HybridArray{Tuple{4}, Int, 1}(collect(3:6))) == collect(3:6)
    @test Matrix(SMatrix{2,2}((1,2,3,4))) == [1 3; 2 4]
    @test convert(Matrix, SMatrix{2,2}((1,2,3,4))) == [1 3; 2 4]
    @test convert(Array, HybridArray{Tuple{2,2,2,2}, Int}(ones(2,2,2,2))) == ones(2,2,2,2)
    # Conversion after reshaping
    @test_broken Array(HybridMatrix{2,2,Int,1,Vector{Int}}([1,2,3,4])) == [1 3; 2 4]
end

@testset "promotion" begin
    @test @inferred(promote_type(HybridVector{1,Float64,1,Vector{Float64}}, HybridVector{1,BigFloat,1,Vector{BigFloat}})) == HybridVector{1,BigFloat,1,Vector{BigFloat}}
    @test @inferred(promote_type(HybridVector{2,Int,1,Vector{Int}}, HybridVector{2,Float64,1,Vector{Float64}})) === HybridVector{2,Float64,1,Vector{Float64}}
    @test @inferred(promote_type(HybridMatrix{2,3,Float32,2,Matrix{Float32}}, HybridMatrix{2,3,Complex{Float64},2,Matrix{Complex{Float64}}})) === HybridMatrix{2,3,Complex{Float64},2,Matrix{Complex{Float64}}}
end

@testset "dynamically sized axes" begin
    A = rand(Int, 2, 3, 4)
    B = HybridArray{Tuple{2,3,StaticArrays.Dynamic()}, Int, 3}(A)
    C = rand(Int, 2, 3, 4, 5)
    D = HybridArray{Tuple{2,3,StaticArrays.Dynamic(),StaticArrays.Dynamic()}, Int, 4}(C)
    @test size(B) == size(A)
    @test size(D) == size(C)
    @test axes(B) == (SOneTo(2), SOneTo(3), axes(A, 3))
    @test axes(B, 1) == SOneTo(2)
    @test axes(B, 2) == SOneTo(3)

    @test B[1,2,3] == A[1,2,3]
    @test B[1,:,:] == A[1,:,:]
    inds = @SVector [2, 1]
    @test B[1,inds,:] == A[1,inds,:]
    @test B[:,:,2] == A[:,:,2]
    @test B[:,:,@SVector [2, 3]] == A[:,:,[2, 3]]

    B[1,2,3] = 42
    @test B[1,2,3] == 42
    B[:,2,3] = @SVector [10, 11]
    @test B[:,2,3] == @SVector [10, 11]
    B[:,:,1] = @SMatrix [1 2 3; 4 5 6]
    @test B[:,:,1] == @SMatrix [1 2 3; 4 5 6]
    B[1,2,:] = [10, 11, 12, 13]
    @test B[1,2,:] == @SVector [10, 11, 12, 13]
    B[:,2,:] = @SMatrix [1 2 3 4; 5 6 7 8]
    @test B[:,2,:] == @SMatrix [1 2 3 4; 5 6 7 8]
    B[:,2,:] = [11 12 13 14; 15 16 17 18]
    @test B[:,2,:] == [11 12 13 14; 15 16 17 18]
    B[1,2,:] = @SVector [10, 11, 12, 13]
    @test B[1,2,:] == @SVector [10, 11, 12, 13]
    B[1,SA[2,1],SA[2,3]] = @SMatrix [10 11; 12 13]
    @test B[1,SA[2,1],SA[2,3]] == @SMatrix [10 11; 12 13]
    B[1,SA[2,1],SA[2,3]] = [14 15; 16 17]
    @test B[1,SA[2,1],SA[2,3]] == @SMatrix [14 15; 16 17]

    D[:,2,3,4] = @SVector [10, 11]
    @test D[:,2,3,4] == @SVector [10, 11]
    D[:,:,1,2] = @SMatrix [1 2 3; 4 5 6]
    @test D[:,:,1,2] == @SMatrix [1 2 3; 4 5 6]
    D[1,2,:,1] = [10, 11, 12, 13]
    @test D[1,2,:,1] == @SVector [10, 11, 12, 13]

    E = HybridArray{Tuple{}}(fill(10))
    E[] = 12
    @test E[] == 12

    @test_throws DimensionMismatch (D[:,2,3,4] = @SVector [10, 11, 11])
    @test_throws DimensionMismatch (D[:,2,3,4] = [10, 11, 11])
    @test_throws DimensionMismatch (B[1,2,:] = [10, 11])

    @test_throws ErrorException HybridArray.all_dynamic_fixed_val(typeof(B))
    @test HybridArrays.all_dynamic_fixed_val(Tuple{}) === Val(:dynamic_fixed_true)
end

include("abstractarray.jl")
include("arraymath.jl")
include("broadcast.jl")
include("linalg.jl")
include("ssubarray.jl")
if VERSION >= v"1.2"
    include("array_interface_compat.jl")
end


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

#=
    LinearSolver's interface
=#
"""
    AbstractLinearSolver

Abstract type for linear solver targeting
the resolution of the linear system ``Ax=b``.
"""
abstract type AbstractLinearSolver end

"""
    introduce(::AbstractLinearSolver)

Print the name of the linear solver.
"""
function introduce end

"""
    factorize!(::AbstractLinearSolver)

Factorize the matrix ``A`` and updates the factors
inside the `AbstractLinearSolver` instance.
"""
function factorize! end

"""
    solve!(::AbstractLinearSolver, x::AbstractVector)

Solve the linear system ``Ax = b``.

This function assumes the linear system has been
factorized previously with [`factorize!`](@ref).
"""
function solve! end

"""
    is_inertia(::AbstractLinearSolver)

Return `true` if the linear solver supports the
computation of the inertia of the linear system.
"""
function is_inertia end

"""
    inertia(::AbstractLinearSolver)

Return the inertia `(n, m, p)` of the linear system
as a tuple.

### Note
The inertia is defined as a tuple ``(n, m, p)``,
with
- ``n``: number of positive eigenvalues
- ``m``: number of negative eigenvalues
- ``p``: number of zero eigenvalues
"""
function inertia end

function improve! end

# Default function for AbstractKKTVector
function solve!(s::AbstractLinearSolver, x::AbstractKKTVector)
    solve!(s, full(x))
end

#=
    Iterator's interface
=#
abstract type AbstractIterator end

"""
    solve_refine!(x, ::AbstractIterator, b)

Solve the linear system ``Ax = b`` using iterative
refinement. The object `AbstractIterator` stores an instance
of a [`AbstractLinearSolver`](@ref) for the backsolve
operations.

### Notes
This function assumes the matrix stores in the linear solver
has been factorized previously.

"""
function solve_refine! end

# LinearSolverExceptions
struct SymbolicException <: Exception end
struct FactorizationException <: Exception end
struct SolveException <: Exception end
struct InertiaException <: Exception end
LinearSolverException=Union{SymbolicException,FactorizationException,SolveException,InertiaException}

# iterative solvers
include("backsolve.jl")

#=
    DEFAULT SOLVERS
=#

# dense solvers
include("lapack.jl")

# direct solvers
include("umfpack.jl")


text:
#####################################################################################################
## Comparison noise-modeling vs. no noise modeling. Pore model taken from:
##   https://github.com/jts/nanopolish/blob/master/etc/r9-models/r9.4_450bps.cpg.6mer.template.model
#####################################################################################################
using Random
using StatsPlots
using Distributed
using Distributions
# @everywhere using Pkg
# @everywhere Pkg.activate("./")
@everywhere using CpelNano

#####################################################################################################
# Functions
#####################################################################################################
function cos_sim(α1::Vector{Float64}, β1::Vector{Float64}, α2::Vector{Float64}, β2::Vector{Float64})::Float64

    # Create aux vars
    aux_1 = vcat(α1, β1)
    aux_2 = vcat(α2, β2)

    # Return cosine similarity
    return aux_1' * aux_2 / (sqrt(sum(aux_1.^2)) * sqrt(sum(aux_2.^2)))

end
function mse(α1::Vector{Float64}, β1::Vector{Float64}, α2::Vector{Float64}, β2::Vector{Float64})::Float64

    # Create aux vars
    aux_1 = vcat(α1, β1)
    aux_2 = vcat(α2, β2)

    # Return cosine similarity
    return sum((aux_1 .- aux_2).^2) / length(aux_1)

end
function plot_dists(μ_m::Float64, μ_u::Float64, σs::Vector{Float64})

    # Loop over σs
    ps = []
    for (i, σ) in enumerate(σs)
        # Sigmas
        σ_m = σ_u = σ

        # Plot distributions
        if i == 1
            p = plot(Normal(μ_m, σ_m), lw=3, label="methylated", title="\\sigma=$(σ)", xlabel="Current (pA)")
            plot!(p, Normal(μ_u, σ_u), lw=3, label="unmethylated")
            plot!(p, Uniform(μ_m - σ_m, μ_m + σ_m), label="")
            plot!(p, Uniform(μ_u - σ_u, μ_u + σ_u), label="")
        else
            p = plot(Normal(μ_m, σ_m), lw=3, label="", title="\\sigma=$(σ)", xlabel="Current (pA)")
            plot!(p, Normal(μ_u, σ_u), lw=3, label="")
            plot!(p, Uniform(μ_m - σ_m, μ_m + σ_m), label="")
            plot!(p, Uniform(μ_u - σ_u, μ_u + σ_u), label="")
        end

        # Push plot
        push!(ps, p)

    end

    # Return plot array
    return ps

end
function sim(μ_m::Float64, μ_u::Float64, σs::Vector{Float64}, seed::Int64, reps::Int64)

    # Seed
    Random.seed!(seed)

    # Parameter vector
    a = 1.0; b = -5.0; c = 2.0; ϕ = [a,b,c]; L = 20;

    # Configuration
    config = CpelNano.CpelNanoConfig(); config.max_em_iters = 100;
    config.max_em_init = 20; config.verbose = false;

    # Other parameters
    x = CpelNano.RegStruct(); x.L = L; x.Nl = fill(1.0, x.L);
    x.ρl = rand(0.001:0.001:0.5, x.L); x.dl = fill(10.0, x.L - 1);
    α, β = CpelNano.get_αβ_from_ϕ(ϕ, x);

    # We take 6-mers AACGAA vs AAMGAA
    M = 10; pobs = 0.8;

    # Loop over σs
    nse_mat = fill(NaN, (length(σs), reps));
    no_nse_mat = fill(NaN, (length(σs), reps));
    for (i, σ) in enumerate(σs)

        CpelNano.print_log("Working on $(σ)")

        # Sigmas
        σ_m = σ_u = σ

        # Perform reps times
        Threads.@threads for j = 1:reps

            CpelNano.print_log("Working on iteration j=$(j) (thread $(Threads.threadid()) of out $(Threads.nthreads()))")

            # Simulate data
            rs = CpelNano.cpel_samp_ont(M, α, β, pobs, μ_m, σ_m, μ_u, σ_u);
            rs.L = x.L; rs.Nl = x.Nl; rs.ρl = x.ρl; rs.dl = x.dl;
            rs_bin = deepcopy(rs); CpelNano.bin_calls!(rs_bin);

            # Run EM algorithm
            CpelNano.get_ϕhat!(rs, config)
            # CpelNano.print_log("$(rs.ϕhat)")
            CpelNano.get_ϕhat!(rs_bin, config)
            # CpelNano.print_log("$(rs_bin.ϕhat)")

            # Compute cos similarity
            if length(rs.ϕhat) > 0
                nse_α, nse_β = CpelNano.get_αβ_from_ϕ(rs.ϕhat, rs)
            end
            if length(rs_bin.ϕhat) > 0
                no_nse_α, no_nse_β = CpelNano.get_αβ_from_ϕ(rs_bin.ϕhat, rs_bin)
            end

            # Filter out problematic runs
            if mse(α, β, nse_α, nse_β) > 10.0
                CpelNano.print_log("MSE issue: $(rs.ϕhat)")
                continue
            end

            # Record
            if length(rs.ϕhat) > 0
                nse_mat[i,j] = mse(α, β, nse_α, nse_β)
            end
            if length(rs_bin.ϕhat) > 0
                no_nse_mat[i,j] = mse(α, β, no_nse_α, no_nse_β)
            end

        end

    end

    # Return
    return nse_mat, no_nse_mat

end
function bxplt(nse_mat, no_nse_mat, σs, reps)

    # Cat matrices
    nse_mat_vcat = vcat(nse_mat'...)
    nse_sig_vcat = vcat([fill("$(σ)", reps) for σ in σs]...)
    nse_grp_vcat = fill("M1", length(nse_sig_vcat))
    no_nse_mat_vcat = vcat(no_nse_mat'...)
    no_nse_sig_vcat = vcat([fill("$(σ)", reps) for σ in σs]...)
    no_nse_grp_vcat = fill("M2", length(no_nse_sig_vcat))
    data_vcat = vcat(nse_mat_vcat, no_nse_mat_vcat)
    sig_vcat = vcat(nse_sig_vcat, nse_sig_vcat)
    grp_vcat = vcat(nse_grp_vcat, no_nse_grp_vcat)

    # Remove NaNs
    non_nan_ind = .!isnan.(data_vcat)
    data_vcat = data_vcat[non_nan_ind]
    sig_vcat = sig_vcat[non_nan_ind]
    grp_vcat = grp_vcat[non_nan_ind]

    # Grouped boxplots
    p = groupedboxplot(sig_vcat,data_vcat,group=grp_vcat,outliers=false,
        xlabel="\\sigma",ylabel="MSE",size=(700, 600))

    # Return plot
    return p

end
#####################################################################################################
# Calls
#####################################################################################################

# AAAAACGAAAAA
# AAAAAC => log p(e1|AAAAAC) / AAAAAM => log p(e1|AAAAAM)
# AAAACG => log p(e2|AAAACG) / AAAAMG => log p(e2|AAAAMG)
# ...

# Parameters from pore model (AACGAA - AAMGAA)
μ_m = 101.93
μ_u = 103.13

# Plot current level distribution
# σs = [0.2,0.4,0.6,0.8]
σs = collect(0.1:0.1:0.8)
ps = plot_dists(μ_m, μ_u, σs)
plot(ps...,layout=(4, 2),size=(800, 1000))

# Run simulations
seed = 123
reps = 1000
nse_mat, no_nse_mat = sim(μ_m, μ_u, σs, seed, reps)

# Boxplot
p = bxplt(nse_mat, no_nse_mat, σs, reps);
plot(p)

# sum(nse_mat[1,:] .<=  no_nse_mat[1,:])
# sum(nse_mat[2,:] .<=  no_nse_mat[2,:])
# sum(nse_mat[3,:] .<=  no_nse_mat[3,:])
# sum(nse_mat[4,:] .<=  no_nse_mat[4,:])


text:
#
# CAVEAT: src/all_plots.jl generates all plots and save them into the scr/figures/ dir.
#
using GnuplotScripting
using Documenter

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

makedocs(;
    modules=[GnuplotScripting],
    authors="Vincent Picaud",
    repo="https://github.com/vincent-picaud/GnuplotScripting.jl/blob/{commit}{path}#{line}",
    sitename="GnuplotScripting.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://vincent-picaud.github.io/GnuplotScripting.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "New things to learn" => "new_things_to_learn.md",
        "Beside `free_form()`" => "beside_free_form.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/vincent-picaud/GnuplotScripting.jl",
    devbranch="main",
)


text:
export Gtk,
    GtkButton,
    GtkScale,
    GtkLabel,
    GtkEntry,
    GtkSpinButton,
    GtkSpinner,
    GtkTextView,
    GtkStatusbar,
    GtkSwitch,
    GtkLinkButton,
    GtkFontButton,
    GtkExpander,
    GtkAppChooser,
    GtkColorButton,
    GtkProgressBar,
    GtkFileChooser,
    GtkCheckButton,
    GtkToggleButton,
    GtkComboBoxText,
    GtkVolumeButton,
    GtkCanvas,
    GtkBox,
    GtkFrame,
    GtkPaned,
    GtkWindow,
    GtkMenuBar,
    GtkToolbar,
    GtkTreeView,
    GtkNotebook,
    GtkScrolledWindow,
    GtkWidget,
    GdkRGBA,
    GtkAlign,
    set_gtk_property!,
    get_gtk_property,
    waitforsignal,
    destroy,
    signal_connect,
    signal_handler_disconnect,
    bytestring, CssProvider,
    GAccessor, GConstants,
    show, showall,
    width, height

    #================ gtk.jl =================#
export set!,
    getprop,
    onevent,
    offevent,
    @inline_widget,
    @container_widget,
    value,
    setvalue!,
    add!,
    #================ widgets.jl =================#
    Scale,
    Slider,
    SpinButton,
    Button,
    @button_str,
    ColorButton,
    @colorbtn_str,
    Label,
    @label_str,
    Entry,
    @entry_str,
    Spinner,
    TextView,
    Canvas,
    StatusBar,
    Switch,
    LinkButton,
    FontButton,
    AppChooser,
    ProgressBar,
    FileChooser,
    CheckBox,
    ToggleButton,
    TextList,
    VolumeButton,
    Box,
    Frame,
    Paned,
    Window,
    MenuBar,
    Toolbar,
    TreeView,
    Notebook,
    ScrolledWindow,
    Expander,
    #================= utils.jl =================#
    Maybe,
    set!,
    align,
    keyboard,
    @margin,
    @valign,
    @halign,
    @align,
    @spacing,
    @expand,
    @hexpand,
    @vexpand,
    @homogeneous,
    @on,
    #================= math.jl =================#
    random,
    point,
    mapr,
    ⟶,
    #================= vector.jl =================#
    Vec,
    VECT_LETTER_INDICIES,
    set_vector_lindex!,
    extend!,
    marg,
    mag2,
    mag²,
    mag,
    rot,
    limit,
    randv,
    withangle,
    withmag,
    #================= app.jl =================#
    AlexyaApp,
    emit_event,
    add_listener,
    CURRENT_ALEXYA_APP,
    DEFAULT_LISTENERS,
    current_app,
    addwidget,
    @create,
    init, @init,
    use, @use,
    @window, @canvas,
    @width, @height,
    @framecount,
    @framerate,
    start, @canvasfocus,
    framerate, noLoop,
    #================= layout.jl =================#
    canvasonly,
    aside,
    uselayout,
    @layout,
    #================= grid.jl =================#
    span,
    →,
    ↓,
    Grid,
    #================= svg.jl =================#
    SVGCmd,
    SVGPath,
    getpoints,
    #================= sprites.jl =================#
    SpriteImage,
    loadsprite,
    drawsprite,
    #================ controls.jl =================#
    AbstractOption,
    Option,
    format_option_name,
    valuetype,
    create_option,
    Controls,
    NoLabel,
    @options,
    @nolabel


text:
Network = NetLib.ResNet

# netparams = NetLib.SimpleNetHP(
#   width=45,
#   depth_common=6,
#   use_batch_norm=true,
#   batch_norm_momentum=1.)

netparams = NetLib.ResNetHP(
  num_filters=128,
  num_blocks=5,
  conv_kernel_size=(3, 3),
  num_policy_head_filters=32,
  num_value_head_filters=32,
  batch_norm_momentum=0.1)

self_play = SelfPlayParams(
  sim=SimParams(
    num_games=5000,
    num_workers=128,
    use_gpu=true,
    reset_every=4,
    flip_probability=0.,
    alternate_colors=false),
  mcts=MctsParams(
    num_iters_per_turn=600,
    cpuct=1.0,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.2,
    dirichlet_noise_α=1.0))

arena = ArenaParams(
  sim=SimParams(
    num_games=100,
    num_workers=100,
    use_gpu=true,
    reset_every=1,
    flip_probability=0.,
    alternate_colors=false),
  mcts = MctsParams(
    self_play.mcts,
    dirichlet_noise_ϵ=0.1),
  update_threshold=0.05)

learning = LearningParams(
  use_gpu=true,
  samples_weighing_policy=LOG_WEIGHT,
  l2_regularization=1e-4,
  optimiser=Adam(lr=2e-3),
  batch_size=1024,
  loss_computation_batch_size=1024,
  nonvalidity_penalty=1.,
  min_checkpoints_per_epoch=1,
  max_batches_per_checkpoint=2000,
  num_checkpoints=1)

params = Params(
  arena=arena,
  self_play=self_play,
  learning=learning,
  num_iters=10,
  memory_analysis=MemAnalysisParams(
    num_game_stages=3),
  ternary_rewards=true,
  use_symmetries=false,
  mem_buffer_size=PLSchedule(
  [      0,        15],
  [400_000, 1_000_000]))

#####
##### Evaluation benchmark
#####

mcts_baseline =
  Benchmark.MctsRollouts(
    MctsParams(
      arena.mcts,
      num_iters_per_turn=1000,
      cpuct=1.))

minmax_baseline = Benchmark.MinMaxTS(
  depth=5,
  τ=0.2,
  amplify_rewards=true)

alphazero_player = Benchmark.Full(arena.mcts)

network_player = Benchmark.NetworkOnly(τ=0.5)

benchmark_sim = SimParams(
  arena.sim;
  num_games=128,
  num_workers=128)

benchmark = [
    Benchmark.Duel(alphazero_player, mcts_baseline,   benchmark_sim),
    Benchmark.Duel(alphazero_player, minmax_baseline, benchmark_sim),
    Benchmark.Duel(network_player,   mcts_baseline,   benchmark_sim),
    Benchmark.Duel(network_player,   minmax_baseline, benchmark_sim)
]

experiment = Experiment(
  "fanorona", GameSpec(), params, Network, netparams, benchmark)


text:
"""Q-learning Functions"""

module qlearning

    using JSON3

    export simulate_game

    function get_actions(game::Main.init.model, s::Array{Int8,1}, t::Int64)::Array{Int8,1}
        """Get actions"""
        a = Int8.(zeros(game.n));
        pr_explore::Float64 = exp(- t * game.beta);
        e = pr_explore .> rand(1,3);
        for n = 1:game.n
            if e[n]
                a[n] = rand(1:game.k);
            else
                a[n] = argmax(game.Q[n, s[1], s[2], :]);
            end
        end
        return a
    end

    function update_Q(game::Main.init.model, profits::Array{Float32,1}, s::Array{Int8,1}, s1::Array{Int8,1}, a::Array{Int8,1}, stable::Int64)::Tuple{Array{Float32,4},Int64}
        """Update Q function"""
        for n=1:game.n
            old_value = game.Q[n, s[1], s[2], a[n]];
            max_q1 = max(game.Q[n, s1[1], s1[2], :]...);
            new_value = profits[n] + game.delta * max_q1;
            old_argmax = argmax(game.Q[n, s[1], s[2], :])
            game.Q[n, s[1], s[2], a[n]] = (1 - game.alpha) * old_value + game.alpha * new_value;
            # Check stability
            new_argmax = argmax(game.Q[n, s[1], s[2], :])
            same_argmax = Int8(old_argmax == new_argmax);
            stable = (stable + same_argmax) * same_argmax
        end
        return game.Q, stable
    end

    function check_convergence(game::Main.init.model, t::Int64, stable::Int64)::Bool
        """Check if game converged"""
        if rem(t, game.tstable)==0
            print("\nstable = ", stable)
        end
        if stable > game.tstable
            print("\nConverged!")
            return true;
        elseif t==game.tmax
            print("\nERROR! Not Converged!")
            return true;
        else
            return false;
        end
    end

    function export_game(game::Main.init.model)
        """Export game"""
        open("output/games/game1.json", "w") do io
            JSON3.write(io, game)
        end
    end

    # Simulate game
    function simulate_game(game::Main.init.model)::Main.init.model
        s = game.s0;
        stable = 0
        # Iterate until convergence
        for t=1:game.tmax
            a = get_actions(game, s, t)
            profits = game.PI[a[1], a[2], :]
            s1 = a
            game.Q, stable = update_Q(game, profits, s, s1, a, stable)
            s = s1;
            if check_convergence(game, t, stable)
                break
            end
        end
        return game
    end

end


text:
using LightXML
xdoc = parse_file("data/books.xml");
xtop = root(xdoc);
println(name(xtop));

for c in child_nodes(xtop)
    if is_elementnode(c)
        e = XMLElement(c)
        t = find_element(e, "title")
        title = content(t)
        genre = attribute(t, "genre")
        @printf "%20s :%s\n" title genre
    end
end

for c in child_nodes(xtop)
   if is_elementnode(c)
      e = XMLElement(c)
      t = find_element(e, "title")
      genre = attribute(t, "genre")
      if genre == "Computing"
         a = find_element(e,"author")
         p = find_element(e,"price")
         curr = attribute(p, "currency")
         ds = find_element(e,"publish_date")
         desc = find_element(e,"description")
         println("Title:      ", content(t))
         println("Author:     " ,content(a))
         println("Date:       " ,ds)
         println("Price:      " ,p ," (", curr, ")")
         println(content(desc),"\n");
      end
   end
end

max_stars_repo_path: Module 3/Chapter05/light-xml.jl

max_issues_repo_path: Module 3/Chapter05/light-xml.jl

max_forks_repo_path: Module 3/Chapter05/light-xml.jl


text:
using IsApprox
using IsApprox: isone, iszero, isreal, isinteger, ispossemidef, isposdef, isdiag
using Test
import LinearAlgebra

@testset "isone, iszero" begin
    @test isone(1.0)
    @test ! isone(0.0)
    @test iszero(0.0)
    @test ! iszero(1.0)
    @test isone(1.0, Equal())
    @test ! isone(1.0 + 1e-10, Equal())
    @test isone(1.0 + 1e-10, Approx())
    @test isone(1.0 + 1e-7, Approx(atol=1e-6))
    @test ! isone(1.0 + 1e-7, Approx())

    # maximum element is 9.9e-11
    # norm(m) == 1.8892148837710255e-10
    m = [9.914830625828892e-11 4.473996666408231e-11 3.066184129948095e-11;
         9.914231937308404e-12 7.07270051593276e-11 6.487445728681865e-11;
         9.587109787642197e-11 6.520784660158512e-11 1.2903560723094866e-11]

    @test iszero(m, Approx(atol=2e-10)) # in norm
    @test ! iszero(m, Approx(atol=1e-10))
    @test iszero(m, EachApprox(atol=2e-10)) # element-wise
    @test iszero(m, EachApprox(atol=1e-10))
    @test ! iszero(m) # exact
    @test ! iszero(m, Equal()) # exact
end

@testset "isreal, isinteger" begin
    @test isreal(1)
    @test isreal(1.0)
    @test isreal(1.0, Approx())
    @test ! isreal(1.0 + 1e-10im)
    @test ! isreal(1.0 + 1e-10im, Approx())
    @test isreal(1.0 + 1e-10im, Approx(atol=1e-9))
    @test isinteger(1)
    @test isinteger(1, Approx())
end

@testset "isdiag" begin
    m0 = [1 0; 0 1]
    m = [1.0 6.108385298833888e-20; 7.691926633708195e-20 1.0]
    @test isdiag(m0)
    @test ! isdiag(m)
    @test isdiag(m; approx=Approx(atol=1e-10))
end

@testset "isposdef" begin
    n = 3
    m = one(rand(n, n))
    m2 = m + rand(n, n) * 1e-10
    @test ispossemidef(m)
    @test ! ispossemidef(m2)
    @test ispossemidef(m2, Approx(atol=1e-9))
    @test ispossemidef(m2, EachApprox(atol=1e-9))
    @test ispossemidef(0)
    @test ispossemidef(1)
    @test ! ispossemidef(-1)
    @test isposdef(1 + 0im)
    @test ispossemidef(-1e-5, Approx(atol=1e-3))

    @test isposdef(1)
    @test ! isposdef(0)
    @test ! isposdef(-1e5)
    @test isposdef(-1e-10, Approx(atol=1e-8))
    @test isposdef(1 + 1e-10im, Approx(atol=1e-8))
    @test ! isposdef(1 + 1e-10im, Equal())

    m3 = [1 0; 0 0]
    @test ispossemidef(m3)
    @test ! isposdef(m3)
    @test isposdef(m3, Approx())
end

@testset "isunitary isinvolution" begin
    s0 = [1 0; 0 1]
    s1 = [0 1; 1 0]
    s2 = [0 -im; im 0]
    s3 = [1 0; 0 -1]

    @test isunitary(s3)
    @test isunitary(s2)
    @test isinvolution(s3)
    @test isinvolution(s2)

    # merr = rand(0:100, 8,8) * 1e-12
    errmat = [5.3e-11 9.4e-11 3.5e-11 6.6e-11 8.6e-11 4.9e-11 5.2e-11 1.6e-11;
              2.8e-11 3.9e-11 4.0e-11 6.0e-12 3.1e-11 2.9e-11 9.7e-11 7.1e-11;
              8.0e-11 2.5e-11 8.5e-11 7.2e-11 5.7e-11 4.8e-11 5.8e-11 4.1e-11;
              5.5e-11 9.8e-11 4.7e-11 8.7e-11 3.6e-11 9.0e-11 4.3e-11 8.0e-11;
              7.0e-11 1.5e-11 2.1e-11 9.8e-11 9.7e-11 2.1e-11 7.6e-11 4.0e-11;
              4.8e-11 3.9e-11 9.8e-11 0.0 4.0e-11 4.0e-12 4.0e-11 5.3e-11;
              6.8e-11 8.3e-11 4.0e-12 1.5e-11 6.7e-11 3.4e-11 1.0e-11 9.0e-12;
              7.0e-11 5.2e-11 2.2e-11 7.2e-11 9.0e-12 5.4e-11 4.9e-11 5.3e-11]

    m = kron(s1, s2, s3)
    merr = m + errmat

    @test isunitary(m)
    @test isunitary(m, Equal())
    @test ! isunitary(merr)
    @test ! isunitary(merr, Equal())
    @test isunitary(merr, Approx())
    @test ! isunitary(merr, EachApprox())
    @test ! isunitary(merr, EachApprox(atol=1e-10))
    @test isunitary(merr, EachApprox(atol=5e-10))
    @test ! isunitary(merr, Approx(atol=1e-10))

    m1 = [1 0; 1 -1]
    _m1err = (rand(2, 2) .- 0.5) * 1e-10
    m1err = m1 + _m1err

    @test isinvolution(m1)
    @test ! isunitary(m1)
    @test isinvolution(m1, Equal())
    @test ! isinvolution(m1err)
    @test isinvolution(m1err, EachApprox(atol=1e-9))
    @test isinvolution(m1err, Approx(atol=1e-9))
    @test ! isunitary(m1err, Approx(atol=1e-9))
    @test ! isinvolution(m1err, Approx(atol=1e-12))

    @test isunitary(1)
    @test ! isunitary(2)
    @test isunitary((1e-10 + exp(im * 2.5)) * LinearAlgebra.I, Approx())

end

@testset "UpToPhase" begin
    @test isapprox(UpToPhase(), 2, 2)
    @test isapprox(UpToPhase(), 2.1, 2.1 + 1e-10)
    @test ! isapprox(UpToPhase(atol=1e-12), 2.1, 2.1 + 1e-10)
    @test isapprox(UpToPhase(), 2.1, 2.1 * cis(2*pi*3.1))
    m = rand(5, 5)
    @test isapprox(UpToPhase(), m, m .* cis(2*pi*1.3))
    @test isapprox(UpToPhase(), m, m .* (1+1e-10)*cis(2*pi*1.3))
    @test ! isapprox(UpToPhase(), m, m .* (1+1e-7)*cis(2*pi*1.3))
end


text:
# To run the code, use the command: include(filename)
# Do this in the julia prompt
# Example: include("hello_world.jl")

example_string = "Hello, world"
println(example_string)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
module SSLTestCHELSA
using SimpleSDMLayers
using Test

lc1 = SimpleSDMPredictor(CHELSA, BioClim, 1; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test typeof(lc1) <: SimpleSDMPredictor

lc1 = SimpleSDMPredictor(CHELSA, BioClim, 1:2; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test eltype(lc1) <: SimpleSDMPredictor

lc1 = SimpleSDMPredictor(CHELSA, BioClim, [1,2,3]; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test eltype(lc1) <: SimpleSDMPredictor

lc1 = SimpleSDMPredictor(CHELSA, BioClim, CanESM2, RCP26, 1; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test typeof(lc1) <: SimpleSDMPredictor

lc1 = SimpleSDMPredictor(CHELSA, BioClim, CanESM2, RCP26, 1:2; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test eltype(lc1) <: SimpleSDMPredictor

lc1 = SimpleSDMPredictor(CHELSA, BioClim, CanESM2, RCP26, [1,2,3]; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test eltype(lc1) <: SimpleSDMPredictor

end


text:
# Configuration for running simple neutral copy error simulation
#=
Recommended command line to run:
>  julia neutral.jl examples/sn_example1
=#
export simtype
@everywhere simtype = 3
@everywhere const N = 1        # population size
const mutstddev = 0.04
const ngens = 1000
#const initial_value = exp(5.0)
const initial_value = 1.0
const num_trials=100000
const record_interval = 50
const use_population = true
const log_error = true


text:
using HPAT

HPAT.setBlockSize(2)

@acc hpat function mult(file1,file2)
    @partitioned(M,HPAT_2D);
    M = DataSource(Matrix{Float64},HDF5,"/M", file1)
    DataSink(M,HDF5,"/y", file2)
    return 0
end

using MPI
rank = MPI.Comm_rank(MPI.COMM_WORLD)
using HDF5
if rank==0 h5write("dtest.hdf5","/M", reshape(convert(Vector{Float64},collect(1:16)),4,4)) end

mult("dtest.hdf5","testout.hdf5")
if rank==0
    A = h5read("testout.hdf5","/y")
    @test_approx_eq A reshape(convert(Vector{Float64},collect(1:16)),4,4)
    rm("dtest.hdf5")
    rm("testout.hdf5")
end


text:
using SpectralLibraries
using Test

@testset "SpectralLibraries.jl" begin
    # Write your own tests here.
end


text:
# ---------------------------------------------------------------------------
# Writers
# ---------------------------------------------------------------------------
using Printf

# Note: the client will read this data for display.
function samples_writer(samples::SpanArray, model::ModelData, f::IOStream)
    synapses = Model.synapses(model)
    span_time = Model.span_time(model)

    # Now write each stream/synpase-input
    # Format:
    # id 1010001011010...::
    for id in 1:synapses
        # Write stream id
        print(f, @sprintf("%03d ", id))

        # write all spikes (1) and non-spikes (0)
        # print("($id) ")
        for t in 1:span_time
            # print(samples.poi_samples[id, t])
            print(f, samples[id, t])
        end

        # Terminate synapse stream with "::" marker
        # println("::")
        println(f, "::")
    end
end

function cell_samples_writer(samples::Array{UInt8,1}, model::ModelData, f::IOStream)
    span_time = Model.span_time(model)

    # Now write each stream/synpase-input
    # Format:
    # id 1010001011010...::

    print(f, @sprintf("%03d ", 1))

    # write all spikes (1) and non-spikes (0)
    for t in 1:span_time
        print(f, samples[t])
    end

    # Terminate synapse stream with "::" marker
    println(f, "::")
end

function float_samples_writer(samples::Array{Float64,1}, model::ModelData, f::IOStream)
    span_time = Model.span_time(model)

    # Format:
    # float
    # float
    # ...

    for t in 1:span_time
        println(f, samples[t])
    end
end

function synapse_writer(syn_samples::SynapticSamples, model::ModelData, f::IOStream)
    span_time = Model.span_time(model)

    # Format:
    # id float,float,float...
    # id float,float,float...
    # Nth float,float,float...

    synapses = Model.synapses(model)

    for syn_id in 1:synapses
        syn_data = syn_samples.data[syn_id]

        print(f, @sprintf("%03d ", syn_id))

        for t in 1:span_time - 1
            print(f, syn_data.samples[t], " ")
        end

        println(f, syn_data.samples[span_time])
    end

end

# Write the current samples to disk, then reset.
function write_samples_out(samples::Array, file::String, model::ModelData, writer)
    open(file, "w") do f
        writer(samples, model, f)
    end
end

function write_samples_out(data::SampleData, file::String, model::ModelData, writer)
    open(file, "w") do f
        writer(data.samples, model, f)
    end
end

function write_samples_out(samples::SynapticSamples, file::String, model::ModelData, writer)
    open(file, "w") do f
        writer(samples, model, f)
    end
end


text:
@info "adaptive_loss_logging_tests"
using DiffEqFlux
using ModelingToolkit
using Test, NeuralPDE
using Optimization, OptimizationFlux
import ModelingToolkit: Interval, infimum, supremum
using Random
#using Plots
@info "Starting Soon!"

nonadaptive_loss = NeuralPDE.NonAdaptiveLoss(pde_loss_weights=1, bc_loss_weights=1)
gradnormadaptive_loss = NeuralPDE.GradientScaleAdaptiveLoss(100, pde_loss_weights=1e3, bc_loss_weights=1)
adaptive_loss = NeuralPDE.MiniMaxAdaptiveLoss(100; pde_loss_weights=1, bc_loss_weights=1)
adaptive_losses = [nonadaptive_loss, gradnormadaptive_loss,adaptive_loss]
maxiters=800
seed=60

## 2D Poisson equation
function test_2d_poisson_equation_adaptive_loss(adaptive_loss, run, outdir, haslogger; seed=60, maxiters=800)
    logdir = joinpath(outdir, string(run))
    if haslogger
        logger = TBLogger(logdir)
    else
        logger = nothing
    end
    Random.seed!(seed)
    hid = 40
    chain_ = FastChain(FastDense(2,hid,Flux.σ),FastDense(hid,hid,Flux.σ),FastDense(hid,1))
    strategy_ =  NeuralPDE.StochasticTraining(256)
    @info "adaptive reweighting test logdir: $(logdir), maxiters: $(maxiters), 2D Poisson equation, adaptive_loss: $(nameof(typeof(adaptive_loss))) "
    @parameters x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2

    # 2D PDE
    eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

    # Initial and boundary conditions
    bcs = [u(0,y) ~ 0.0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
           u(x,0) ~ 0.0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
    # Space and time domains
    domains = [x ∈ Interval(0.0,1.0),
               y ∈ Interval(0.0,1.0)]

    initθ = Float64.(DiffEqFlux.initial_params(chain_))
    iteration = [0]
    discretization = NeuralPDE.PhysicsInformedNN(chain_,
                                                 strategy_;
                                                 init_params = initθ,
                                                 adaptive_loss = adaptive_loss,
                                                 logger = logger,
                                                 iteration=iteration)

    @named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
    prob = NeuralPDE.discretize(pde_system,discretization)
    phi = discretization.phi
    sym_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)

    xs,ys = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
    analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)
    u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))

    callback = function (p,l)
        iteration[1] += 1
        if iteration[1] % 100 == 0
            @info "Current loss is: $l, iteration is $(iteration[1])"
        end
        if haslogger
            log_value(logger, "outer_error/loss", l, step=iteration[1])
            if iteration[1] % 30 == 0
                u_predict = reshape([first(phi([x,y],p)) for x in xs for y in ys],(length(xs),length(ys)))
                diff_u = abs.(u_predict .- u_real)
                total_diff = sum(diff_u)
                log_value(logger, "outer_error/total_diff", total_diff, step=iteration[1])
                total_u = sum(abs.(u_real))
                total_diff_rel = total_diff / total_u
                log_value(logger, "outer_error/total_diff_rel", total_diff_rel, step=iteration[1])
                total_diff_sq = sum(diff_u .^ 2)
                log_value(logger, "outer_error/total_diff_sq", total_diff_sq, step=iteration[1])
            end
        end
        return false
    end
    res = Optimization.solve(prob, ADAM(0.03); maxiters=maxiters, callback=callback)

    u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
    diff_u = abs.(u_predict .- u_real)
    total_diff = sum(diff_u)
    total_u = sum(abs.(u_real))
    total_diff_rel = total_diff / total_u

    #p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
    #p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
    #p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
    #(plot=plot(p1,p2,p3), error=total_diff, total_diff_rel=total_diff_rel)
    (error=total_diff, total_diff_rel=total_diff_rel)
end

possible_logger_dir = mktempdir()
if ENV["LOG_SETTING"] == "NoImport"
    haslogger = false
    expected_log_folders = 0
elseif ENV["LOG_SETTING"] == "ImportNoUse"
    using NeuralPDELogging
    haslogger = false
    expected_log_folders = 0
elseif ENV["LOG_SETTING"] == "ImportUse"
    using NeuralPDELogging
    using TensorBoardLogger
    haslogger = true
    expected_log_folders = 3
end

@info "has logger: $(haslogger), expected log folders: $(expected_log_folders)"

test_2d_poisson_equation_adaptive_loss_run_seediters(adaptive_loss, run) = test_2d_poisson_equation_adaptive_loss(adaptive_loss, run, possible_logger_dir, haslogger; seed=seed, maxiters=maxiters)
error_results = map(test_2d_poisson_equation_adaptive_loss_run_seediters, adaptive_losses, 1:length(adaptive_losses))

@test length(readdir(possible_logger_dir)) == expected_log_folders
if expected_log_folders > 0
    @info "dirs at $(possible_logger_dir): $(string(readdir(possible_logger_dir)))"
    for logdir in readdir(possible_logger_dir)
        @test length(readdir(joinpath(possible_logger_dir, logdir))) > 0
    end
end